You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2019/08/06 10:42:24 UTC

[commons-rng] 01/07: RNG-85: Added Middle-Square Weyl Sequence generator.

This is an automated email from the ASF dual-hosted git repository.

aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-rng.git

commit 7314a5af2dabc45303fc602f69ab054c53c92e04
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Mon Aug 5 11:00:41 2019 +0100

    RNG-85: Added Middle-Square Weyl Sequence generator.
---
 .../core/source32/MiddleSquareWeylSequence.java    | 122 +++++++++++++++++++++
 .../org/apache/commons/rng/core/ProvidersList.java |   3 +
 .../source32/MiddleSquareWeylSequenceTest.java     |  81 ++++++++++++++
 3 files changed, 206 insertions(+)

diff --git a/commons-rng-core/src/main/java/org/apache/commons/rng/core/source32/MiddleSquareWeylSequence.java b/commons-rng-core/src/main/java/org/apache/commons/rng/core/source32/MiddleSquareWeylSequence.java
new file mode 100644
index 0000000..b84cd47
--- /dev/null
+++ b/commons-rng-core/src/main/java/org/apache/commons/rng/core/source32/MiddleSquareWeylSequence.java
@@ -0,0 +1,122 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.rng.core.source32;
+
+import org.apache.commons.rng.core.util.NumberFactory;
+
+/**
+ * Middle Square Weyl Sequence Random Number Generator.
+ *
+ * <p>A fast all-purpose 32-bit generator. Memory footprint is 192 bits and the period is at least
+ * {@code 2^64}.</p>
+ *
+ * <p>Implementation is based on the paper
+ * <a href="https://arxiv.org/abs/1704.00358v3">Middle Square Weyl Sequence RNG</a>.</p>
+ *
+ * @see <a href="https://en.wikipedia.org/wiki/Middle-square_method">Middle Square Method</a>
+ *
+ * @since 1.3
+ */
+public class MiddleSquareWeylSequence extends IntProvider {
+    /** Size of the seed array. */
+    private static final int SEED_SIZE = 3;
+
+    /** State of the generator. */
+    private long x;
+    /** State of the Weyl sequence. */
+    private long w;
+    /**
+     * Increment for the Weyl sequence. This must be odd to ensure a full period.
+     *
+     * <p>This is not final to support the restore functionality.</p>
+     */
+    private long s;
+
+    /**
+     * Creates a new instance.
+     *
+     * <p>Note: The generator output quality is highly dependent on the initial seed.
+     * If the generator is seeded poorly (for example with all zeros) it is possible the
+     * generator may output zero for many cycles before the internal state recovers randomness.
+     * The seed elements are used to set:</p>
+     *
+     * <ol>
+     *   <li>The state of the generator
+     *   <li>The state of the Weyl sequence
+     *   <li>The increment of the Weyl sequence
+     * </ol>
+     *
+     * <p>The third element is set to odd to ensure a period of at least 2<sup>64</sup>. If the
+     * increment is of low complexity then the Weyl sequence does not contribute high quality
+     * randomness. It is recommended to use a permutation of 8 hex characters for the upper
+     * and lower 32-bits of the increment.</p>
+     *
+     * <p>The state of the generator is squared during each cycle. There is a possibility that
+     * different seeds can produce the same output, for example 0 and 2<sup>32</sup> produce
+     * the same square. This can be avoided by using the high complexity Weyl increment for the
+     * state seed element.</p>
+     *
+     * @param seed Initial seed.
+     * If the length is larger than 3, only the first 3 elements will
+     * be used; if smaller, the remaining elements will be automatically set.
+     */
+    public MiddleSquareWeylSequence(long[] seed) {
+        if (seed.length < SEED_SIZE) {
+            final long[] tmp = new long[SEED_SIZE];
+            fillState(tmp, seed);
+            setSeedInternal(tmp);
+        } else {
+            setSeedInternal(seed);
+        }
+    }
+
+    /**
+     * Seeds the RNG.
+     *
+     * @param seed Seed.
+     */
+    private void setSeedInternal(long[] seed) {
+        x = seed[0];
+        w = seed[1];
+        // Ensure the increment is odd to provide a maximal period Weyl sequence.
+        this.s = seed[2] | 1L;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected byte[] getStateInternal() {
+        return composeStateInternal(NumberFactory.makeByteArray(new long[] {x, w, s}),
+                                    super.getStateInternal());
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected void setStateInternal(byte[] state) {
+        final byte[][] c = splitStateInternal(state, SEED_SIZE * 8);
+        setSeedInternal(NumberFactory.makeLongArray(c[0]));
+        super.setStateInternal(c[1]);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public int next() {
+        x *= x;
+        x += w += s;
+        return (int) (x = (x >>> 32) | (x << 32));
+    }
+}
diff --git a/commons-rng-core/src/test/java/org/apache/commons/rng/core/ProvidersList.java b/commons-rng-core/src/test/java/org/apache/commons/rng/core/ProvidersList.java
index 1c00bfb..9d4b6b5 100644
--- a/commons-rng-core/src/test/java/org/apache/commons/rng/core/ProvidersList.java
+++ b/commons-rng-core/src/test/java/org/apache/commons/rng/core/ProvidersList.java
@@ -34,6 +34,7 @@ import org.apache.commons.rng.core.source32.Well44497a;
 import org.apache.commons.rng.core.source32.Well44497b;
 import org.apache.commons.rng.core.source32.ISAACRandom;
 import org.apache.commons.rng.core.source32.MersenneTwister;
+import org.apache.commons.rng.core.source32.MiddleSquareWeylSequence;
 import org.apache.commons.rng.core.source32.MultiplyWithCarry256;
 import org.apache.commons.rng.core.source32.KISSRandom;
 import org.apache.commons.rng.core.source32.PcgXshRr32;
@@ -103,6 +104,8 @@ public final class ProvidersList {
             add(LIST32, new PcgXshRs32(new long[] {g.nextLong()}));
             add(LIST32, new PcgMcgXshRr32(g.nextLong()));
             add(LIST32, new PcgMcgXshRs32(g.nextLong()));
+            // Ensure a high complexity increment is used for the Weyl sequence
+            add(LIST32, new MiddleSquareWeylSequence(new long[] {g.nextLong(), g.nextLong(), 0xb5ad4eceda1ce2a9L}));
             // ... add more here.
 
             // "long"-based RNGs.
diff --git a/commons-rng-core/src/test/java/org/apache/commons/rng/core/source32/MiddleSquareWeylSequenceTest.java b/commons-rng-core/src/test/java/org/apache/commons/rng/core/source32/MiddleSquareWeylSequenceTest.java
new file mode 100644
index 0000000..7a6c590
--- /dev/null
+++ b/commons-rng-core/src/test/java/org/apache/commons/rng/core/source32/MiddleSquareWeylSequenceTest.java
@@ -0,0 +1,81 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.rng.core.source32;
+
+import org.apache.commons.rng.core.RandomAssert;
+import org.junit.Test;
+
+public class MiddleSquareWeylSequenceTest {
+    @Test
+    public void testReferenceCode() {
+        /*
+         * The data was generated using the following program based on the author's C code:
+         *     https://mswsrng.wixsite.com/rand
+         *
+         * #include <stdio.h>
+         * #include <stdint.h>
+         *
+         * uint64_t x, w, s;
+         *
+         * inline static uint32_t msws() {
+         *     x *= x; x += (w += s); return x = (x>>32) | (x<<32);
+         * }
+         *
+         * int main() {
+         *     x = 0x012de1babb3c4104;
+         *     w = 0xc8161b4202294965;
+         *     s = 0xb5ad4eceda1ce2a9;
+         *     for (int i=0; i<10; i++) {
+         *         for (int j=0; j<4; j++) {
+         *             printf("0x%08x, ", msws());
+         *         }
+         *         printf("\n");
+         *     }
+         * }
+         */
+        final long[] seed = {0x012de1babb3c4104L, 0xc8161b4202294965L, 0xb5ad4eceda1ce2a9L};
+
+        final int[] expectedSequence = {
+            0xe7f4010b, 0x37bdb1e7, 0x05d8934f, 0x22970c75,
+            0xe7432a9f, 0xd157c60f, 0x26e9b5ae, 0x3dd91250,
+            0x8dbf85f1, 0x99e3aa17, 0xcb90322b, 0x29a007e2,
+            0x25a431fb, 0xcc612768, 0x510db5cd, 0xeb0aec2f,
+            0x05f88c18, 0xcdb79066, 0x5222c513, 0x9075045c,
+            0xf11a0e0e, 0x0106ab1d, 0xe2546700, 0xdf0a7656,
+            0x170e7908, 0x17a7b775, 0x98d69720, 0x74da3b78,
+            0x410ea18e, 0x4f708277, 0x471853e8, 0xa2cd2587,
+            0x16238d96, 0x57653154, 0x7ecbf9c8, 0xc5dd75bf,
+            0x32ed82a2, 0x4700e664, 0xb0ad77c9, 0xfb87df7b,
+        };
+
+        RandomAssert.assertEquals(expectedSequence, new MiddleSquareWeylSequence(seed));
+    }
+
+    /**
+     * Test the self-seeding functionality. The seeding of the generator requires a high complexity
+     * increment for the Weyl sequence. The standard test for the statistical quality of the
+     * generator uses a good increment. This test ensures the generator can be created with a seed
+     * smaller than the seed size without exception; the statistical quality of the output is not
+     * assessed and expected to be poor.
+     */
+    @Test
+    public void testSelfSeeding() {
+        // Ensure this does not throw. The output quality will be poor.
+        final MiddleSquareWeylSequence rng = new MiddleSquareWeylSequence(new long[0]);
+        rng.nextInt();
+    }
+}