MathgeomGLS

1
unit Velthuis.Numerics;
2

3
interface
4

5
(*
6
  Compare
7
  CompareTo
8
  HashCode
9
  HighestOneBit
10
  LowestOneBit
11
  Reverse
12
  ReverseBytes
13
  RotateLeft
14
  RotateRight
15
  Sign
16
  ToBinaryString
17
  ToHexString
18
  ToOctalString
19
  ToString
20
  ToString(base)
21
*)
22

23
// For Delphi XE3 and up:
24
{$IF CompilerVersion >= 24.0 }
25
  {$LEGACYIFEND ON}
26
{$IFEND}
27

28
// For Delphi XE and up:
29
{$IF CompilerVersion >= 22.0}
30
  {$CODEALIGN 16}
31
  {$ALIGN 16}
32
{$IFEND}
33

34
uses
35
  Math, Types;
36

37
function BitCount(U: UInt8): Integer; overload;
38
function BitCount(U: UInt16): Integer; overload;
39
function BitCount(S: Int32): Integer; overload;
40
function BitCount(U: UInt32): Integer; overload;
41

42
function BitLength(S: Int32): Integer; overload;
43
function BitLength(U: UInt32): Integer; overload;
44

45
function HighestOneBit(S: Int32): Int32; overload;
46
function HighestOneBit(U: UInt32): UInt32; overload;
47

48
function IsPowerOfTwo(S: Int32): Boolean; overload;
49
function IsPowerOfTwo(U: UInt32): Boolean; overload;
50

51
function LowestOneBit(S: Int32): Int32; overload;
52
function LowestOneBit(U: UInt32): UInt32; overload;
53

54
function NumberOfLeadingZeros(U: UInt16): Integer; overload;
55
function NumberOfLeadingZeros(S: Int32): Integer; overload;
56
function NumberOfLeadingZeros(U: UInt32): Integer; overload;
57
function NumberOfTrailingZeros(U: UInt32): Integer; overload;
58

59
function Reverse(U: UInt8): UInt8; overload;
60
function Reverse(U: UInt16): UInt16; overload;
61
function Reverse(S: Int32): Int32; overload;
62
function Reverse(U: UInt32): UInt32; overload;
63

64
function ReverseBytes(S: Int32): Int32; overload;
65
function ReverseBytes(U: UInt32): UInt32; overload;
66

67
function RotateLeft(S: Int32; Distance: Integer): Int32; overload;
68
function RotateLeft(U: UInt32; Distance: Integer): UInt32; overload;
69

70
function RotateRight(S: Int32; Distance: Integer): Int32; overload;
71
function RotateRight(U: UInt32; Distance: Integer): UInt32; overload;
72

73
function Sign(S: Int32): TValueSign;
74

75
function ToBinaryString(S: Int32): string; overload;
76
function ToBinaryString(U: UInt32): string; overload;
77

78
function ToHexString(S: Int32): string; overload;
79
function ToHexString(U: UInt32): string; overload;
80

81
function ToOctalString(S: Int32): string; overload;
82
function ToOctalString(U: UInt32): string; overload;
83

84
function ToString(S: Int32; Base: Byte): string; overload;
85
function ToString(U: UInt32; Base: Byte): string; overload;
86

87
implementation
88

89
// https://en.wikipedia.org/wiki/Find_first_set
90

91
uses
92
  SysUtils;
93

94
const
95
  // Currently not used.
96
  NLZDeBruijn32Mult = $07C4ACDD;
97
  NLZDeBruijn32: array[0..31] of Byte =
98
  (
99
    31, 22, 30, 21, 18, 10, 29,  2, 20, 17, 15, 13,  9,  6, 28,  1,
100
    23, 19, 11,  3, 16, 14,  7, 24, 12,  4,  8, 25,  5, 26, 27,  0
101
  );
102

103
  NTZDeBruijn32Mult = $077CB531;
104
  NTZDeBruijn32: array[0..31] of Byte =
105
  (
106
     0,  1, 28,  2, 29, 14, 24,  3, 30, 22, 20, 15, 25, 17,  4,  8,
107
    31, 27, 13, 23, 21, 19, 16,  7, 26, 12, 18,  6, 11,  5, 10,  9
108
  );
109

110
  BitCounts: array[0..15] of Byte = (0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4);
111

112
function BitCount(U: UInt8): Integer;
113
begin
114
  Result := BitCounts[U and $0F] + BitCounts[U shr 4];
115
end;
116

117
function BitCount(U: UInt16): Integer;
118
{$IF DEFINED(WIN32)}
119
asm
120
        MOV     DX,AX
121
        SHR     DX,1
122
        AND     DX,$5555
123
        SUB     AX,DX
