1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 package org.apache.commons.rng.core.source64; 18 19 import java.util.Arrays; 20 import org.apache.commons.rng.core.util.NumberFactory; 21 22 /** 23 * This class provides the 64-bits version of the originally 32-bits 24 * {@link org.apache.commons.rng.core.source32.MersenneTwister 25 * Mersenne Twister}. 26 * 27 * <p> 28 * This class is mainly a Java port of 29 * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt64.html"> 30 * the 2014/2/23 version of the generator 31 * </a> written in C by Takuji Nishimura and Makoto Matsumoto. 32 * </p> 33 * 34 * <p> 35 * Here is their original copyright: 36 * </p> 37 * 38 * <table border="0" width="80%" cellpadding="10" style="background-color: #E0E0E0" summary="Mersenne Twister licence"> 39 * <tr><td>Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura, 40 * All rights reserved.</td></tr> 41 * 42 * <tr><td>Redistribution and use in source and binary forms, with or without 43 * modification, are permitted provided that the following conditions 44 * are met: 45 * <ol> 46 * <li>Redistributions of source code must retain the above copyright 47 * notice, this list of conditions and the following disclaimer.</li> 48 * <li>Redistributions in binary form must reproduce the above copyright 49 * notice, this list of conditions and the following disclaimer in the 50 * documentation and/or other materials provided with the distribution.</li> 51 * <li>The names of its contributors may not be used to endorse or promote 52 * products derived from this software without specific prior written 53 * permission.</li> 54 * </ol></td></tr> 55 * 56 * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND 57 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, 58 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 59 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 60 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS 61 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, 62 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 63 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 64 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY 65 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 66 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE 67 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH 68 * DAMAGE.</strong></td></tr> 69 * </table> 70 * 71 * @see <a href="https://en.wikipedia.org/wiki/Mersenne_Twister">Mersenne Twister (Wikipedia)</a> 72 * @since 1.0 73 */ 74 public class MersenneTwister64 extends LongProvider { 75 /** Size of the bytes pool. */ 76 private static final int NN = 312; 77 /** Period second parameter. */ 78 private static final int MM = 156; 79 /** X * MATRIX_A for X = {0, 1}. */ 80 private static final long[] MAG01 = { 0x0, 0xb5026f5aa96619e9L }; 81 /** Most significant 33 bits. */ 82 private static final long UM = 0xffffffff80000000L; 83 /** Least significant 31 bits. */ 84 private static final long LM = 0x7fffffffL; 85 /** Bytes pool. */ 86 private long[] mt = new long[NN]; 87 /** Current index in the bytes pool. */ 88 private int mti; 89 90 /** 91 * Creates a new random number generator. 92 * 93 * @param seed Initial seed. 94 */ 95 public MersenneTwister64(long[] seed) { 96 setSeedInternal(seed); 97 } 98 99 /** {@inheritDoc} */ 100 @Override 101 protected byte[] getStateInternal() { 102 final long[] s = Arrays.copyOf(mt, NN + 1); 103 s[NN] = mti; 104 105 return composeStateInternal(NumberFactory.makeByteArray(s), 106 super.getStateInternal()); 107 } 108 109 /** {@inheritDoc} */ 110 @Override 111 protected void setStateInternal(byte[] s) { 112 final byte[][] c = splitStateInternal(s, (NN + 1) * 8); 113 114 final long[] tmp = NumberFactory.makeLongArray(c[0]); 115 System.arraycopy(tmp, 0, mt, 0, NN); 116 mti = (int) tmp[NN]; 117 118 super.setStateInternal(c[1]); 119 } 120 121 /** 122 * Initializes the generator with the given seed. 123 * 124 * @param seed Initial seed. 125 */ 126 private void setSeedInternal(long[] seed) { 127 if (seed.length == 0) { 128 // Accept empty seed. 129 seed = new long[1]; 130 } 131 132 initState(19650218L); 133 int i = 1; 134 int j = 0; 135 136 for (int k = Math.max(NN, seed.length); k != 0; k--) { 137 final long mm1 = mt[i - 1]; 138 mt[i] = (mt[i] ^ ((mm1 ^ (mm1 >>> 62)) * 0x369dea0f31a53f85L)) + seed[j] + j; // non linear 139 i++; 140 j++; 141 if (i >= NN) { 142 mt[0] = mt[NN - 1]; 143 i = 1; 144 } 145 if (j >= seed.length) { 146 j = 0; 147 } 148 } 149 for (int k = NN - 1; k != 0; k--) { 150 final long mm1 = mt[i - 1]; 151 mt[i] = (mt[i] ^ ((mm1 ^ (mm1 >>> 62)) * 0x27bb2ee687b0b0fdL)) - i; // non linear 152 i++; 153 if (i >= NN) { 154 mt[0] = mt[NN - 1]; 155 i = 1; 156 } 157 } 158 159 mt[0] = 0x8000000000000000L; // MSB is 1; assuring non-zero initial array 160 } 161 162 /** 163 * Initialize the internal state of this instance. 164 * 165 * @param seed Seed. 166 */ 167 private void initState(long seed) { 168 mt[0] = seed; 169 for (mti = 1; mti < NN; mti++) { 170 final long mm1 = mt[mti - 1]; 171 mt[mti] = 0x5851f42d4c957f2dL * (mm1 ^ (mm1 >>> 62)) + mti; 172 } 173 } 174 175 /** {@inheritDoc} */ 176 @Override 177 public long next() { 178 long x; 179 180 if (mti >= NN) { // generate NN words at one time 181 for (int i = 0; i < NN - MM; i++) { 182 x = (mt[i] & UM) | (mt[i + 1] & LM); 183 mt[i] = mt[i + MM] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)]; 184 } 185 for (int i = NN - MM; i < NN - 1; i++) { 186 x = (mt[i] & UM) | (mt[i + 1] & LM); 187 mt[i] = mt[ i + (MM - NN)] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)]; 188 } 189 190 x = (mt[NN - 1] & UM) | (mt[0] & LM); 191 mt[NN - 1] = mt[MM - 1] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)]; 192 193 mti = 0; 194 } 195 196 x = mt[mti++]; 197 198 x ^= (x >>> 29) & 0x5555555555555555L; 199 x ^= (x << 17) & 0x71d67fffeda60000L; 200 x ^= (x << 37) & 0xfff7eee000000000L; 201 x ^= x >>> 43; 202 203 return x; 204 } 205 }