import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.security.NoSuchAlgorithmException;
import java.security.SecureRandom;
import java.util.Random;

public class StatsBigIntegerRandom {
    public static void main(String[] args) throws NoSuchAlgorithmException {
        if (args.length != 3) {
            throw new RuntimeException
                 ("Usage: java StatsBigIntegerRandom bits binSize count");
        }

        final int bits = Integer.valueOf(args[0]);
        final int binSize = Integer.valueOf(args[1]);
        assert bits % 32 == 0 : bits + " is not a multiple of 32";
        assert bits % binSize == 0 : bits + " is not a multiple of " + binSize;
        final int count = Integer.valueOf(args[2]);
        final BigInteger mask =
            BigInteger.ONE.shiftLeft(binSize).subtract(BigInteger.ONE);

        System.out.println("--- Parameters ---");
        System.out.format("Max bits: %d%n", bits);
        System.out.format("Bin size: %d%n", binSize);
        System.out.format("Number of samples: %d%n%n", count);

        Random rnd = SecureRandom.getInstance("SHA1PRNG");

        final int numBins = bits/binSize + 1;
        final int[] hist = new int[numBins];
        final int maxBin = numBins - 1;

        BigInteger sum = BigInteger.ZERO;
        BigInteger sumOfSquares = BigInteger.ZERO;

        for (int i = 0; i < count; ++i) {
            BigInteger big = new BigInteger(bits, rnd);

            sum = sum.add(big);
            sumOfSquares = sumOfSquares.add(big.multiply(big));

            int bitLength = big.bitLength();
            if (bitLength == 0) {
                hist[maxBin]++;
            } else {
                int bins = (bitLength + binSize - 1)/binSize;
                for (int j = maxBin - bins; j < maxBin; j++) {
                    BigInteger n =
                        big.shiftRight(binSize*(maxBin - j - 1)).and(mask);
                    if (n.longValueExact() > 0) {
                        hist[j]++;
                        break;
                    }
                }
            }
        }

        BigDecimal bigCount = BigDecimal.valueOf(count);
        BigDecimal mean = new BigDecimal(sum).divide(bigCount);
        BigDecimal squareMean = new BigDecimal(sumOfSquares).divide(bigCount);
        BigDecimal variance = squareMean.subtract(mean.multiply(mean));
        MathContext mc = new MathContext(3);
        BigDecimal stdev = variance.sqrt(mc);
        BigDecimal cv = stdev.divide(mean, mc);

        System.out.println("--- Statistics ---");
        System.out.format("Bit length of mean: %d%n",
            mean.toBigInteger().bitLength());
        System.out.format("Coefficient of variation: %s%n", cv);
        System.out.println();

        System.out.println("--- Histogram ---");
        System.out.format("Bin #: Number of values for which this is the largest non-zero %d-bit bin.%n%n", binSize);

        long hSum = 0L;
        for (int k = 0; k < numBins; k++) {
            int h = hist[k];
            hSum += h;
            System.out.format("%5d: %d%n", k, h);
        }

        assert hSum == count : "expected " + count + ", got " + hSum;
    }
}