124
        MOV     DX,AX
125
        AND     AX,$3333
126
        SHR     DX,2
127
        AND     DX,$3333
128
        ADD     AX,DX
129
        MOV     DX,AX
130
        SHR     DX,4
131
        ADD     AX,DX
132
        AND     AX,$0F0F
133
        MOV     DX,AX
134
        SHR     AX,8
135
        ADD     AX,DX
136
        AND     EAX,$7F
137
end;
138
{$ELSEIF DEFINED(WIN64)}
139
asm
140
        .NOFRAME
141

142
        MOV     AX,CX
143
        SHR     CX,1
144
        AND     CX,$5555
145
        SUB     AX,CX
146
        MOV     CX,AX
147
        AND     AX,$3333
148
        SHR     CX,2
149
        AND     CX,$3333
150
        ADD     AX,CX
151
        MOV     CX,AX
152
        SHR     CX,4
153
        ADD     AX,CX
154
        AND     AX,$0F0F
155
        MOV     CX,AX
156
        SHR     AX,8
157
        ADD     AX,CX
158
        AND     EAX,$7F
159
end;
160
{$ELSE PUREPASCAL}
161
begin
162
  U := U - ((U shr 1) and $5555);
163
  U := (U and $3333) + ((U shr 2) and $3333);
164
  U := (U + (U shr 4)) and $0F0F;
165
  U := U + (U shr 8);
166
  Result := U and $7F;
167
end;
168
{$IFEND PUREPASCAL}
169

170
function BitCount(S: Int32): Integer;
171
begin
172
  Result := BitCount(UInt32(S));
173
end;
174

175
// Faster than 16 bit table lookups
176
function BitCount(U: UInt32): Integer;
177
{$IF DEFINED(WIN32)}
178
asm
179
        MOV     EDX,EAX
180
        SHR     EDX,1
181
        AND     EDX,$55555555
182
        SUB     EAX,EDX
183
        MOV     EDX,EAX
184
        AND     EAX,$33333333
185
        SHR     EDX,2
186
        AND     EDX,$33333333
187
        ADD     EAX,EDX
188
        MOV     EDX,EAX
189
        SHR     EDX,4
190
        ADD     EAX,EDX
191
        AND     EAX,$0F0F0F0F
192
        MOV     EDX,EAX
193
        SHR     EAX,8
194
        ADD     EAX,EDX
195
        MOV     EDX,EAX
196
        SHR     EDX,16
197
        ADD     EAX,EDX
198
        AND     EAX,$7F
199
end;
200
{$ELSEIF DEFINED(WIN64)}
201
asm
202
        .NOFRAME
203

204
        MOV     EAX,ECX
205
        SHR     ECX,1
206
        AND     ECX,$55555555
207
        SUB     EAX,ECX
208
        MOV     ECX,EAX
209
        AND     EAX,$33333333
210
        SHR     ECX,2
211
        AND     ECX,$33333333
212
        ADD     EAX,ECX
213
        MOV     ECX,EAX
214
        SHR     ECX,4
215
        ADD     EAX,ECX
216
        AND     EAX,$0F0F0F0F
217
        MOV     ECX,EAX
218
        SHR     EAX,8
219
        ADD     EAX,ECX
220
        MOV     ECX,EAX
221
        SHR     ECX,16
222
        ADD     EAX,ECX
223
        AND     EAX,$7F
224
end;
225
{$ELSE PUREPASCAL}
226
begin
227
  U := U - ((U shr 1) and $55555555);
228
  U := (U and $33333333) + ((U shr 2) and $33333333);
229
  U := (U + (U shr 4)) and $0F0F0F0F;
230
  U := U + (U shr 8);
231
  U := U + (U shr 16);
232
  Result := U and $7F;
233
end;
234
{$IFEND PUREPASCAL}
235

236
function BitLength(S: Int32): Integer;
237
begin
238
  Result := BitLength(UInt32(S));
239
end;
240

241
function BitLength(U: UInt32): Integer;
242
begin
243
  Result := 32 - NumberOfLeadingZeros(U);
244
end;
245

246
function IsPowerOfTwo(S: Int32): Boolean;
247
begin
248
  Result := IsPowerofTwo(UInt32(Abs(S)));
249
end;
250

251
function IsPowerOfTwo(U: UInt32): Boolean;
252
begin
253
  Result := (U and (U - 1)) = 0;
254
end;
255

256
function HighestOneBit(S: Int32): Int32;
257
begin
258
  Result := Int32(HighestOneBit(UInt32(S)));
259
end;
260

261
function HighestOneBit(U: UInt32): UInt32;
262
begin
263
  if U = 0 then
264
    Result := 0
265
  else
