001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.rng.core.source32;
019
020import org.apache.commons.rng.core.util.NumberFactory;
021
022import java.util.Arrays;
023
024/**
025 * Middle Square Weyl Sequence Random Number Generator.
026 *
027 * <p>A fast all-purpose 32-bit generator. Memory footprint is 192 bits and the period is at least
028 * {@code 2^64}.</p>
029 *
030 * <p>Implementation is based on the paper
031 * <a href="https://arxiv.org/abs/1704.00358v3">Middle Square Weyl Sequence RNG</a>.</p>
032 *
033 * @see <a href="https://en.wikipedia.org/wiki/Middle-square_method">Middle Square Method</a>
034 * @since 1.3
035 */
036public class MiddleSquareWeylSequence extends IntProvider {
037    /** Size of the seed array. */
038    private static final int SEED_SIZE = 3;
039    /**
040     * The default seed.
041     * This has a high quality Weyl increment (containing many bit state transitions).
042     */
043    private static final long[] DEFAULT_SEED =
044        {0x012de1babb3c4104L, 0xc8161b4202294965L, 0xb5ad4eceda1ce2a9L};
045
046    /** State of the generator. */
047    private long x;
048    /** State of the Weyl sequence. */
049    private long w;
050    /**
051     * Increment for the Weyl sequence. This must be odd to ensure a full period.
052     *
053     * <p>This is not final to support the restore functionality.</p>
054     */
055    private long s;
056
057    /**
058     * Creates a new instance.
059     *
060     * <p>Note: The generator output quality is highly dependent on the initial seed.
061     * If the generator is seeded poorly (for example with all zeros) it is possible the
062     * generator may output zero for many cycles before the internal state recovers randomness.
063     * The seed elements are used to set:</p>
064     *
065     * <ol>
066     *   <li>The state of the generator
067     *   <li>The state of the Weyl sequence
068     *   <li>The increment of the Weyl sequence
069     * </ol>
070     *
071     * <p>The third element is set to odd to ensure a period of at least 2<sup>64</sup>. If the
072     * increment is of low complexity then the Weyl sequence does not contribute high quality
073     * randomness. It is recommended to use a permutation of 8 hex characters for the upper
074     * and lower 32-bits of the increment.</p>
075     *
076     * <p>The state of the generator is squared during each cycle. There is a possibility that
077     * different seeds can produce the same output, for example 0 and 2<sup>32</sup> produce
078     * the same square. This can be avoided by using the high complexity Weyl increment for the
079     * state seed element.</p>
080     *
081     * @param seed Initial seed.
082     * If the length is larger than 3, only the first 3 elements will
083     * be used; if smaller, the remaining elements will be automatically set.
084     */
085    public MiddleSquareWeylSequence(long[] seed) {
086        if (seed.length < SEED_SIZE) {
087            // Complete the seed with a default to avoid
088            // low complexity Weyl increments.
089            final long[] tmp = Arrays.copyOf(seed, SEED_SIZE);
090            System.arraycopy(DEFAULT_SEED, seed.length, tmp, seed.length, SEED_SIZE - seed.length);
091            setSeedInternal(tmp);
092        } else {
093            setSeedInternal(seed);
094        }
095    }
096
097    /**
098     * Seeds the RNG.
099     *
100     * @param seed Seed.
101     */
102    private void setSeedInternal(long[] seed) {
103        x = seed[0];
104        w = seed[1];
105        // Ensure the increment is odd to provide a maximal period Weyl sequence.
106        this.s = seed[2] | 1L;
107    }
108
109    /** {@inheritDoc} */
110    @Override
111    protected byte[] getStateInternal() {
112        return composeStateInternal(NumberFactory.makeByteArray(new long[] {x, w, s}),
113                                    super.getStateInternal());
114    }
115
116    /** {@inheritDoc} */
117    @Override
118    protected void setStateInternal(byte[] state) {
119        final byte[][] c = splitStateInternal(state, SEED_SIZE * 8);
120        setSeedInternal(NumberFactory.makeLongArray(c[0]));
121        super.setStateInternal(c[1]);
122    }
123
124    /** {@inheritDoc} */
125    @Override
126    public int next() {
127        x *= x;
128        x += w += s;
129        x = (x >>> 32) | (x << 32);
130        return (int) x;
131    }
132
133    /** {@inheritDoc} */
134    @Override
135    public long nextLong() {
136        // Avoid round trip from long to int to long by performing two iterations inline
137        x *= x;
138        x += w += s;
139        final long i1 = x & 0xffffffff00000000L;
140        x = (x >>> 32) | (x << 32);
141        x *= x;
142        x += w += s;
143        final long i2 = x >>> 32;
144        x = i2 | x << 32;
145        return i1 | i2;
146    }
147}