001/*
002 * Ported to Java from the PCG library. Its copyright header follows:
003 *
004 * PCG Random Number Generation for C++
005 *
006 * Copyright 2014 Melissa O'Neill <oneill@pcg-random.org>
007 *
008 * Licensed under the Apache License, Version 2.0 (the "License");
009 * you may not use this file except in compliance with the License.
010 * You may obtain a copy of the License at
011 *
012 *     http://www.apache.org/licenses/LICENSE-2.0
013 *
014 * Unless required by applicable law or agreed to in writing, software
015 * distributed under the License is distributed on an "AS IS" BASIS,
016 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
017 * See the License for the specific language governing permissions and
018 * limitations under the License.
019 *
020 * For additional information about the PCG random number generation scheme,
021 * including its license and other licensing options, visit
022 *
023 *     http://www.pcg-random.org
024 */
025package squidpony.squidmath;
026
027import squidpony.StringKit;
028
029/**
030 * This is a RandomnessSource in the PCG-Random family. It performs pseudo-
031 * random modifications to the output based on the techniques from the
032 * Permuted Congruential Generators created by M.E. O'Neill.
033 * Specifically, this variant is:
034 * RXS M XS -- random xorshift, mcg multiply, fixed xorshift
035 *
036 * The most statistically powerful generator, but all those steps
037 * make it slower than some of the others.
038 *
039 * Because it's usually used in contexts where the state type and the
040 * result type are the same, it is a permutation and is thus invert-able.
041 * We thus provide a (protected) function to invert it.
042 * <br>
043 * Even though benchmarks on similar programs in C would lead you to
044 * believe this should be somewhat faster than LightRNG, benchmarking
045 * with JMH seems to show LightRNG being between 2 and 3 times faster
046 * than PermutedRNG, and both drastically faster than java.util.Random .
047 * @author Melissa E. O'Neill (Go HMC!)
048 * @author Tommy Ettinger
049 */
050public class PermutedRNG implements RandomnessSource, StatefulRandomness
051{
052        /** 2 raised to the 53, - 1. */
053    private static final long DOUBLE_MASK = ( 1L << 53 ) - 1;
054    /** 2 raised to the -53. */
055    private static final double NORM_53 = 1. / ( 1L << 53 );
056    /** 2 raised to the 24, -1. */
057    private static final long FLOAT_MASK = ( 1L << 24 ) - 1;
058    /** 2 raised to the -24. */
059    private static final double NORM_24 = 1. / ( 1L << 24 );
060    /**
061     * The state can be seeded with any value.
062     */
063    public long state;
064
065        private static final long serialVersionUID = 3637443966125527620L;
066
067    /** Creates a new generator seeded using Math.random. */
068    public PermutedRNG() {
069        this((long)Math.floor(Math.random() * Long.MAX_VALUE));
070    }
071
072    public PermutedRNG(final long seed) {
073        state = (seed + 1442695040888963407L) * 6364136223846793005L + 1442695040888963407L;
074    }
075
076    @Override
077    public int next( int bits ) {
078        return (int)( nextLong() & ( 1L << bits ) - 1 );
079    }
080
081
082    /**
083     * From the PCG-Random source:
084     * XorShifts are invert-able, but they are something of a pain to invert.
085     * This function backs them out.
086     * @param n a XorShift-ed value
087     * @param bits the number of bits we still need to invert, not constant
088     * @param shift the crazy one; the wild-card; it's some weird value every time it's used
089     * @return a long that inverts the shift done to n
090     */
091    private static long unxorshift(long n, int bits, int shift)
092    {
093        if (2*shift >= bits) {
094            return n ^ (n >>> shift);
095        }
096        long lowmask1 = (1L << (bits - shift*2)) - 1;
097        long highmask1 = ~lowmask1;
098        long top1 = n;
099        long bottom1 = n & lowmask1;
100        top1 ^= top1 >>> shift;
101        top1 &= highmask1;
102        n = top1 | bottom1;
103        long lowmask2 = (1L << (bits - shift)) - 1;
104        long bottom2 = n & lowmask2;
105        bottom2 = unxorshift(bottom2, bits - shift, shift);
106        bottom2 &= lowmask1;
107        return top1 | bottom2;
108    }
109
110    /**
111     * Can return any int, positive or negative, of any size permissible in a 32-bit signed integer.
112     * Calls nextLong() exactly one time.
113     * @return any int, all 32 bits are random
114     */
115    public int nextInt() {
116        return (int)nextLong();
117    }
118    /**
119     * Can return any long, positive or negative, of any size permissible in a 64-bit signed integer.
120     *
121     * @return any long, all 64 bits are random
122     */
123    @Override
124    public long nextLong()
125    {
126        // increment  = 1442695040888963407L;
127        // multiplier = 6364136223846793005L;
128
129        long p = state;
130        p ^= p >>> (5 + ((p >>> 59) & 31));
131        p *= -5840758589994634535L;
132        state = state * 6364136223846793005L + 1442695040888963407L;
133        return p ^ (p >>> 43);
134    }
135
136    /**
137     * Produces a copy of this RandomnessSource that, if next() and/or nextLong() are called on this object and the
138     * copy, both will generate the same sequence of random numbers from the point copy() was called. This just need to
139     * copy the state so it isn't shared, usually, and produce a new value with the same exact state.
140     *
141     * @return a copy of this RandomnessSource
142     */
143    @Override
144    public RandomnessSource copy() {
145        PermutedRNG next = new PermutedRNG(state);
146        next.setState(state);
147        return next;
148    }
149
150    private static long permute(long p)
151    {
152        p ^= p >>> (5 + ((p >>> 59) & 31));
153        p *= -5840758589994634535L;
154        return p ^ (p >>> 43);
155
156    }
157
158    protected static long invert(long internal)
159    {
160        internal = unxorshift(internal, 64, 43);
161        internal *= -3437190434928431767L;
162        return unxorshift(internal, 64, 5 + ((int)(internal >>> 59) & 31));
163    }
164
165    /**
166     * Exclusive on the upper bound n.  The lower bound is 0.
167     * Will call nextLong() with no arguments at least 1 time, possibly more.
168     * @param n the upper bound; should be positive
169     * @return a random int less than n and at least equal to 0
170     */
171    public int nextInt( final int n ) {
172        if ( n <= 0 ) return 0;
173        int threshold = (0x7fffffff - n + 1) % n;
174        for (;;) {
175            int bits = (int)(nextLong() & 0x7fffffff);
176            if (bits >= threshold)
177                return bits % n;
178        }
179    }
180
181    /**
182     * Inclusive lower, exclusive upper.
183     * Will call nextLong() with no arguments at least 1 time, possibly more.
184     * @param lower the lower bound, inclusive, can be positive or negative
185     * @param upper the upper bound, exclusive, should be positive, must be greater than lower
186     * @return a random int at least equal to lower and less than upper
187     */
188    public int nextInt( final int lower, final int upper ) {
189        if ( upper - lower <= 0 ) throw new IllegalArgumentException();
190        return lower + nextInt(upper - lower);
191    }
192
193    /**
194     * Exclusive on the upper bound n. The lower bound is 0.
195     *
196     * Will call nextLong() with no arguments at least 1 time, possibly more.
197     * @param n the upper bound; should be positive
198     * @return a random long less than n
199     */
200    public long nextLong( final long n ) {
201        if ( n <= 0 ) return 0;
202        long threshold = (0x7fffffffffffffffL - n + 1) % n;
203        for (;;) {
204            long bits = nextLong() & 0x7fffffffffffffffL;
205            if (bits >= threshold)
206                return bits % n;
207        }
208    }
209
210    /**
211     * Exclusive on the upper bound n. The lower bound is 0.
212     *
213     * Will call nextLong() at least 1 time, possibly more.
214     * @param lower the lower bound, inclusive, can be positive or negative
215     * @param upper the upper bound, exclusive, should be positive, must be greater than lower
216     * @return a random long at least equal to lower and less than upper
217     */
218    public long nextLong( final long lower, final long upper ) {
219        if ( upper - lower <= 0 ) return 0;
220        return lower + nextLong(upper - lower);
221    }
222
223    /**
224     * Gets a uniform random double in the range [0.0,1.0)
225     *
226     * Calls nextLong() exactly one time.
227     *
228     * @return a random double at least equal to 0.0 and less than 1.0
229     */
230    public double nextDouble() {
231        return ( nextLong() & DOUBLE_MASK ) * NORM_53;
232    }
233
234    /**
235     * Gets a uniform random double in the range [0.0,outer) given a positive parameter outer. If outer
236     * is negative, it will be the (exclusive) lower bound and 0.0 will be the (inclusive) upper bound.
237     *
238     * Calls nextLong() exactly one time.
239     *
240     *  @param outer the exclusive outer bound, can be negative
241     * @return a random double between 0.0 (inclusive) and outer (exclusive)
242     */
243    public double nextDouble(final double outer) {
244        return nextDouble() * outer;
245    }
246
247    /**
248     * Gets a uniform random float in the range [0.0,1.0)
249     *
250     * Calls nextLong() exactly one time.
251     *
252     * @return a random float at least equal to 0.0f and less than 1.0f
253     */
254    public float nextFloat() {
255        return (float)( ( nextLong() & FLOAT_MASK ) * NORM_24 );
256    }
257
258    /**
259     * Gets a random value, true or false.
260     * Calls nextLong() once.
261     * @return a random true or false value.
262     */
263    public boolean nextBoolean() {
264        return ( nextLong() & 1L ) != 0L;
265    }
266
267    /**
268     * Given a byte array as a parameter, this will fill the array with random bytes (modifying it
269     * in-place). Calls nextLong() {@code Math.ceil(bytes.length / 8.0)} times.
270     * @param bytes a byte array that will have its contents overwritten with random bytes.
271     */
272    public void nextBytes( final byte[] bytes ) {
273        int i = bytes.length, n = 0;
274        while( i != 0 ) {
275            n = Math.min(i, 8 );
276            for ( long bits = nextLong(); n-- != 0; bits >>>= 8 ) bytes[ --i ] = (byte)bits;
277        }
278    }
279
280
281    /**
282     * Sets the seed of this generator (which is also the current state).
283     * @param seed the seed to use for this PermutedRNG, as if it was constructed with this seed.
284     */
285    public void setSeed( final long seed ) {
286        state = seed;
287    }
288    /**
289     * Sets the seed (also the current state) of this generator.
290     * @param seed the seed to use for this PermutedRNG, as if it was constructed with this seed.
291     */
292    @Override
293        public void setState( final long seed ) {
294        state = seed;
295    }
296    /**
297     * Gets the current state of this generator.
298     * @return the current seed of this PermutedRNG, changed once per call to nextLong()
299     */
300    @Override
301        public long getState( ) {
302        return state;
303    }
304
305    /**
306     * Advances or rolls back the PermutedRNG's state without actually generating numbers. Skip forward
307     * or backward a number of steps specified by advance, where a step is equal to one call to nextLong().
308     * @param advance Number of future generations to skip past. Can be negative to backtrack.
309     * @return the state after skipping.
310     */
311    public long skip(long advance)
312    {
313        // The method used here is based on Brown, "Random Number Generation
314        // with Arbitrary Stride,", Transactions of the American Nuclear
315        // Society (Nov. 1994).  The algorithm is very similar to fast
316        // exponentiation.
317        //
318        // Even though advance is a signed long, it is treated as unsigned, effectively, for the purposes
319        // of how many iterations it goes through (at most 63 for forwards, 64 for "backwards").
320        if(advance == 0)
321            return state;
322        long acc_mult = 1, acc_plus = 0, cur_mult = 6364136223846793005L, cur_plus = 1442695040888963407L;
323
324        do {
325            if ((advance & 1L) != 0L) {
326                acc_mult *= cur_mult;
327                acc_plus = acc_plus*cur_mult + cur_plus;
328            }
329            cur_plus = (cur_mult+1L)*cur_plus;
330            cur_mult *= cur_mult;
331            advance >>>= 1;
332        }while (advance > 0L);
333        return acc_mult * state + acc_plus;
334    }
335
336    @Override
337    public String toString() {
338        return "PermutedRNG with state 0x" + StringKit.hex(state) + 'L';
339    }
340
341}