266
    Result := UInt32(1) shl (31 - NumberOfLeadingZeros(U));
267
end;
268

269
function LowestOneBit(S: Int32): Int32;
270
begin
271
  Result := Int32(LowestOneBit(UInt32(S)));
272
end;
273

274
function LowestOneBit(U: UInt32): UInt32;
275
begin
276
  Result := U and -Int32(U);
277
end;
278

279
function NumberOfLeadingZeros(U: UInt16): Integer;
280
{$IF DEFINED(WIN32)}
281
asm
282
        MOVZX   EAX,AX
283
        BSR     EDX,EAX
284
        JNZ     @Invert
285
        MOV     EAX,16
286
        RET
287

288
@Invert:
289

290
        MOV     EAX,15
291
        SUB     EAX,EDX
292
end;
293
{$ELSEIF DEFINED(WIN64)}
294
asm
295
        .NOFRAME
296

297
        MOVZX   EAX,CX
298
        BSR     ECX,EAX
299
        JNZ     @Invert
300
        MOV     EAX,16
301
        RET
302

303
@Invert:
304

305
        MOV     EAX,15
306
        SUB     EAX,ECX
307
end;
308
{$ELSE PUREPASCAL}
309
begin
310
  if U = 0 then
311
    Result := 16
312
  else
313
  begin
314
    Result := 0;
315
    if U <= High(Word) shr 8 then
316
    begin
317
      Result := Result + 8;
318
      U := U shl 8;
319
    end;
320
    if U <= High(Word) shr 4 then
321
    begin
322
      Result := Result + 4;
323
      U := U shl 4;
324
    end;
325
    if U <= High(Word) shr 2 then
326
    begin
327
      Result := Result + 2;
328
      U := U shl 2;
329
    end;
330
    if U <= High(Word) shr 1 then
331
      Result := Result + 1;
332
  end;
333
end;
334
{$IFEND PUREPASCAL}
335

336
function NumberOfLeadingZeros(S: Int32): Integer;
337
begin
338
  Result := NumberOfLeadingZeros(UInt32(S));
339
end;
340

341
function NumberOfLeadingZeros(U: UInt32): Integer;
342
{$IF DEFINED(WIN32)}
343
asm
344
        BSR     EDX,EAX
345
        JNZ     @Invert
346
        MOV     EAX,32
347
        RET
348

349
@Invert:
350

351
        MOV     EAX,31
352
        SUB     EAX,EDX
353

354
@Exit:
355
end;
356
{$ELSEIF DEFINED(WIN64)}
357
asm
358
         .NOFRAME
359

360
         BSR    EDX,ECX
361
         JNZ    @Invert
362
         MOV    EAX,32
363
         RET
364

365
@Invert:
366

367
         MOV    EAX,31
368
         SUB    EAX,EDX
369

370
@Exit:
371
end;
372
{$ELSE PUREPASCAL}
373

374
// Faster than X := X or X shr 1..16; Result := NLZDeBruijn32[...];
375

376
begin
377
  if U = 0 then
378
    Result := 32
379
  else
380
  begin
381
    Result := 0;
382
    if U <= High(Cardinal) shr 16 then
383
    begin
384
      Result := Result + 16;
385
      U := U shl 16;
386
    end;
387
    if U <= High(Cardinal) shr 8 then
388
    begin
389
      Result := Result + 8;
390
      U := U shl 8;
391
    end;
392
    if U <= High(Cardinal) shr 4 then
393
    begin
394
      Result := Result + 4;
395
      U := U shl 4;
396
    end;
397
    if U <= High(Cardinal) shr 2 then
398
    begin
399
      Result := Result + 2;
400
      U := U shl 2;
401
    end;
402
    if U <= High(Cardinal) shr 1 then
403
      Result := Result + 1;
404
  end;
405
end;
406
{$IFEND PUREPASCAL}
407

408
// Faster than NumberOfTrailingZeros2().
409
function NumberOfTrailingZeros(U: UInt32): Integer;
410
{$IF DEFINED(WIN32)}
411
asm
412
        BSF     EAX,EAX
413
        JNZ     @Exit
414
        MOV     EAX,32
415

416
@Exit:
417
end;
418
{$ELSEIF DEFINED(WIN64)}
419
asm
420
        .NOFRAME
421

422
        BSF     EAX,ECX
423
        JNZ     @Exit
424
        MOV     EAX,32
425

426
@Exit:
427
end;
428
{$ELSE PUREPASCAL}
429
begin
430
  if U = 0 then
431
    Result := 32
432
  else
433
    Result := NTZDeBruijn32[((U and (-Integer(U))) * NTZDeBruijn32Mult) shr 27];
434
end;
435
{$IFEND PUREPASCAL}
436

