sfmt source code

1 module sfmt;
2 
3 import std.stdio;
4 static import sfmt.internal;
5 alias recursion = sfmt.internal.recursion!(18, 1, 11, 1, masks);
6 import sfmt.internal : func1, func2, idxof, ucent_;
7 
8 import std.algorithm : max, min;
9 
10 version (MT19937)
11 {
12     enum SFMT_MEXP = 19937;
13     enum SFMT_N = (SFMT_MEXP >> 7) + 1;
14     enum SFMT_N64 = SFMT_N << 1;
15     enum SFMT_N32 = SFMT_N << 2;
16     enum SFMT_POS1 = 122;
17     enum masks = [0xdfffffefU, 0xddfecb7fU, 0xbffaffffU, 0xbffffff6U];
18     enum parity = [0x00000001U, 0x00000000U, 0x00000000U, 0x13c9e684U];
19     enum id = "SFMT-19937:122-18-1-11-1:dfffffef-ddfecb7f-bffaffff-bffffff6";
20 
21     struct SFMT
22     {
23         enum id = .id;
24         this (uint seed)
25         {
26             this.seed(seed);
27         }
28         this (uint[] seed)
29         {
30             this.seed(seed);
31         }
32         void printState()
33         {
34             import std.stdio;
35             "state:".writeln;
36             foreach (row; state[0..2]~state[$-2..$])
37                 "%(%08x %)".writefln(row.u32);
38         }
39         void fillState(ubyte b)
40         {
41             ucent_ x;
42             x.u32[0] = b;
43             x.u32[0] = x.u32[0] << 8 | x.u32[0];
44             x.u32[0] = x.u32[0] << 16 | x.u32[0];
45             x.u32[1..$] = x.u32[0];
46             state[] = x;
47         }
48         // checked
49         void seed(uint seed)
50         {
51             uint* psfmt32 = &(state[0].u32[0]);
52             psfmt32[idxof!0] = seed;
53             foreach (i; 1..SFMT_N32)
54                 psfmt32[i.idxof] = 1812433253U * (psfmt32[(i - 1).idxof] ^ (psfmt32[(i - 1).idxof] >> 30)) + i;
55             idx = SFMT_N32;
56             assureLongPeriod;
57         }
58         void seed(uint[] seed)
59         {
60             enum size = SFMT_N * 4;
61             static if (size >= 623)
62                 enum lag = 11;
63             else static if (size >= 68)
64                 enum lag = 7;
65             else static if (size >= 39)
66                 enum lag = 5;
67             else
68                 enum lag = 3;
69             enum mid = (size - lag) / 2;
70             fillState(0x8b);
71             immutable count = seed.length.max(SFMT_N32 - 1);
72             uint* psfmt32 = &(state[0].u32[0]);
73             uint r = func1(psfmt32[idxof!0] ^ psfmt32[idxof!mid] ^ psfmt32[idxof!(SFMT_N32 - 1)]);
74             psfmt32[idxof!mid] += r;
75             r += seed.length;
76             psfmt32[idxof!(mid+lag)] += r;
77             psfmt32[idxof!0] = r;
78 
79             size_t i = 1;
80             foreach (j; 0..count.min(seed.length))
81             {
82                 r = func1(
83                         psfmt32[i.idxof]
84                       ^ psfmt32[((i+mid)%SFMT_N32).idxof]
85                       ^ psfmt32[((i+SFMT_N32-1)%SFMT_N32).idxof]);
86                 psfmt32[((i+mid)%SFMT_N32).idxof] += r;
87                 r += seed[j] + i;
88                 psfmt32[((i+mid+lag)%SFMT_N32).idxof] += r;
89                 psfmt32[i.idxof] = r;
90                 i = (i + 1) % SFMT_N32;
91             }
92             foreach (j; count.min(seed.length)..count)
93             {
94                 r = func1(
95                         psfmt32[i.idxof]
96                       ^ psfmt32[((i+mid)%SFMT_N32).idxof]
97                       ^ psfmt32[((i+SFMT_N32-1)%SFMT_N32).idxof]);
98                 psfmt32[((i+mid)%SFMT_N32).idxof] += r;
99                 r += i;
100                 psfmt32[((i+mid+lag)%SFMT_N32).idxof] += r;
101                 psfmt32[i.idxof] = r;
102                 i = (i + 1) % SFMT_N32;
103             }
104             foreach (j; 0..SFMT_N32)
105             {
106                 r = func2(
107                         psfmt32[i.idxof]
108                       + psfmt32[((i+mid)%SFMT_N32).idxof]
109                       + psfmt32[((i+SFMT_N32-1)%SFMT_N32).idxof]);
110                 psfmt32[((i+mid)%SFMT_N32).idxof] ^= r;
111                 r -= i;
112                 psfmt32[((i+mid+lag)%SFMT_N32).idxof] ^= r;
113                 psfmt32[i.idxof] = r;
114                 i = (i + 1) % SFMT_N32;
115             }
116             idx = SFMT_N32;
117             assureLongPeriod;
118         }
119         version (Big32){} else
120         T next(T)()
121             if (is (T == ulong))
122         {
123             ulong* psfmt64 = &(state[0].u64[0]);
124             assert (idx % 2 == 0, "out of alignment");
125             if (SFMT_N32 <= idx)
126                 generateAll;
127             immutable r = psfmt64[idx / 2];
128             idx += 2;
129             return r;
130         }
131         version (Big64){} else
132         T next(T)()
133             if (is (T == uint))
134         {
135             uint* psfmt32 = &(state[0].u32[0]);
136             if (SFMT_N32 <= idx)
137                 generateAll;
138             immutable r = psfmt32[idx];
139             idx += 1;
140             return r;
141         }
142         T next(T)(size_t size)
143             if (is (T == ulong[]) || is (T == uint[]))
144         {
145             return cast(T)fill(cast(ucent_[])(new T(size)));
146         }
147         private auto fill(ucent_[] array)
148         {
149             immutable size = array.length;
150             recursion(
151                 array[0], state[0],
152                 state[0 + SFMT_POS1],
153                 state[SFMT_N - 2], state[SFMT_N - 1]);
154             recursion(
155                 array[1], state[1],
156                 state[1 + SFMT_POS1],
157                 state[SFMT_N - 1], array[0]);
158 
159             foreach (i; 2 .. SFMT_N-SFMT_POS1)
160             {
161                 recursion(
162                     array[i], state[i],
163                     state[i + SFMT_POS1],
164                     array[i - 2], array[i - 1]);
165             }
166             foreach (i; SFMT_N-SFMT_POS1 .. SFMT_N)
167             {
168                 recursion(
169                     array[i], state[i],
170                     array[i + SFMT_POS1 - SFMT_N],
171                     array[i - 2], array[i - 1]);
172             }
173             foreach (i; SFMT_N .. size-SFMT_N)
174             {
175                 recursion(
176                     array[i], array[i - SFMT_N],
177                     array[i + SFMT_POS1 - SFMT_N],
178                     array[i - 2], array[i - 1]);
179             }
180             foreach (j; 0..ptrdiff_t(2*SFMT_N-size).max(0))
181             {
182                 state[j] = array[j + size - SFMT_N];
183             }
184             size_t j = ptrdiff_t(2*SFMT_N-size).max(0);
185             foreach (i; size-SFMT_N..size)
186             {
187                 recursion(
188                     array[i], array[i - SFMT_N],
189                     array[i + SFMT_POS1 - SFMT_N],
190                     array[i - 2], array[i - 1]);
191                 state[j] = array[i];
192                 j += 1;
193             }
194             return array;
195         }
196         private void generateAll()
197         {
198             recursion(
199                 state[0], state[0],
200                 state[0+SFMT_POS1],
201                 state[SFMT_N - 2], state[SFMT_N - 1]);
202             recursion(
203                 state[1], state[1],
204                 state[1+SFMT_POS1],
205                 state[SFMT_N - 1], state[0]);
206             foreach (i; 2..SFMT_N-SFMT_POS1)
207             {
208                 recursion(
209                     state[i], state[i],
210                     state[i+SFMT_POS1],
211                     state[i - 2], state[i - 1]);
212             }
213             foreach (i; SFMT_N-SFMT_POS1..SFMT_N)
214             {
215                 recursion(
216                     state[i], state[i],
217                     state[i+SFMT_POS1-SFMT_N],
218                     state[i - 2], state[i - 1]);
219             }
220             idx = 0;
221         }
222         ucent_[SFMT_N] state;
223         int idx;
224         /// returns true if modification is done
225         bool assureLongPeriod()
226         {
227             uint inner;
228             uint* psfmt32 = &(state[0].u32[0]);
229             foreach (i; 0..4)
230                 inner ^= psfmt32[i.idxof] ^ parity[i];
231             foreach (i; [16, 8, 4, 2, 1])
232                 inner ^= inner >> i;
233             inner ^= 1;
234             if (inner == 1)
235                 return false;
236             foreach (i; 0..4)
237             {
238                 uint working = 1;
239                 foreach (j; 0..32)
240                 {
241                     if (working & parity[i])
242                     {
243                         psfmt32[i.idxof] ^= working;
244                         return true;
245                     }
246                     working <<= 1;
247                 }
248             }
249             assert (false, "unreachable?");
250         }
251     }
252 }
253 else
254 {
255     static assert (false, "Not supported");
256 }
257 
258 version (BigEndian)
259 {
260     pragma (msg, "not tested");
261     version (Only64bit)
262         version = Big64;
263     else version (With32bit)
264         version = Big32;
265     else static assert (false, "Specify Only64bit or With32bit in BigEndian environment");
266 }
267 version (LittleEndian)
268 {
269     pragma (msg, "supported");
270 }
271 
272 version (Only64bit)
273     version (With32bit)
274         static assert (false, "Specify (at most) one of Only64bit or With32bit");
275