MathgeomGLS

1
{---------------------------------------------------------------------------}
2
{                                                                           }
3
{ File:       Velthuis.XorShifts.pas                                        }
4
{ Function:   Simple xorshift random number generators, implementing        }
5
{             IRandom interface from Velthuis.RandomNumbers.                }
6
{ Language:   Delphi version XE3 or later                                   }
7
{ Author:     Rudy Velthuis                                                 }
8
{ Copyright:  (c) 2018 Rudy Velthuis                                        }
9
{                                                                           }
10
{ Literature: https://de.wikipedia.org/wiki/Xorshift                        }
11
{             https://en.wikipedia.org/wiki/Xorshift                        }
12
{                                                                           }
13
{ Acknowledgement:                                                          }
14
{             Several of the algorithms below were developed by             }
15
{             Sebastiano Vigna and released to the public domain,           }
16
{             see http://vigna.di.unimi.it/                                 }
17
{                                                                           }
18
{ License:    Redistribution and use in source and binary forms, with or    }
19
{             without modification, are permitted provided that the         }
20
{             following conditions are met:                                 }
21
{                                                                           }
22
{             * Redistributions of source code must retain the above        }
23
{               copyright notice, this list of conditions and the following }
24
{               disclaimer.                                                 }
25
{             * Redistributions in binary form must reproduce the above     }
26
{               copyright notice, this list of conditions and the following }
27
{               disclaimer in the documentation and/or other materials      }
28
{               provided with the distribution.                             }
29
{                                                                           }
30
{ Disclaimer: THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER "AS IS"     }
31
{             AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT     }
32
{             LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND     }
33
{             FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO        }
34
{             EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE     }
35
{             FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,     }
36
{             OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,      }
37
{             PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,     }
38
{             DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED    }
39
{             AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT   }
40
{             LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)        }
41
{             ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF   }
42
{             ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                    }
43
{                                                                           }
44
{---------------------------------------------------------------------------}
45

46
unit Velthuis.XorShifts;
47

48
interface
49

50
uses
51
  Velthuis.RandomNumbers;
52

53
type
54
  TXorShift32 = class(TRandomBase)
55
  private
56
    FSeed: UInt32;
57
  protected
58
    function GetSeed: Int64; override;
59
    function Next(Bits: Integer): UInt32; override;
60
    procedure SetSeed(Seed: Int64); override;
61
  public
62
    constructor Create;
63
  end;
64

65
  TXorShift64 = class(TRandomBase64)
66
  private
67
    FSeed: UInt64;
68
  protected
69
    function GetSeed: Int64; override;
70
    function Next64(Bits: Integer): UInt64; override;
71
    procedure SetSeed(Seed: Int64); override;
72
  public
73
    constructor Create;
74
  end;
75

76
  TXorShift128 = class(TRandomBase)
77
  private
78
    FSeed: array[0..3] of UInt32;
79
    FSeedIndex: Integer;
80
  protected
81
    function GetSeed: Int64; override;
82
    function Next(Bits: Integer): UInt32; override;
83
    procedure SetSeed(Seed: Int64); override;
84
  public
85
    constructor Create;
86
  end;
87

88
  TXorWowState = array[0..4] of UInt32;
89

90
  TXorWow = class(TRandomBase)
91
  private
92
    FSeed: TXorWowState;
93
    FSeedIndex: Integer;
94
  protected
95
    function GetSeed: Int64; override;
96
    function Next(Bits: Integer): UInt32; override;
97
    procedure SetSeed(Seed: Int64); override;
98
  public
99
    constructor Create(const State: TXorWowState);
100
  end;
101

102
  TXorShift64Star = class(TRandomBase64)
103
  private
104
    FSeed: UInt64;
105
  public
106
    constructor Create(const State: UInt64);
107
    function GetSeed: Int64; override;
108
    function Next64(Bits: Integer): UInt64; override;
109
    procedure SetSeed(Seed: Int64); override;
110
  end;
111

112
  TXorShift1024Star = class(TRandomBase64)
113
  private
114
    FSeed: array[0..15] of UInt64;
115
    FSeedIndex: Integer;
116
    FNextIndex: Integer;
117
  protected
118
    function GetSeed: Int64; override;
119
    function Next64(Bits: Integer): UInt64; override;
120
    procedure SetSeed(Seed: Int64); override;
121
  public
122
    constructor Create(State: array of UInt64);
123
  end;
124

125
  TXorShift128Plus = class(TRandomBase64)
126
  private
127
    FSeed: array[0..1] of UInt64;
128
    FSeedIndex: Integer;
129
  protected
130
    function GetSeed: Int64; override;
131
    function Next64(Bits: Integer): UInt64; override;
132
    procedure SetSeed(Seed: Int64); override;
133
  public
134
    constructor Create(State0, State1: UInt64);
135
  end;
136

137
implementation
138

139
uses
140
  System.Math, Winapi.Windows;
141

142
{$RANGECHECKS OFF}
143
{$OVERFLOWCHECKS OFF}
144

145
function SplitMix64(var X: UInt64) : UInt64;
146
var
147
  Z: UInt64;
148
begin
149
  Inc(X, UInt64($9E3779B97F4A7C15));
150
  Z := (X xor (X shr 30)) * UInt64($BF58476D1CE4E5B9);
151
  Z := (Z xor (Z shr 27)) * UInt64($94D049BB133111EB);
152
  Result := Z xor (Z shr 31);
153
end;
154

155
{ TXorShift32 }
156

157
constructor TXorShift32.Create;
158
var
159
  C: Int64;
160
begin
161
  if QueryPerformanceCounter(C) then
162
    FSeed := UInt32(C)
163
  else
164
    FSeed := GetTickCount;
165
end;
166

167
function TXorShift32.GetSeed: Int64;
168
begin
169
  Result := FSeed;
170
end;
171

172
function TXorShift32.Next(Bits: Integer): UInt32;
173
begin
174
  FSeed := FSeed xor (FSeed shl 13);
175
  FSeed := FSeed xor (FSeed shr 17);
176
  FSeed := FSeed xor (FSeed shl 5);
177
  Result := FSeed shr (32 - Bits);
178
end;
179

180
procedure TXorShift32.SetSeed(Seed: Int64);
181
begin
182
  FSeed := UInt32(Seed);
183
end;
184

185
{ TXorShift64 }
186

187
constructor TXorShift64.Create;
188
var
189
  C: Int64;
190
begin
191
  if QueryPerformanceCounter(C) then
192
    FSeed := C
193
  else
194
    FSeed := 88172645463325252 + GetTickCount;
195
end;
196

197
function TXorShift64.GetSeed: Int64;
198
begin
199
  Result := Int64(FSeed);
200
end;
201

202
function TXorShift64.Next64(Bits: Integer): UInt64;
203
begin
204
  FSeed := FSeed xor (FSeed shl 13);
205
  FSeed := FSeed xor (FSeed shr 7);
206
  FSeed := FSeed xor (FSeed shl 17);
207
  Result := FSeed shr (64 - Bits);
208
end;
209

210
procedure TXorShift64.SetSeed(Seed: Int64);
211
begin
212
  FSeed := UInt64(Seed);
213
end;
214

215
{ TXorShift128 }
216

217
constructor TXorShift128.Create;
218
begin
219
  FSeed[0] := 123456789;
220
  FSeed[1] := 362436069;
221
  FSeed[2] := 521288629;
222
  FSeed[3] := 88675123;
223
  FSeedIndex := 0;
224
end;
225

226
function TXorShift128.GetSeed: Int64;
227
begin
228
  Result := Int64(FSeed[1]) shl 32 + FSeed[0];
229
end;
230

231
function TXorShift128.Next(Bits: Integer): UInt32;
232
const
233
  X = 0;
234
  y = 1;
235
  z = 2;
236
  w = 3;
237
var
238
  T: UInt32;
239
begin
240
  T := FSeed[x] xor (FSeed[x] shl 11);
241
  FSeed[x] := FSeed[y];
242
  FSeed[y] := FSeed[z];
243
  FSeed[z] := FSeed[w];
244
  FSeed[w] := FSeed[w] xor ((FSeed[w] shr 19) xor T xor (T shr 8));
245

246
  Result := FSeed[w] shr (32 - Bits);
247
end;
248

249
// Call twice to set full seed.
250
procedure TXorShift128.SetSeed(Seed: Int64);
251
begin
252
  FSeed[FSeedIndex] := UInt32(Seed);
253
  FSeed[FSeedIndex + 1] := UInt32(Seed shr 32);
254
  FSeedIndex := (FSeedIndex + 2) and 3;
255
end;
256

257
{ TXorWow }
258

259
constructor TXorWow.Create(const State: TXorWowState);
260
begin
261
  FSeed := State;
262
  FSeedIndex := 0;
263
end;
264

265
function TXorWow.GetSeed: Int64;
266
begin
267
  Result := Int64(FSeed[1]) shr 32 + FSeed[0];
268
end;
269

270
function TXorWow.Next(Bits: Integer): UInt32;
271
var
272
  S, T: UInt32;
273
begin
274
  T := FSeed[3];
275
  T := T xor (T shr 2);
276
  T := T xor (T shl 1);
277
  FSeed[3] := FSeed[2];
278
  FSeed[2] := FSeed[1];
279
  FSeed[1] := FSeed[0];
280
  S := FSeed[0];
281
  T := T xor S;
282
  T := T xor (S shl 4);
283
  FSeed[0] := T;
284
  FSeed[4] := FSeed[4] + 362437;
285
  Result := (T + FSeed[4]) shr (32 - Bits);
286
end;
287

288
// Call thrice to set full seed.
289
procedure TXorWow.SetSeed(Seed: Int64);
290
begin
291
  if FSeedIndex = 4 then
292
  begin
293
    FSeed[4] := UInt32(Seed);
294
    FSeedIndex := 0;
295
  end
296
  else
297
  begin
298
    FSeed[FSeedIndex] := UInt32(Seed);
299
    FSeed[FSeedIndex + 1] := UInt32(Seed shr 32);
300
    FSeedIndex := FSeedIndex + 2;
301
  end;
302
end;
303

304
{ TXorShift64Star }
305

306
constructor TXorShift64Star.Create(const State: UInt64);
307
begin
308
  FSeed := State;
309
end;
310

311
function TXorShift64Star.GetSeed: Int64;
312
begin
313
  Result := Int64(FSeed);
314
end;
315

316
function TXorShift64Star.Next64(Bits: Integer): UInt64;
317
var
318
  X: UInt64;
319
begin
320
  X := FSeed;
321
  X := X xor (X shr 12);
322
  X := X xor (X shl 25);
323
  X := X xor (X shr 27);
324
  FSeed := X;
325
  Result := (X * UInt64($2545F4914F6CDD1D)) shr (64 - Bits);
326
end;
327

328
procedure TXorShift64Star.SetSeed(Seed: Int64);
329
begin
330
  FSeed := UInt64(Seed);
331
end;
332

333
{ TXorShift1024Star }
334

335
constructor TXorShift1024Star.Create(State: array of UInt64);
336
var
337
  I: Integer;
338
begin
339
  for I := 0 to Max(High(State), High(FSeed)) do
340
    FSeed[I] := State[I];
341
  FSeedIndex := 0;
342
  FNextIndex := 0;
343
end;
344

345
function TXorShift1024Star.GetSeed: Int64;
346
begin
347
  Result := Int64(FSeed[FSeedIndex]);
348
end;
349

350
function TXorShift1024Star.Next64(Bits: Integer): UInt64;
351
var
352
  S0, S1: UInt64;
353
begin
354
  S0 := FSeed[FNextIndex];
355
  FNextIndex := (FNextIndex + 1) and 15;
356
  S1 := FSeed[FNextIndex];
357
  S1 := S1 xor (S1 shl 31);
358
  S1 := S1 xor (S1 shr 11);
359
  S1 := S1 xor (S0 xor (S0 shr 30));
360
  FSeed[FNextIndex] := S1;
361

362
  Result := (S1 * UInt64(1181783497276652981)) shr (64 - Bits);
363
end;
364

365
procedure TXorShift1024Star.SetSeed(Seed: Int64);
366
begin
367
  FSeed[FSeedIndex] := UInt64(Seed);
368
  FSeedIndex := (FSeedIndex + 1) and 15;
369
end;
370

371
{ TXorShift182Plus }
372

373
constructor TXorShift128Plus.Create(State0, State1: UInt64);
374
begin
375
  FSeed[0] := State0;
376
  FSeed[1] := State1;
377
  FSeedIndex := 0;
378
end;
379

380
function TXorShift128Plus.GetSeed: Int64;
381
begin
382
  Result := FSeed[0];
383
end;
384

385
function TXorShift128Plus.Next64(Bits: Integer): UInt64;
386
var
387
  X, Y: UInt64;
388
begin
389
  X := FSeed[0];
390
  Y := FSeed[1];
391
  FSeed[0] := Y;
392
  X := X xor (X shl 23);
393
  FSeed[1] := X xor Y xor (X shr 17) xor (Y shr 26);
394

395
  Result := (FSeed[1] + Y) shr (64 - Bits);
396
end;
397

398
procedure TXorShift128Plus.SetSeed(Seed: Int64);
399
begin
400
  FSeed[FSeedIndex] := Seed;
401
  FSeedIndex := FSeedIndex xor 1;
402
end;
403

404
end.
405