src.it.unimi.dsi.stat.Ziggurat Maven / Gradle / Ivy

Go to download

Show more of this group Show more artifacts with this name
Show all versions of dsiutils Show documentation

The DSI utilities are a mishmash of classes accumulated during the last twenty years in projects developed at the DSI (Dipartimento di Scienze dell'Informazione, i.e., Information Sciences Department), now DI (Dipartimento di Informatica, i.e., Informatics Department), of the Universita` degli Studi di Milano.

There is a newer version: 2.7.3

Show newest version

package it.unimi.dsi.stat;

/*
 * DSI utilities
 *
 * Copyright (C) 2012-2020 Sebastiano Vigna
 *
 *  This library is free software; you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License as published by the Free
 *  Software Foundation; either version 3 of the License, or (at your option)
 *  any later version.
 *
 *  This library is distributed in the hope that it will be useful, but
 *  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 *  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License
 *  for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with this program; if not, see .
 *
 */

import it.unimi.dsi.Util;
import it.unimi.dsi.bits.Fast;
import it.unimi.dsi.util.SplitMix64Random;

import java.util.Random;

public class Ziggurat {
	private final Random random;
	private static double[] x = new double[257];
	private static double[] x53 = new double[257];
	private static double[] y = new double[257];
	private static double[] h = new double[257];
	private static double[] h53 = new double[257];
	private static long[] threshold = new long[257];

	static {
		final double area = 3.949659822581572e-3;

		/* Warm-up. This is the point defining the bottom (tail) part.
		 * Note that 7.697117470131487 is the target for x[1]
		 * (from the original code). */
		x[0] = 8.69711747013488555703;
		y[0] = 0;
		h[0] = area / x[0]; // Never used, but necessary to compute the remaining data.

		for (int i = 1; i < 256; i++) {
			y[i] = y[i - 1] + h[i - 1];
			x[i] = -Math.log(y[i]);
			h[i] = area / x[i];
		}

		for (int i = 1; i < 256; i++) {
			threshold[i] = Math.round((1L << 53) * x[i + 1] / x[i]);
			// System.err.println(x[i] + " " + y[i] + " " + threshold[i - 1]);
		}

		for (int i = 0; i < 256; i++) {
			x53[i] = x[i] / (1L << 53);
			h53[i] = h[i] / (1L << 53);
		}

		y[256] = 1; // x[256] is 0.
	}

	public Ziggurat() {
		this.random = new SplitMix64Random();
	}

	public Ziggurat(final Random random) {
		this.random = random;
	}

	public double nextDouble() {
		for(;;) {
			final long r = random.nextLong();
			// We never extract 256 because it would fail anyway.
			final int block = (int)(r & 0xFF);
			// This is a discrete representation of a uniformly chosen random number in [0..x[block]).
			final long discreteCandidate = r >>> 11;
			if (discreteCandidate < threshold[block]) {
				return discreteCandidate * x53[block];
			}
			if (block == 0) {
				// Slow but very rare.
				return x[1] - Math.log(random.nextDouble());
			}
			final double candidate = discreteCandidate * x53[block];
			if (y[block] + (random.nextLong() >>> 11) * h53[block] < Math.exp(-candidate)) return candidate;
		}
	}

	public static void main(String arg[]) {
		Ziggurat ziggurat = new Ziggurat(new SplitMix64Random());
		final long time = - System.currentTimeMillis();
		final int n = Integer.parseInt(arg[0]);
		double sum = 0;
		int numExp = 100;
		final int exponentBits = 6;
		final int maxExponent = 2;
		final int minExponent = - (1 << exponentBits) + maxExponent + 1;
		for(int t = numExp; t-- != 0;) {
			int x = Integer.MAX_VALUE;
			for(int k = n; k-- != 0;) {
				x = Math.min(x, Math.max(Math.min(maxExponent, Fast.approximateLog2(ziggurat.nextDouble())), minExponent));
			}
			sum += Fast.pow2(x);
		}

		System.err.println(1000000.0 * (time + System.currentTimeMillis()) / (n * numExp) + " ns/gen");
		System.out.println(numExp / sum + " (" + Util.format(100 * (numExp / sum - n) / n) + "%)");
	}
}