001package squidpony.squidmath;
002
003import squidpony.StringKit;
004
005import java.io.Serializable;
006
007/**
008 * A very-high-quality StatefulRandomness that is meant to be reasonably fast, but also to be robust against frequent
009 * state changes, and is built around a strong determine() method. It passes 32TB of PractRand with no anomalies, and
010 * calling {@link #determine(long)} on two different sequences of inputs (one that added 3 each time, and one that added
011 * 5 each time) showed no failures on 32TB of data produced by XORing calls to determine() on both sequences. PulleyRNG
012 * is one-dimensionally equidistributed across all 64-bit outputs, has 64 bits of state, natively outputs 64 bits at a
013 * time, can have its state set and skipped over as a {@link StatefulRandomness} and {@link SkippingRandomness}. It has
014 * a 1-to-1 correspondence between inputs and outputs for {@link #determine(long)}, and you can get the input to
015 * determine() that produced a particular long output by passing that output to {@link #inverseNextLong(long)}. It is
016 * largely the work of Pelle Evensen, who discovered that where a unary hash (called a determine() method here) can
017 * start with the XOR of the input and two rotations of that input, and that sometimes acts as a better randomization
018 * procedure than multiplying by a large constant (which is what {@link LightRNG#determine(long)},
019 * {@link LinnormRNG#determine(long)}, and even {@link ThrustAltRNG#determine(long)} do). Evensen also crunched the
020 * numbers to figure out that {@code n ^ n >>> A ^ n >>> B} is a bijection for all distinct non-zero values for A and B,
021 * though this wasn't used in his unary hash rrxmrrxmsx_0.
022 * <br>
023 * The algorithm for {@link #determine(long)} looks like this ({@link #nextLong()} just calls determine() on a counter):
024 * <ol>
025 *     <li>XOR the input with two different bitwise rotations: {@code n ^ (n << 41 | n >>> 23) ^ (n << 17 | n >>> 47)}</li>
026 *     <li>Multiply by a large constant, {@code 0x369DEA0F31A53F85L}, and store it in n</li>
027 *     <li>XOR n with two different unsigned bitwise right shifts: {@code n ^ n >>> 25 ^ n >>> 37}</li>
028 *     <li>Multiply by a large constant, {@code 0xDB4F0B9175AE2165L}, and store it in n</li>
029 *     <li>XOR n with n right-shifted by 28, and return</li>
030 * </ol>
031 * This is the result of some simplifications on PelicanRNG (present here as {@link DiverRNG#randomize(long)}), which is
032 * ridiculously strong (passing the full battery of tests that rrxmrrxmsx_0 only narrowly failed) but not especially
033 * fast. PulleyRNG is an effort to speed up PelicanRNG just a little, but without doing the extensive testing that
034 * ensure virtually any bit pattern given to PelicanRNG will produce pseudo-random outputs. PulleyRNG does well in tough
035 * tests. Other than the input stream correlation test mentioned earlier, this also passes tests if the inputs are
036 * incremented by what is normally one of the worst-case scenarios for other generators -- using an increment that is
037 * the multiplicative inverse (mod 2 to the 64 in this case) of one of the fixed constants in the generator. The first
038 * multiplication performed here is by {@code 0x369DEA0F31A53F85L}, and
039 * {@code 0xBE21F44C6018E14DL * 0x369DEA0F31A53F85L == 1L}, so testing determine() with inputs that change by
040 * 0xBE21F44C6018E14DL should stress the generator, but instead it does fine through 32TB, with only one "unusual"
041 * anomaly rather early on.
042 * <br>
043 * @author Pelle Evensen
044 * @author Tommy Ettinger
045 */
046public final class PulleyRNG implements StatefulRandomness, SkippingRandomness, Serializable {
047
048    private static final long serialVersionUID = 1L;
049
050    private long state; /* The state can be seeded with any value. */
051
052    /**
053     * Creates a new generator seeded using Math.random.
054     */
055    public PulleyRNG() {
056        this((long) ((Math.random() - 0.5) * 0x10000000000000L)
057                ^ (long) (((Math.random() - 0.5) * 2.0) * 0x8000000000000000L));
058    }
059
060    public PulleyRNG(final long seed) {
061        state = seed;
062    }
063
064    public PulleyRNG(final String seed) {
065        state = CrossHash.hash64(seed);
066    }
067
068    @Override
069    public final int next(int bits)
070    {
071        long z = state++;
072        z = (z ^ (z << 41 | z >>> 23) ^ (z << 17 | z >>> 47)) * 0x369DEA0F31A53F85L;
073        z = (z ^ z >>> 25 ^ z >>> 37) * 0xDB4F0B9175AE2165L;
074        return (int)((z ^ z >>> 28) >>> (64 - bits));
075    }
076
077    /**
078     * Can return any long, positive or negative, of any size permissible in a 64-bit signed integer.
079     *
080     * @return any long, all 64 bits are random
081     */
082    @Override
083    public final long nextLong() {
084        long z = state++;
085        z = (z ^ (z << 41 | z >>> 23) ^ (z << 17 | z >>> 47)) * 0x369DEA0F31A53F85L;
086        z = (z ^ z >>> 25 ^ z >>> 37) * 0xDB4F0B9175AE2165L;
087        return (z ^ z >>> 28);
088    }
089
090
091    /**
092     * Produces a copy of this RandomnessSource that, if next() and/or nextLong() are called on this object and the
093     * copy, both will generate the same sequence of random numbers from the point copy() was called. This just need to
094     * copy the state so it isn't shared, usually, and produce a new value with the same exact state.
095     *
096     * @return a copy of this RandomnessSource
097     */
098    @Override
099    public PulleyRNG copy() {
100        return new PulleyRNG(state);
101    }
102
103    @Override
104    public final long skip(final long advance) {
105        long z = state += advance;
106        z = (z ^ (z << 41 | z >>> 23) ^ (z << 17 | z >>> 47)) * 0x369DEA0F31A53F85L;
107        z = (z ^ z >>> 25 ^ z >>> 37) * 0xDB4F0B9175AE2165L;
108        return (z ^ z >>> 28);
109    }
110
111    /**
112     * Can return any int, positive or negative, of any size permissible in a 32-bit signed integer.
113     *
114     * @return any int, all 32 bits are random
115     */
116    public final int nextInt() {
117        long z = state++;
118        z = (z ^ (z << 41 | z >>> 23) ^ (z << 17 | z >>> 47)) * 0x369DEA0F31A53F85L;
119        z = (z ^ z >>> 25 ^ z >>> 37) * 0xDB4F0B9175AE2165L;
120        return (int)(z ^ z >>> 28);
121    }
122
123    /**
124     * Exclusive on the outer bound.  The inner bound is 0.
125     * The bound can be negative, which makes this produce either a negative int or 0.
126     *
127     * @param bound the upper bound; should be positive
128     * @return a random int between 0 (inclusive) and bound (exclusive)
129     */
130    public final int nextInt(final int bound) {
131        long z = state++;
132        z = (z ^ (z << 41 | z >>> 23) ^ (z << 17 | z >>> 47)) * 0x369DEA0F31A53F85L;
133        z = (z ^ z >>> 25 ^ z >>> 37) * 0xDB4F0B9175AE2165L;
134        return (int)((bound * ((z ^ z >>> 28) & 0xFFFFFFFFL)) >> 32);
135    }
136
137    /**
138     * Inclusive inner, exclusive outer.
139     *
140     * @param inner the inner bound, inclusive, can be positive or negative
141     * @param outer the outer bound, exclusive, can be positive or negative, usually greater than inner
142     * @return a random int between inner (inclusive) and outer (exclusive)
143     */
144    public final int nextInt(final int inner, final int outer) {
145        return inner + nextInt(outer - inner);
146    }
147
148    /**
149     * Exclusive on bound (which may be positive or negative), with an inner bound of 0.
150     * If bound is negative this returns a negative long; if bound is positive this returns a positive long. The bound
151     * can even be 0, which will cause this to return 0L every time.
152     * <br>
153     * Credit for this method goes to <a href="https://oroboro.com/large-random-in-range/">Rafael Baptista's blog</a>
154     * for the original idea, and the JDK10 Math class' usage of Hacker's Delight code for the current algorithm. 
155     * This method is drastically faster than the previous implementation when the bound varies often (roughly 4x
156     * faster, possibly more). It also always gets exactly one random long, so by default it advances the state as much
157     * as {@link #nextLong()}.
158     *
159     * @param bound the outer exclusive bound; can be positive or negative
160     * @return a random long between 0 (inclusive) and bound (exclusive)
161     */
162    public long nextLong(long bound) {
163        long rand = state++;
164        rand = (rand ^ (rand << 41 | rand >>> 23) ^ (rand << 17 | rand >>> 47)) * 0x369DEA0F31A53F85L;
165        rand = (rand ^ rand >>> 25 ^ rand >>> 37) * 0xDB4F0B9175AE2165L;
166        rand ^= rand >>> 28;
167        final long randLow = rand & 0xFFFFFFFFL;
168        final long boundLow = bound & 0xFFFFFFFFL;
169        rand >>>= 32;
170        bound >>= 32;
171        final long t = rand * boundLow + (randLow * boundLow >>> 32);
172        return rand * bound + (t >> 32) + (randLow * bound + (t & 0xFFFFFFFFL) >> 32);
173    }
174    /**
175     * Inclusive inner, exclusive outer; lower and upper can be positive or negative and there's no requirement for one
176     * to be greater than or less than the other.
177     *
178     * @param lower the lower bound, inclusive, can be positive or negative
179     * @param upper the upper bound, exclusive, can be positive or negative
180     * @return a random long that may be equal to lower and will otherwise be between lower and upper
181     */
182    public final long nextLong(final long lower, final long upper) {
183        return lower + nextLong(upper - lower);
184    }
185    
186    /**
187     * Gets a uniform random double in the range [0.0,1.0)
188     *
189     * @return a random double at least equal to 0.0 and less than 1.0
190     */
191    public final double nextDouble() {
192        long z = state++;
193        z = (z ^ (z << 41 | z >>> 23) ^ (z << 17 | z >>> 47)) * 0x369DEA0F31A53F85L;
194        z = (z ^ z >>> 25 ^ z >>> 37) * 0xDB4F0B9175AE2165L;
195        return ((z ^ z >>> 28) & 0x1FFFFFFFFFFFFFL) * 0x1p-53;
196
197    }
198
199    /**
200     * Gets a uniform random double in the range [0.0,outer) given a positive parameter outer. If outer
201     * is negative, it will be the (exclusive) lower bound and 0.0 will be the (inclusive) upper bound.
202     *
203     * @param outer the exclusive outer bound, can be negative
204     * @return a random double between 0.0 (inclusive) and outer (exclusive)
205     */
206    public final double nextDouble(final double outer) {
207        long z = state++;
208        z = (z ^ (z << 41 | z >>> 23) ^ (z << 17 | z >>> 47)) * 0x369DEA0F31A53F85L;
209        z = (z ^ z >>> 25 ^ z >>> 37) * 0xDB4F0B9175AE2165L;
210        return ((z ^ z >>> 28) & 0x1FFFFFFFFFFFFFL) * 0x1p-53 * outer;
211    }
212
213    /**
214     * Gets a uniform random float in the range [0.0,1.0)
215     *
216     * @return a random float at least equal to 0.0 and less than 1.0
217     */
218    public final float nextFloat() {
219        long z = state++;
220        z = (z ^ (z << 41 | z >>> 23) ^ (z << 17 | z >>> 47)) * 0x369DEA0F31A53F85L;
221        return ((z ^ z >>> 25 ^ z >>> 37) * 0xDB4F0B9175AE2165L >>> 40) * 0x1p-24f;
222    }
223
224    /**
225     * Gets a random value, true or false.
226     * Calls nextLong() once.
227     *
228     * @return a random true or false value.
229     */
230    public final boolean nextBoolean() {
231        long z = state++;
232        z = (z ^ (z << 41 | z >>> 23) ^ (z << 17 | z >>> 47)) * 0x369DEA0F31A53F85L;
233        return ((z ^ z >>> 25 ^ z >>> 37) * 0xDB4F0B9175AE2165L) < 0;
234    }
235
236    /**
237     * Given a byte array as a parameter, this will fill the array with random bytes (modifying it
238     * in-place). Calls nextLong() {@code Math.ceil(bytes.length / 8.0)} times.
239     *
240     * @param bytes a byte array that will have its contents overwritten with random bytes.
241     */
242    public final void nextBytes(final byte[] bytes) {
243        int i = bytes.length, n;
244        while (i != 0) {
245            n = Math.min(i, 8);
246            for (long bits = nextLong(); n-- != 0; bits >>>= 8) bytes[--i] = (byte) bits;
247        }
248    }
249
250    /**
251     * Sets the seed (also the current state) of this generator.
252     *
253     * @param seed the seed to use for this PulleyRNG, as if it was constructed with this seed.
254     */
255    @Override
256    public final void setState(final long seed) {
257        state = seed;
258    }
259
260    /**
261     * Gets the current state of this generator.
262     *
263     * @return the current seed of this PulleyRNG, changed once per call to nextLong()
264     */
265    @Override
266    public final long getState() {
267        return state;
268    }
269
270    /**
271     * Given the output of a call to {@link #nextLong()} as {@code out}, this finds the state of the PulleyRNG that
272     * produce that output. If you set the state of a PulleyRNG with {@link #setState(long)} to the result of this
273     * method and then call {@link #nextLong()} on it, you should get back {@code out}. This can also reverse
274     * {@link #determine(long)}; it uses the same algorithm as nextLong().
275     * <br>
276     * This isn't as fast as {@link #nextLong()}, but both run in constant time.
277     * <br>
278     * This will not necessarily work if out was produced by a generator other than a PulleyRNG, or if it was produced
279     * with the bounded {@link #nextLong(long)} method by any generator.
280     * @param out a long as produced by {@link #nextLong()}, without changes
281     * @return the state of the RNG that will produce the given long
282     */
283    public static long inverseNextLong(long out)
284    {
285        out ^= out >>> 28 ^ out >>> 56; // inverts xorshift by 28
286        out *= 0xF179F93568D4286DL; // inverse of 0xDB4F0B9175AE2165L
287        out ^= out >>> 25 ^ out >>> 50 ^ out >>> 37; // inverts paired xorshift by 25, 37
288        out *= 0xBE21F44C6018E14DL; // inverse of 0x369DEA0F31A53F85L
289        // follow the steps from http://marc-b-reynolds.github.io/math/2017/10/13/XorRotate.html
290        // this is the inverse of (0, 17, 41), working on 64-bit numbers.
291        out ^= (out << 17 | out >>> 47) ^ (out << 41 | out >>> 23);
292        out ^= (out << 34 | out >>> 30) ^ (out << 18 | out >>> 46);
293        out ^= (out << 4  | out >>> 60) ^ (out << 36 | out >>> 28);
294        return out;
295    }
296
297
298    @Override
299    public String toString() {
300        return "PulleyRNG with state 0x" + StringKit.hex(state) + 'L';
301    }
302
303    @Override
304    public boolean equals(Object o) {
305        if (this == o) return true;
306        if (o == null || getClass() != o.getClass()) return false;
307        return state == ((PulleyRNG) o).state;
308    }
309
310    @Override
311    public int hashCode() {
312        return (int) (state ^ state >>> 32);
313    }
314
315    /**
316     * Static randomizing method that takes its state as a parameter; state is expected to change between calls to this.
317     * It is recommended that you use {@code PulleyRNG.determine(++state)} or {@code PulleyRNG.determine(--state)} to
318     * produce a sequence of different numbers, but you can also use {@code PulleyRNG.determine(state += 12345L)} or
319     * any odd-number increment. All longs are accepted by this method, and all longs can be produced; passing 0 here
320     * will return 0.
321     * @param state any long; subsequent calls should change by an odd number, such as with {@code ++state}
322     * @return any long
323     */
324    public static long determine(long state)
325    {
326        return (state = ((state = (state ^ (state << 41 | state >>> 23) ^ (state << 17 | state >>> 47)) * 0x369DEA0F31A53F85L) ^ state >>> 37 ^ state >>> 25) * 0xDB4F0B9175AE2165L) ^ state >>> 28;
327    }
328    
329    /**
330     * Static randomizing method that takes its state as a parameter and limits output to an int between 0 (inclusive)
331     * and bound (exclusive); state is expected to change between calls to this. It is recommended that you use
332     * {@code PulleyRNG.determineBounded(++state, bound)} or {@code PulleyRNG.determineBounded(--state, bound)} to
333     * produce a sequence of different numbers, but you can also use
334     * {@code PulleyRNG.determineBounded(state += 12345L, bound)} or any odd-number increment. All longs are accepted
335     * by this method, but not all ints between 0 and bound are guaranteed to be produced with equal likelihood (for any
336     * odd-number values for bound, this isn't possible for most generators). The bound can be negative.
337     * @param state any long; subsequent calls should change by an odd number, such as with {@code ++state}
338     * @param bound the outer exclusive bound, as an int
339     * @return an int between 0 (inclusive) and bound (exclusive)
340     */
341    public static int determineBounded(long state, final int bound)
342    {
343        return (int)((bound * (((state = ((state = (state ^ (state << 41 | state >>> 23) ^ (state << 17 | state >>> 47)) * 0x369DEA0F31A53F85L) ^ state >>> 37 ^ state >>> 25) * 0xDB4F0B9175AE2165L) ^ state >>> 28) & 0xFFFFFFFFL)) >> 32);
344    }
345
346    /**
347     * Returns a random float that is deterministic based on state; if state is the same on two calls to this, this will
348     * return the same float. This is expected to be called with a changing variable, e.g. {@code determine(++state)},
349     * where the increment for state should be odd but otherwise doesn't really matter. This should tolerate just about
350     * any increment as long as it is odd. The period is 2 to the 64 if you increment or decrement by 1, but there are
351     * only 2 to the 30 possible floats between 0 and 1.
352     * @param state a variable that should be different every time you want a different random result;
353     *              using {@code determine(++state)} is recommended to go forwards or {@code determine(--state)} to
354     *              generate numbers in reverse order
355     * @return a pseudo-random float between 0f (inclusive) and 1f (exclusive), determined by {@code state}
356     */
357    public static float determineFloat(long state) {
358        return (((state = (state ^ (state << 41 | state >>> 23) ^ (state << 17 | state >>> 47)) * 0x369DEA0F31A53F85L) ^ state >>> 37 ^ state >>> 25) * 0xDB4F0B9175AE2165L >>> 40) * 0x1p-24f;
359    }
360
361    /**
362     * Returns a random double that is deterministic based on state; if state is the same on two calls to this, this
363     * will return the same float. This is expected to be called with a changing variable, e.g.
364     * {@code determine(++state)}, where the increment for state should be odd but otherwise doesn't really matter. This
365     * should tolerate just about any increment, as long as it is odd. The period is 2 to the 64 if you increment or
366     * decrement by 1, but there are only 2 to the 62 possible doubles between 0 and 1.
367     * @param state a variable that should be different every time you want a different random result;
368     *              using {@code determine(++state)} is recommended to go forwards or {@code determine(--state)} to
369     *              generate numbers in reverse order
370     * @return a pseudo-random double between 0.0 (inclusive) and 1.0 (exclusive), determined by {@code state}
371     */
372    public static double determineDouble(long state) {
373        return (((state = ((state = (state ^ (state << 41 | state >>> 23) ^ (state << 17 | state >>> 47)) * 0x369DEA0F31A53F85L) ^ state >>> 37 ^ state >>> 25) * 0xDB4F0B9175AE2165L) ^ state >>> 28) & 0x1FFFFFFFFFFFFFL) * 0x1p-53;
374    }
375}