437
function Reverse(U: UInt8): UInt8;
438
begin
439
  U := ((U shr 1) and $55) or ((U and $55) shl 1);
440
  U := ((U shr 2) and $33) or ((U and $33) shl 2);
441
  U := (U shr 4) or (U shl 4);
442
  Result := U;
443
end;
444

445
function Reverse(U: UInt16): UInt16;
446
begin
447
  U := ((U shr 1) and $5555) or ((U and $5555) shl 1);
448
  U := ((U shr 2) and $3333) or ((U and $3333) shl 2);
449
  U := ((U shr 4) and $0F0F) or ((U and $0F0F) shl 4);
450
  U := Swap(U);
451
  Result := U;
452
end;
453

454
function Reverse(S: Int32): Int32;
455
begin
456
  Result := Int32(Reverse(UInt32(S)));
457
end;
458

459
// See http://stackoverflow.com/questions/746171/best-algorithm-for-bit-reversal-from-msb-lsb-to-lsb-msb-in-c too.
460
// http://stackoverflow.com/a/9144870/95954
461
function Reverse(U: UInt32): UInt32;
462
begin
463
  U := ((U shr 1) and $55555555) or ((U and $55555555) shl 1);  // Swap adjacent bits.
464
  U := ((U shr 2) and $33333333) or ((U and $33333333) shl 2);  // Swap adjacent bit pairs.
465
  U := ((U shr 4) and $0F0F0F0F) or ((U and $0F0F0F0F) shl 4);  // Swap nibbles.
466
  U := ((U shr 8) and $00FF00FF) or ((U and $00FF00FF) shl 8);  // Swap bytes.
467
  U := (U shr 16) or (U shl 16);                                // Swap words.
468
  Result := U;
469
end;
470

471
function ReverseBytes(S: Int32): Int32;
472
begin
473
  Result := Int32(ReverseBytes(UInt32(S)));
474
end;
475

476
// Byte and word swaps of Reverse(U).
477
function ReverseBytes(U: UInt32): UInt32;
478
begin
479
  U := ((U shr 8) and $00FF00FF) or ((U and $00FF00FF) shl 8);  // Swap bytes.
480
  U := (U shr 16) or (U shl 16);                                // Swap words.
481
  Result := U;
482
end;
483

484
function RotateLeft(S: Int32; Distance: Integer): Int32;
485
begin
486
  Result := Int32(RotateLeft(UInt32(S), Distance));
487
end;
488

489
function RotateLeft(U: UInt32; Distance: Integer): UInt32;
490
begin
491
  Distance := Distance and 31;
492
  Result := (U shl Distance) or (U shr (32 - Distance));
493
end;
494

495
function RotateRight(S: Int32; Distance: Integer): Int32;
496
begin
497
  Result := Int32(RotateRight(UInt32(S), Distance));
498
end;
499

500
function RotateRight(U: UInt32; Distance: Integer): UInt32;
501
begin
502
  Distance := Distance and 31;
503
  Result := (U shr Distance) or (U shl (32- Distance));
504
end;
505

506
function Sign(S: Int32): TValueSign;
507
begin
508
  Result := Math.Sign(S);
509
end;
510

511
function ToBinaryString(S: Int32): string;
512
begin
513
  Result := ToString(S, 2);
514
end;
515

516
function ToBinaryString(U: UInt32): string;
517
begin
518
  Result := ToString(U, 2);
519
end;
520

521
function ToHexString(S: Int32): string;
522
begin
523
  Result := ToString(S, 16);
524
end;
525

526
function ToHexString(U: UInt32): string;
527
begin
528
  Result := ToString(U, 16);
529
end;
530

531
function ToOctalString(S: Int32): string;
532
begin
533
  Result := ToString(S, 8);
534
end;
535

536
function ToOctalString(U: UInt32): string;
537
begin
538
  Result := ToString(U, 8);
539
end;
540

541
const
542
  Digits: array[0..35] of Char = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ';
543

544
function ToString(S: Int32; Base: Byte): string;
545
begin
546
  if S < 0 then
547
    Result := '-' + ToString(UInt32(Abs(S)), Base)
548
  else
549
    Result := ToString(UInt32(S), Base);
550
end;
551

552
function ToString(U: UInt32; Base: Byte): string;
553
begin
554
  if not (Base in [2..36]) then
555
    raise Exception.Create('Error Message');  // convert error? argument error?
556

557
  if U = 0 then
558
    Result := '0'
559
  else
560
  begin
561
    Result := '';
562
    while U > 0 do
563
    begin
564
      Result := Digits[U mod Base] + Result;
565
      U := U div Base;
566
    end;
567
  end;
568
end;
569

570
end.
571