MathgeomGLS

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

37
unit Velthuis.RandomNumbers;
38

39
// TODO: streamline this.
40
// TODO: better random number algorithms
41

42
interface
43

44
type
45
  IRandom = interface
46
    function NextInteger: Integer; overload;
47
    function NextInteger(MaxValue: Integer): Integer; overload;
48
    function NextInteger(MinValue, MaxValue: Integer): Integer; overload;
49
    function NextDouble: Double;
50
    function NextInt64: Int64;
51
    procedure NextBytes(var Bytes: array of Byte);
52
    procedure SetSeed(Seed: Int64);
53
    function GetSeed: Int64;
54
    property Seed: Int64 read GetSeed write SetSeed;
55
  end;
56

57
  /// <summary>Base for 32 bit random number generators implementing IRandom</summary>
58
  TRandomBase = class(TInterfacedObject, IRandom)
59
  protected
60
    // Abstract. Must be overridden.
61
    function Next(Bits: Integer): UInt32; virtual; abstract;
62
    procedure SetSeed(ASeed: Int64); virtual; abstract;
63
    function GetSeed: Int64; virtual; abstract;
64
  public
65
    // Generates exception.
66
    constructor Create;
67

68
    // default implementations.
69
    function NextInteger: Integer; overload;
70
    function NextInteger(MaxValue: Integer): Integer; overload;
71
    function NextInteger(MinValue, MaxValue: Integer): Integer; overload;
72
    procedure NextBytes(var Bytes: array of Byte);
73
    function NextDouble: Double;
74
    function NextInt64: Int64;
75
  end;
76

77
  /// <summary>Base for 64 bit random number generators implementing IRandom</summary>
78
  TRandomBase64 = class(TRandomBase, IRandom)
79
  protected
80
    function Next64(Bits: Integer): UInt64; virtual; abstract;
81
    function Next(Bits: Integer): UInt32; override;
82
  public
83
    procedure NextBytes(var Bytes: array of Byte);
84
    function NextDouble: Double;
85
    function NextInt64: Int64;
86
  end;
87

88
  TRandom = class(TRandomBase, IRandom)
89
  private
90
    FSeed: Int64;       // Only 48 bits are used.
91
  protected
92
    function Next(Bits: Integer): UInt32; override;
93
    procedure SetSeed(ASeed: Int64); override;
94
    function GetSeed: Int64; override;
95
  public
96
    constructor Create(Seed: Int64 = 0);
97
  end;
98

99
  TDelphiRandom = class(TRandomBase, IRandom)
100
  protected
101
    function Next(Bits: Integer): UInt32; override;
102
    procedure SetSeed(ASeed: Int64); override;
103
    function GetSeed: Int64; override;
104
  public
105
    constructor Create; overload;
106
    constructor Create(Seed: Int64); overload;
107
  end;
108

109
implementation
110

111
uses
112
  System.SysUtils, Velthuis.Numerics;
113

114
{ TRandom }
115

116
const
117
  CMultiplier = Int64(6364136223846793005);
118
  CIncrement  = Int64(1442695040888963407);
119
  CSeedSize   = 64 div 8;
120

121
constructor TRandom.Create(Seed: Int64);
122
begin
123
  FSeed := Seed;
124
end;
125

126
function TRandom.Next(Bits: Integer): UInt32;
127
begin
128
{$IFOPT Q+}
129
{$DEFINE HasRangeChecks}
130
{$ENDIF}
131
  FSeed := (FSeed * CMultiplier + CIncrement);
132
  Result := UInt32(FSeed shr (64 - Bits)); // Use the highest bits; Lower bits have lower period.
133
{$IFDEF HasRangeChecks}
134
{$RANGECHECKS ON}
135
{$ENDIF}
136
end;
137

138
function TRandom.GetSeed: Int64;
139
begin
140
  Result := FSeed;
141
end;
142

143
procedure TRandom.SetSeed(ASeed: Int64);
144
begin
145
  FSeed := ASeed;
146
end;
147

148
{ TDelphiRandom }
149

150
constructor TDelphiRandom.Create(Seed: Int64);
151
begin
152
  System.RandSeed := Integer(Seed);
153
end;
154

155
constructor TDelphiRandom.Create;
156
begin
157
  Randomize;
158
end;
159

160
function TDelphiRandom.GetSeed: Int64;
161
begin
162
  Result := System.RandSeed;
163
end;
164

165
function TDelphiRandom.Next(Bits: Integer): UInt32;
166
begin
167
  Result := UInt32(System.RandSeed) shr (32 - Bits);
168
  System.Random;
169
end;
170

171
procedure TDelphiRandom.SetSeed(ASeed: Int64);
172
begin
173
  System.RandSeed := Integer(ASeed);
174
end;
175

176
{ TRandomBase }
177

178
constructor TRandomBase.Create;
179
begin
180
  raise EArgumentException.Create('Seed needs initialization');
181
end;
182

183
procedure TRandomBase.NextBytes(var Bytes: array of Byte);
184
var
185
  Head, Tail: Integer;
186
  N, Rnd, I: Integer;
187
begin
188
  Head := Length(Bytes) div SizeOf(Int32);
189
  Tail := Length(Bytes) mod SizeOf(Int32);
190
  N := 0;
191
  for I := 1 to Head do
192
  begin
193
    Rnd := Next(32);
194
    Bytes[N] := Byte(Rnd);
195
    Bytes[N + 1] := Byte(Rnd shr 8);
196
    Bytes[N + 2] := Byte(Rnd shr 16);
197
    Bytes[N + 3] := Byte(Rnd shr 24);
198
    Inc(N, 4);
199
  end;
200
  Rnd := Next(32);
201
  for I := 1 to Tail do
202
  begin
203
    Bytes[N] := Byte(Rnd);
204
    Rnd := Rnd shr 8;
205
    Inc(N);
206
  end;
207
end;
208

209
function TRandomBase.NextDouble: Double;
210
const
211
  Divisor = UInt64(1) shl 53;
212
begin
213
  Result := (UInt64(Next(26) shl 27) + Next(27)) / Divisor;
214
end;
215

216
function TRandomBase.NextInteger: Integer;
217
begin
218
  Result := Next(32);
219
end;
220

221
function TRandomBase.NextInt64: Int64;
222
begin
223
  Result := Int64(Next(32)) shl 32 + Next(32);
224
end;
225

226
function TRandomBase.NextInteger(MinValue, MaxValue: Integer): Integer;
227
begin
228
  if MinValue < 0 then
229
    raise EArgumentException.Create('MinValue must be positive or 0');
230
  Result := MinValue + NextInteger(MaxValue - MinValue);
231
end;
232

233
function TRandomBase.NextInteger(MaxValue: Integer): Integer;
234
var
235
  Bits: Integer;
236
begin
237
  if MaxValue = 0 then
238
    raise EArgumentException.Create('MaxValue not be 0');
239

240
  if IsPowerOfTwo(MaxValue) then
241
  begin
242
    Bits := Next(31);
243
    Exit((Int64(MaxValue) * Bits) shr 31);
244
  end;
245

246
  repeat
247
    Bits := Next(31);
248
    Result := Bits mod MaxValue;
249
  until (Bits - Result + (MaxValue - 1) >= 0);
250
end;
251

252
{ TRandomBase64 }
253

254
function TRandomBase64.Next(Bits: Integer): UInt32;
255
begin
256
  Result := Next64(Bits + 32) shr 32;
257
end;
258

259
procedure TRandomBase64.NextBytes(var Bytes: array of Byte);
260
var
261
  Head, Tail: Integer;
262
  N, I: Integer;
263
  Rnd: UInt64;
264
begin
265
  Head := Length(Bytes) div SizeOf(UInt64);
266
  Tail := Length(Bytes) mod SizeOf(UInt64);
267
  N := 0;
268
  for I := 1 to Head do
269
  begin
270
    Rnd := Next64(64);
271
    Bytes[N] := Byte(Rnd);
272
    Bytes[N + 1] := Byte(Rnd shr 8);
273
    Bytes[N + 2] := Byte(Rnd shr 16);
274
    Bytes[N + 3] := Byte(Rnd shr 24);
275
    Bytes[N + 4] := Byte(Rnd shr 32);
276
    Bytes[N + 5] := Byte(Rnd shr 40);
277
    Bytes[N + 6] := Byte(Rnd shr 48);
278
    Bytes[N + 7] := Byte(Rnd shr 56);
279
    Inc(N, 8);
280
  end;
281
  Rnd := Next64(64);
282
  for I := 1 to Tail do
283
  begin
284
    Bytes[N] := Byte(Rnd);
285
    Rnd := Rnd shr 8;
286
    Inc(N);
287
  end;
288
end;
289

290
function TRandomBase64.NextDouble: Double;
291
const
292
  Divisor = UInt64(1) shl 53;
293
begin
294
  Result := Next64(53) / Divisor;
295
end;
296

297
function TRandomBase64.NextInt64: Int64;
298
begin
299
  Result := Int64(Next64(64));
300
end;
301

302
end.
303

304

305

306

307