-
Notifications
You must be signed in to change notification settings - Fork 14
/
bigdecimalmath.pas
2502 lines (2201 loc) · 86.7 KB
/
bigdecimalmath.pas
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
{
Copyright (C) 2013 Benito van der Zander (BeniBela)
benito@benibela.de
www.benibela.de
This file is distributed under under the same license as Lazarus and the LCL itself:
This file is distributed under the Library GNU General Public License
with the following modification:
As a special exception, the copyright holders of this library give you
permission to link this library with independent modules to produce an
executable, regardless of the license terms of these independent modules,
and to copy and distribute the resulting executable under terms of your choice,
provided that you also meet, for each linked independent module, the terms
and conditions of the license of that module. An independent module is a
module which is not derived from or based on this library. If you modify this
library, you may extend this exception to your version of the library, but
you are not obligated to do so. If you do not wish to do so, delete this
exception statement from your version.
}
(***
A unit for arbitrary precision arithmetics on BCD floats
See @link(BigDecimal)
*)
unit bigdecimalmath;
{$mode objfpc}{$H+}
{$ModeSwitch advancedrecords}
{$COperators on}
interface
uses
Classes, SysUtils, math;
const PACKAGE_VERSION = '0.9.0.repo';
{$DEFINE USE_9_DIGITS}
{$IF defined(USE_1_DIGIT) or defined(USE_1_DIGITS)}
const DIGITS_PER_ELEMENT = 1;
const ELEMENT_OVERFLOW = 10;
type BigDecimalBin = shortint; BigDecimalBinSquared = longint; //must be large enough to store ELEMENT_OVERFLOW*ELEMENT_OVERFLOW
{$ELSEIF defined(USE_2_DIGITS)}
const DIGITS_PER_ELEMENT = 2;
const ELEMENT_OVERFLOW = 100;
type BigDecimalBin = smallint {shortint is to small to store overflow during addition}; BigDecimalBinSquared = longint; //must be large enough to store ELEMENT_OVERFLOW*ELEMENT_OVERFLOW
{$ELSEIF defined(USE_3_DIGITS)}
const DIGITS_PER_ELEMENT = 3;
const ELEMENT_OVERFLOW = 1000;
type BigDecimalBin = smallint; BigDecimalBinSquared = longint; //must be large enough to store ELEMENT_OVERFLOW*ELEMENT_OVERFLOW
{$ELSEIF defined(USE_4_DIGITS)}
const DIGITS_PER_ELEMENT = 4;
const ELEMENT_OVERFLOW = 10000;
type BigDecimalBin = smallint; BigDecimalBinSquared = longint; //must be large enough to store ELEMENT_OVERFLOW*ELEMENT_OVERFLOW
{$ELSEIF defined(USE_5_DIGITS)}
const DIGITS_PER_ELEMENT = 5;
const ELEMENT_OVERFLOW = 100000;
type BigDecimalBin = integer; BigDecimalBinSquared = int64; //must be large enough to store ELEMENT_OVERFLOW*ELEMENT_OVERFLOW
{$ELSEIF defined(USE_6_DIGITS)}
const DIGITS_PER_ELEMENT = 6;
const ELEMENT_OVERFLOW = 1000000;
type BigDecimalBin = integer; BigDecimalBinSquared = int64; //must be large enough to store ELEMENT_OVERFLOW*ELEMENT_OVERFLOW
{$ELSEIF defined(USE_7_DIGITS)}
const DIGITS_PER_ELEMENT = 7;
const ELEMENT_OVERFLOW = 10000000;
type BigDecimalBin = integer; BigDecimalBinSquared = int64; //must be large enough to store ELEMENT_OVERFLOW*ELEMENT_OVERFLOW
{$ELSEIF defined(USE_8_DIGITS)}
const DIGITS_PER_ELEMENT = 8;
const ELEMENT_OVERFLOW = 100000000;
type BigDecimalBin = integer; BigDecimalBinSquared = int64; //must be large enough to store ELEMENT_OVERFLOW*ELEMENT_OVERFLOW
{$ELSEIF defined(USE_9_DIGITS)}
const DIGITS_PER_ELEMENT = 9;
const ELEMENT_OVERFLOW = 1000000000;
type BigDecimalBin = integer; BigDecimalBinSquared = int64; //must be large enough to store ELEMENT_OVERFLOW*ELEMENT_OVERFLOW
{$ELSE}
Invalid digit count
{$ENDIF}
type TBigDecimalFormat = (bdfExact, bdfExponent);
type TBigDecimalRoundingMode = (bfrmTrunc, bfrmCeil, bfrmFloor, bfrmRound, bfrmRoundHalfUp, bfrmRoundHalfToEven);
type
//** @abstract(Big Decimal type). @br
//** Consisting of an bcd integer times a decimal exponent ([integer digits] * 10 ^ (DIGITS_PER_ELEMENT * exponent)) @br
//** It can be used like a normal floating point number. E.g: @longCode(#
//** var bd: BigDecimal;
//** bd := 12.34;
//** bd := bd * 1000 - 42; // bd = 12298
//** bd := bd / 7.0; // bd = 1756.85714285714286
//** bd := StrToBigDecimal('123456789012345678901234567890123456789') + 1; // bd = 123456789012345678901234567890123456790
//** #) @br@br
//** It has an arbitrary precision (up to 18 billion digits), and can be converted to a decimal string without loss of precision, since
//** it stores decimal digits (up to 9 digits / array element, depending on compiler define). @br
BigDecimal = record
digits: array of BigDecimalBin;
exponent: integer;
signed, lastDigitHidden: ByteBool;
//** Returns true iff the bigdecimal is zero
function isZero(): boolean;
//** Returns true iff v has no fractional digits
function isIntegral(): boolean;
//** Returns true iff v has no fractional digits and can be stored within an longint (32 bit integer)
function isLongint(): boolean;
//** Returns true iff v has no fractional digits and can be stored within an int64
function isInt64(): boolean;
//** Checks if v is odd. A number with fractional digits is never odd (only weird)
function isOdd(): boolean;
//** Checks if v is even. A number with fractional digits is never even (and neither odd, which is odd)
function isEven(): boolean;
//** How many non-zero digits the number contains
function precision(): integer;
//** The index of the leading, most significant digit @br
//** That is, the exponent of number when it is written in scientific notation @br
//** That is, 10^result <= v < 10^(result+1)
function mostSignificantExponent(): integer;
//** Returns the digit-th digit of v. @br
//** Last integer digit is digit 0, digits at negative indices are behind the decimal point.
function getDigit(digit: integer): BigDecimalBin;
//** Compares the big decimals. Returns -1, 0 or 1 corresponding to a <, = or > b
class function compare(const a, b: BigDecimal): integer; static;
function toLongint: longint;
function toInt64: int64;
function toSizeInt: sizeint;
function tryToLongint(out v: longint): boolean;
function tryToInt64(out v: int64): boolean;
function tryToSizeInt(out v: sizeint): boolean;
function toString(format: TBigDecimalFormat = bdfExact): string;
{$ifdef FPC_HAS_TYPE_SINGLE}
function toSingle: single;
{$endif}
{$ifdef FPC_HAS_TYPE_Double}
function toDouble: double;
{$endif}
{$ifdef FPC_HAS_TYPE_EXTENDED}
function toExtended: extended;
{$endif}
//** Sets the bigdecimal to 0
procedure setZero();
//** Sets the bigdecimal to 1
procedure setOne();
//** Removes leading (pre .) and trailing (post .) zeros
procedure normalize();
//** Universal rounding function @br
//** Rounds to the precision of a certain digit, subject to a certain rounding mode. @br
//** Positive toDigit will round to an integer with toDigit trailing zeros, negative toDigit will round to a decimal with -toDigit numbers after the decimal point
function rounded(toDigit: integer = 0; roundingMode: TBigDecimalRoundingMode = bfrmRound): BigDecimal;
//** Calculates a decimal shift: @code(self := self * 10^shift)
procedure shift10(shift: integer);
//** Calculates a decimal shift: @code(result := self * 10^shift)
function shifted10(shift: integer): BigDecimal;
end;
PBigDecimal = ^BigDecimal;
type TBigDecimalErrorCode = (bdceNoError, bdceParsingInvalidFormat, bdceParsingTooBig );
PBigDecimalErrorCode = ^TBigDecimalErrorCode;
//** Converts a decimal pchar to a bigdecimal. @br
//** Supports standard decimal notation, like -123.456 or 1E-2 (@code(-?[0-9]+(.[0-9]+)?([eE][-+]?[0-9]+)))
function TryStrToBigDecimal(pstart: pchar; length: SizeInt; res: PBigDecimal; errCode: PBigDecimalErrorCode = nil): boolean;
//** Converts a decimal string to a bigdecimal. @br
//** Supports standard decimal notation, like -123.456 or 1E-2 (@code(-?[0-9]+(.[0-9]+)?([eE][-+]?[0-9]+)))
function TryStrToBigDecimal(const s: string; res: PBigDecimal; errCode: PBigDecimalErrorCode = nil): boolean;
//** Converts a decimal string to a bigdecimal. @br
//** Supports standard decimal notation, like -123.456 or 1E-2 (@code(-?[0-9]+(.[0-9]+)?([eE][-+]?[0-9]+)))
//** Raises an exception on invalid input.
function StrToBigDecimal(const s: string): BigDecimal; inline;
//type TBigDecimalFormat = (bdfExact, bdfExponent); format: TBigDecimalFormat = bdfExact
//** Converts a bigdecimal to a decimal string @br
//** The result will be fixed width format [0-9]+(.[0-9]+)?, even if the input had an exponent
function BigDecimalToStr(const v: BigDecimal; format: TBigDecimalFormat = bdfExact): string;
//**Converts a bigdecimal to a native int (can overflow)
function BigDecimalToLongint(const a: BigDecimal): Longint;
//**Converts a bigdecimal to a native int (can overflow)
function BigDecimalToInt64(const a: BigDecimal): Int64;
{$ifdef FPC_HAS_TYPE_Extended}
//**Converts a bigdecimal to an extended (may introduce rounding errors)
function BigDecimalToExtended(const a: BigDecimal): Extended; deprecated 'Use .toExtended record method';
{$endif}
type TBigDecimalFloatFormat = (bdffExact, bdffShortest);
{$ifdef FPC_HAS_TYPE_Double}
function FloatToBigDecimal(const v: Double; format: TBigDecimalFloatFormat = bdffShortest): BigDecimal; overload;
{$endif FPC_HAS_TYPE_Double}
{$ifdef FPC_HAS_TYPE_Single}
function FloatToBigDecimal(const v: Single; format: TBigDecimalFloatFormat = bdffShortest): BigDecimal; overload;
{$endif FPC_HAS_TYPE_Single}
{$ifdef FPC_HAS_TYPE_Extended}
function FloatToBigDecimal(const v: Extended; format: TBigDecimalFloatFormat = bdffShortest): BigDecimal; overload;
{$endif FPC_HAS_TYPE_Extended}
//operator :=(const a: BigDecimal): Integer;
//** Converts a native integer to a BigDecimal
operator :=(const a: Integer): BigDecimal;
//operator :=(const a: BigDecimal): Int64;
//** Converts a native integer to a BigDecimal
operator :=(const a: Int64): BigDecimal;
//operator :=(const a: BigDecimal): QWord;
//** Converts a native integer to a BigDecimal
operator :=(const a: QWord): BigDecimal;
//operator :=(const a: BigDecimal): Extended; auto conversion of bigdecimal to extended is possible, but it confuses fpc overload resolution. Then e.g. power calls either math or bigdecimalbc depending on the unit order in the uses clause
//** Converts an extended to a BigDecimal @br
//** Marked as deprecated, because it may lead to rounding errors. FloatToBigDecimal is exact, but probably some magnitudes slower. For constant values StrToBigDecimal should be used instead.
operator :=(const a: Extended): BigDecimal; deprecated 'Direct casting of float to bigdecimal might lead to rounding errors. Consider using StrToBigDecimal.';
//** Standard operator unary -
operator -(const a: BigDecimal): BigDecimal;
//** Standard operator binary +
operator +(const a: BigDecimal; const b: BigDecimal): BigDecimal;
//** Standard operator binary -
operator -(const a: BigDecimal; const b: BigDecimal): BigDecimal;
//** Standard operator binary *
operator *(const a: BigDecimal; const b: BigDecimal): BigDecimal;
//** Standard operator binary / @br
//** If the result can not be represented as finite decimal number (e.g. 1/3) it will be calculated with 18 digit precision after the decimal
//** point, with an additional hidden digit for rounding (so 1/3 is 0.333333333333333333, and 0.333333333333333333*3 is 0.999999999999999999, but (1/3) * 3 is 1).
operator /(const a: BigDecimal; const b: BigDecimal): BigDecimal;
//** Standard operator binary div @br
//** The result is an integer, so 1.23E3 / 7 will be 175
operator div(const a: BigDecimal; const b: BigDecimal): BigDecimal;
//** Standard operator binary mod @br
//** Calculates the remainder of an integer division @code(a - (a div b) * b)
operator mod(const a: BigDecimal; const b: BigDecimal): BigDecimal;
//** Standard operator binary ** @br
operator **(const a: BigDecimal; const b: int64): BigDecimal;
type TBigDecimalDivisionFlags = set of (bddfKeepDividentPrecision, bddfKeepDivisorPrecision, bddfAddHiddenDigit, bddfFillIntegerPart, bddfNoFractionalPart);
//** Universal division/modulo function. Calculates the quotient and remainder of a / b. @br
//** @param maximalAdditionalFractionDigits How many digits should be added to the quotient, if the result cannot be represented with the current precision
//** @param flags Division options:
//** bddfKeepDividentPrecision: calculates as least as many non-zero digit of the quotient as the divident (1st arg) has @br
//** bddfKeepDivisorPrecision: calculates as least as many non-zero digit of the quotient as the divisor (2nd arg) has @br
//** bddfAddHiddenDigit: Calculates an additional digit for rounding, which will not be displayed by BigDecimalToStr@br
//** bddfFillIntegerPart: Calculate at least all digits of the integer part of the quotient, independent of the precision of the input @br
//** bddfNoFractionalPart: Do not calculate the fractional part of the quotient (remember that a bigdecimal is a scaled integer. So bfdfFillIntegerPart ensures that the result has not less digits than an integer division (necessary in case of an exponent > 0) and bfdfKillFractions that the result has not more digits than an integer division (in case of an exponent < 0) ) @br
//** not all flag combinations were tested
procedure divideModNoAlias(out quotient, remainder: BigDecimal; const a, b: BigDecimal; targetPrecision: integer = 18; flags: TBigDecimalDivisionFlags = [bddfKeepDividentPrecision, bddfKeepDivisorPrecision, bddfAddHiddenDigit]);
//** Wrapper around divideModNoAlias, ignoring the calculated remainder
function divide(const a, b: BigDecimal; maximalAdditionalFractionDigits: integer = 18; flags: TBigDecimalDivisionFlags = [bddfKeepDividentPrecision, bddfKeepDivisorPrecision, bddfAddHiddenDigit]): BigDecimal;
procedure shift10(var v: BigDecimal; shift: integer); deprecated 'use advanced record method';
function shifted10(const v: BigDecimal; shift: integer): BigDecimal; deprecated 'use advanced record method';
function compareBigDecimals(const a, b: BigDecimal): integer; deprecated 'use advanced record method';
operator <(const a: BigDecimal; const b: BigDecimal): boolean;
operator <=(const a: BigDecimal; const b: BigDecimal): boolean;
operator =(const a: BigDecimal; const b: BigDecimal): boolean;
operator >=(const a: BigDecimal; const b: BigDecimal): boolean;
operator >(const a: BigDecimal; const b: BigDecimal): boolean;
procedure normalize(var x: BigDecimal); deprecated 'use advanced record method';
function precision(const v: BigDecimal): integer; deprecated 'use advanced record method';
function mostSignificantExponent(const v: BigDecimal): integer; deprecated 'use advanced record method';
//** Universal rounding function @br
//** Rounds v to the precision of a certain digit, subject to a certain rounding mode. @br
//** Positive toDigit will round to an integer with toDigit trailing zeros, negative toDigit will round to a decimal with -toDigit numbers after the decimal point
function round(const v: BigDecimal; toDigit: integer = 0; roundingMode: TBigDecimalRoundingMode = bfrmRound): BigDecimal; overload;
//**Given mi < exact < ma, truncate exact to a bigdecimal result, such that @br
//** mi < result < ma @br
//** result has the minimal number of non-zero digits @br
//** | result - exact | is minimized
function roundInRange(mi, exact, ma: BigDecimal): BigDecimal;
function getDigit(const v: BigDecimal; digit: integer): BigDecimalBin;
//** Sets the bigdecimal to 0
procedure setZero(out r: BigDecimal); deprecated 'use advanced record method';
//** Sets the bigdecimal to 1
procedure setOne(out r: BigDecimal); deprecated 'use advanced record method';
function isZero(const v: BigDecimal): boolean; overload; deprecated 'use advanced record method';
function isIntegral(const v: BigDecimal): boolean; deprecated 'use advanced record method';
function isLongint(const v: BigDecimal): boolean; deprecated 'use advanced record method';
function isInt64(const v: BigDecimal): boolean; deprecated 'use advanced record method';
function odd(const v: BigDecimal): boolean; overload; deprecated 'use advanced record method';
function even(const v: BigDecimal): boolean; overload; deprecated 'use advanced record method';
//** Returns the absolute value of v
function abs(const v: BigDecimal): BigDecimal; overload;
//** Calculates v ** exp, with exp being an integer
function power(const v: BigDecimal; const exp: Int64): BigDecimal; overload;
//** Calculates the square root of v, to precision digits after the decimal point @br
//** Not much tested
function sqrt(const v: BigDecimal; precision: integer = 9): BigDecimal; overload;
//** Calculates the greatest common denominator (only makes sense for positive integer input)
function gcd(const a,b: BigDecimal): BigDecimal; overload;
//** Calculates the least common multiple
function lcm(const a,b: BigDecimal): BigDecimal; overload;
//** Calculates 2 ** exp exactly, with exp being an integer (faster than power for negative exp)
function fastpower2to(const exp: Int64): BigDecimal;
//** Calculates 5 ** exp exactly, with exp being an integer (faster than power for negative exp)
function fastpower5to(const exp: Int64): BigDecimal;
type TFloatInformation = class
type
{$ifdef FPC_HAS_TYPE_SINGLE}
TSingleInformation = class
type float = single;
const PowersOf10: array[0..10] of single = (1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000,
1e10);
const MaxExactPowerOf10 = 10; //todo: is this correct? (log to base 5)
const MaxExactMantissa = 16777215; //2^24 - 1
const MaxExactMantissaDigits = 8;
end;
{$endif}
{$ifdef FPC_HAS_TYPE_DOUBLE}
TDoubleInformation = class
type float = double;
//numbers with mantissa <= 2^53 - 1 = 9007199254740991 and -22 <= true exponent <= 22 can be represented exactly as double
const PowersOf10: array[0..22] of double = (1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000,
1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22 );
const MaxExactPowerOf10 = 22;
const MaxExactMantissa = 9007199254740991; //2^53 - 1
const MaxExactMantissaDigits = 16;
end;
{$endif}
{$ifdef FPC_HAS_TYPE_EXTENDED}
TExtendedInformation = class
type float = extended;
const PowersOf10: array[0..27] of extended = (1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000,
1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22,
1e23, 1e24,1e25, 1e26,1e27);
const MaxExactPowerOf10 = 27; //todo: is this correct? (log to base 5)
const MaxExactMantissa: QWord = QWord($FFFFFFFFFFFFFFFF); //2^64 - 1 = 18446744073709551615
const MaxExactMantissaDigits = 20;
end;
{$endif}
end;
implementation
const divisionDefaultPrecision = 18;
divisionDefaultFlags = [bddfKeepDividentPrecision,
bddfKeepDivisorPrecision,
bddfAddHiddenDigit,
bddfFillIntegerPart
];
const powersOf10: array[0..9] of longint = (1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000);
//returns:
//@þaram(intstart Position of the first digit in integer part)
//@þaram(intend Position AFTER the last digit in integer part)
//@þaram(dot Position of the . character, or nil)
//@þaram(exp Position of the [eE] character or nil)
function TryStrDecodeDecimal(const pstart, pend: pchar; out intstart, intend, dot, exp: pchar): boolean;
var
p: PChar;
begin
result := false;
if pend <= pstart then exit();
dot := nil;
exp := nil;
p := pstart;
if p^ in ['+', '-'] then inc(p);
intstart := p;
while p < pend do begin
case p^ of
'0'..'9': ;
'.': if (dot <> nil) or (exp <> nil) then exit() else dot := p;
'e', 'E': if exp <> nil then exit() else exp := p;
'+', '-': if p <> exp + 1 then exit();
else exit();
end;
inc(p);
end;
if exp = pstart then exit;
if exp = nil then intend := pend
else intend := exp;
if intend = dot + 1 then begin intend -= 1; dot := nil; end;
if intend <= intstart then exit;
result := true;
end;
function TryStrToBigDecimal(pstart: pchar; length: SizeInt; res: PBigDecimal; errCode: PBigDecimalErrorCode): boolean;
var dot, exp: pchar;
i: pchar;
intstart: pchar;
intend: pchar;
trueexponent: int64;
p: Integer;
j: Integer;
totalintlength, k, code: Integer;
pend: pchar;
expparsing: shortstring;
begin
pend := pstart + length;
result := TryStrDecodeDecimal(pstart, pstart + length, intstart, intend, dot, exp);
if not result then begin
if Assigned(errCode) then errCode^ := bdceParsingInvalidFormat;
exit;
end else if Assigned(errCode) then errCode^ := bdceNoError;
if exp = nil then trueexponent := 0
else begin
inc(exp);
if (pend - exp {length exp until pend} <= 10) and (res = nil) then exit; //if the exponent is small, we know it is okay. If we do not need res, we can exit, otherwise we need to actually get the exponent
expparsing := '';
SetLength(expparsing, pend - exp);
move(exp^, expparsing[1], pend - exp);
val(expparsing, trueexponent, code); //this is faster than using inttostr
if code <> 0 then trueexponent := high(int64);
if (trueexponent < DIGITS_PER_ELEMENT * int64(low(integer))) or (trueexponent > DIGITS_PER_ELEMENT * int64(high(integer))) then begin
dec(exp);
while pstart < exp do begin
if not (pstart^ in ['0', '.', '-']) then begin
//if there is anything non-zero, the exponent is too big
if assigned(errCode) then errCode^ := bdceParsingTooBig;
exit(false);
end;
inc(pstart);
end;
if res <> nil then res^.setZero(); //but if all digits are 0, the exponent can be ignored
exit;
end;
end;
if res = nil then exit;
with res^ do begin
signed := pstart^ = '-';
lastDigitHidden := false;
if dot <> nil then trueexponent -= intend - dot - 1; //shifting the dot to the left corresponds to divisions by 10, so it reduces the exponent
exponent := trueexponent div DIGITS_PER_ELEMENT;
if (trueexponent < 0) and (int64(exponent) * DIGITS_PER_ELEMENT <> trueexponent) then exponent -= 1; //truncate to negative infinity
totalintlength := (intend - intstart) + (trueexponent - int64(exponent) * DIGITS_PER_ELEMENT);
if dot <> nil then totalintlength -= 1;
//parse digits from string
SetLength(digits, (totalintlength + DIGITS_PER_ELEMENT - 1) div DIGITS_PER_ELEMENT);
p := high(digits);
i := intstart;
//if totalintlength is not divisible by DIGITS_PER_ELEMENT, the first and last bin need additional zeros
if totalintlength mod DIGITS_PER_ELEMENT = 0 then j := 1
else j := DIGITS_PER_ELEMENT + 1 - (totalintlength) mod DIGITS_PER_ELEMENT;
while i < intend do begin
digits[p] := 0;
while (i < intend) and (j <= DIGITS_PER_ELEMENT) do begin
if i <> dot then begin
digits[p] := digits[p] * 10 + ord(i^) - ord('0');
j += 1;
end;
i += 1;
end;
k := j;
j := 1;
p -= 1;
end;
digits[0] := digits[0] * powersOf10[DIGITS_PER_ELEMENT - k + 1];
if signed and res^.isZero() then res^.setZero();
end;
end;
function TryStrToBigDecimal(const s: string; res: PBigDecimal; errCode: PBigDecimalErrorCode = nil): boolean;
begin
result := TryStrToBigDecimal(pchar(s), length(s), res, errCode);
end;
function StrToBigDecimal(const s: string): BigDecimal;
begin
if not TryStrToBigDecimal(s, @result) then
raise EConvertError.Create(s +' is not a valid number');
end;
function digitsInBin(i: integer): integer;
begin
if i < 100000 then begin
if i < 100 then begin
if i >= 10 then exit(2)
else exit(1);
end else begin
if i >= 10000 then exit(5)
else if i >= 1000 then exit(4)
else exit(3);
end;
end else begin
if i >= 1000000000 then exit(10)
else if i >= 100000000 then exit(9)
else if i >= 10000000 then exit(8)
else if i >= 1000000 then exit(7)
else exit(6);
end;
end;
function trailingZeros(i: integer): integer;
var j: integer;
begin
if i = 0 then exit(DIGITS_PER_ELEMENT);
j := i div 10;
result := 0;
while (i <> 0) and (i - 10 * j {= i mod 10} = 0) do begin
result += 1;
i := j {= i div 10};
j := i div 10;
end;
end;
procedure skipZeros(const v: BigDecimal; out highskip, lowskip: integer);
var
i: Integer;
begin
with v do begin
highskip := 0;
for i := high(digits) downto max(0, -exponent + 1) do
if digits[i] = 0 then highskip+=1
else break;
lowskip := 0;
for i := 0 to min(high(digits) - highskip, - exponent - 1) do
if digits[i] = 0 then lowskip += 1
else break;
end;
end;
function BigDecimalToStr(const v: BigDecimal; format: TBigDecimalFormat = bdfExact): string;
const BigDecimalDecimalSeparator = '.';
BigDecimalExponent = 'E';
BigDecimalZeroResult = '0'; //return value on zero input. Might depend on format some day (e.g. 0.0, 0.0E0)
procedure intToStrFixedLength(t: integer; var p: pchar; len: integer = DIGITS_PER_ELEMENT); inline;
var
j: Integer;
begin
for j := 1 to len do begin
p^ := chr(t mod 10 + ord('0'));
t := t div 10;
dec(p);
end;
end;
var
lowskip: integer = 0;
lowBinLength: Integer = 0;
lowBin: BigDecimalBin = 0;
displayed: PBigDecimal;
dotBinPos: Integer = 0;
firstHigh: integer = 0;
procedure lowBinTrimTrailingZeros;
begin
while (lowBin mod 10 = 0) and (lowBinLength > 0) do begin
lowBin := lowBin div 10;
lowBinLength -= 1;
end;
end;
procedure setLowBin;
begin
while (lowskip <= firstHigh) and (displayed^.digits[lowskip] = 0) do lowskip += 1;
if lowskip > firstHigh then exit;
lowBinLength := DIGITS_PER_ELEMENT;
lowBin := displayed^.digits[lowskip];
if lowskip <= dotBinPos then lowBinTrimTrailingZeros;
end;
var
skip: Integer = 0; //leading trimmed zeros
procedure init;
begin
with displayed^ do begin
dotBinPos := -exponent - 1; //first bin after the decimal point (every bin i <= dotBinPos belongs to the fractional part)
skipZeros(displayed^, skip, lowskip);
firstHigh:=high(digits) - skip;
if length(digits) = skip + lowskip then exit();
end;
setLowBin;
end;
var
p: PAnsiChar;
i: Integer;
reslen: integer;
highBin: BigDecimalBin;
highBinLength: Integer;
additionalCarry: Boolean;
tempdecimal: BigDecimal;
explength: Integer;
realexponent: int64;
begin
//Algorithm for bdfExact:
//print all numbers bin starting at the lexical last one (bin nr. 0)
// trim trailing 0 after .
// print bins till .
// print 000 for very low/high exp
// print bins before . trimming leading 0
//
//Three cases:
//a) aaa bbb ccc ddd eee fff 000 000 000
// 5 4 3 2 1 0 <-exponent->
//
//b) aaa bbb ccc ddd . eee fff
// 5 4 3 2 1 0
// -exponent (= 2)
//
//c) 000 . 000 000 aaa bbb ccc ddd eee fff
// 5 4 3 2 1 0
// -exponent (8)
displayed := @v;
init;
with v do begin
if length(digits) = skip + lowskip then exit(BigDecimalZeroResult);
//remove last hidden digit, and increment the number by one if the hidden digit is >= 5
if (lastDigitHidden) then begin
additionalCarry := (lowBin mod 10 >= 5) ;
if additionalCarry and (lowBin div 10 + 1 >= powersOf10[lowBinLength-1]) then begin
tempdecimal := v.rounded((exponent + lowskip + 1) * DIGITS_PER_ELEMENT - (lowBinLength - 1));
displayed := @tempdecimal;
init;
if length(displayed^.digits) = skip + lowskip then exit(BigDecimalZeroResult);
end else if (lowskip <= dotBinPos) then begin
lowBin := lowBin div 10;
lowBinLength -= 1;
if additionalCarry then lowBin+=1;
if lowBinLength = 0 then begin
lowskip += 1;
setLowBin;
end;
end else begin
lowBin := lowBin - lowBin mod 10;
if additionalCarry then lowBin+=10;
end;
end;
end;
with displayed^ do begin
case format of
bdfExact: begin
//calculate the length of the result
if firstHigh > dotBinPos then highBin := digits[firstHigh] else highBin := 0;
highBinLength := digitsInBin(highBin);
if dotBinPos < lowskip then reslen := (firstHigh + exponent) * DIGITS_PER_ELEMENT //integer number
else begin
//(each += corresponds to a for loop below)
reslen := lowBinLength ;
reslen += (min(high(digits), dotBinPos) - lowskip) * DIGITS_PER_ELEMENT;
reslen += max(0, dotBinPos - high(digits) ) * DIGITS_PER_ELEMENT;
reslen += max(0, firstHigh - max(-exponent, 0)) * DIGITS_PER_ELEMENT;
if reslen <> 0 then
reslen += 1; //dot
end;
reslen += highBinLength;
if reslen = 0 then exit('0');
if signed then reslen += 1;
//generate result (last digit bin to first digit bin)
SetLength(result, reslen);
p := @result[length(Result)];
if dotBinPos >= lowskip then begin
//fractional part
intToStrFixedLength(lowBin, p, lowBinLength); //last bin (with trimmed trailing zeros)
for i := lowskip + 1 to min(high(digits), dotBinPos) do //other bins
intToStrFixedLength(digits[i], p, DIGITS_PER_ELEMENT);
for i := high(digits)+1 to dotBinPos do begin //additional zeros given by exponent (after .)
p -= DIGITS_PER_ELEMENT;
FillChar((p + 1)^, DIGITS_PER_ELEMENT, '0');
end;
p^ := BigDecimalDecimalSeparator; dec(p);
end;
//additional zeros given by exponent (before .)
for i := 1 to exponent do begin
p -= DIGITS_PER_ELEMENT;
FillChar(p^, DIGITS_PER_ELEMENT + 1, '0');
end;
if lowskip > 0 then lowskip := 0; //print zeros here (they are only skipped for exponent output format)
if (lowskip > dotBinPos) and (lowskip < firstHigh) then begin
lowBin := displayed^.digits[lowskip];
intToStrFixedLength(lowBin, p, DIGITS_PER_ELEMENT);
i := lowskip + 1;
end else i := max(-exponent, 0);
for i := i to firstHigh - 1 do //other bins
intToStrFixedLength(digits[i], p, DIGITS_PER_ELEMENT);
intToStrFixedLength(highBin, p, highBinLength); //first bin (with trimmed leading zeros)
end;
bdfExponent: begin
while (firstHigh >= 0) and (digits[firstHigh] = 0) do dec(firstHigh);
if firstHigh < 0 then exit(BigDecimalZeroResult);
highBin := digits[firstHigh];
highBinLength := digitsInBin(highBin);
lowBinTrimTrailingZeros;
//calculate the length of the result
if lowskip <> firstHigh then begin
reslen := highBinLength + (firstHigh - lowskip - 1) * DIGITS_PER_ELEMENT + lowBinLength;
end else begin
lowBinLength := highBinLength + lowBinLength - DIGITS_PER_ELEMENT;
if lowBinLength = 1 then begin
lowBin := lowBin * 10;
lowBinLength := 2;
end;
reslen := lowBinLength;
end;
reslen += 1; //dot
// if reslen = 2 then reslen += 1; //always something after the dot
realexponent := int64(exponent + firstHigh) * DIGITS_PER_ELEMENT + highBinLength - 1;
if (realexponent >= low(integer)) and (realexponent <= high(integer)) then explength := digitsInBin(abs(realexponent))
else begin
explength := digitsInBin(abs(realexponent) div 1000000000);
reslen += 9;
end;
reslen += 1 + explength; //E...
if realexponent < 0 then reslen+=1; //E-...
if signed then reslen += 1;
//generate result
SetLength(result, reslen);
p := @result[length(Result)];
if (realexponent >= low(integer)) and (realexponent <= high(integer)) then intToStrFixedLength(abs(realexponent), p, explength)
else begin
intToStrFixedLength(abs(realexponent) mod 1000000000, p, 9);
intToStrFixedLength(abs(realexponent) div 1000000000, p, explength);
end;
if realexponent < 0 then begin p^ := '-'; dec(p); end;
p^ := BigDecimalExponent; dec(p);
if lowskip <> firstHigh then begin
intToStrFixedLength(lowBin, p, lowBinLength);
for i := lowskip+1 to firstHigh - 1 do
intToStrFixedLength(digits[i], p);
intToStrFixedLength(highBin, p, highBinLength);
end else
intToStrFixedLength(lowBin, p, lowBinLength);
p^ := (p+1)^;
(p+1)^ := BigDecimalDecimalSeparator;
dec(p);
end;
{$if FPC_FULLVERSION < 030300} else result := ''; p := nil; {$endif}
end;
if signed then begin p^ := '-'; dec(p); end;
//safety check
if p + 1 <> @result[1] then
raise EInvalidOperation.Create('Expected result length wrong');
end;
end;
function BigDecimalToLongint(const a: BigDecimal): Longint;
var
i: Integer;
begin
result := 0;
for i := high(a.digits) downto max(0, - a.exponent) do
result := result * ELEMENT_OVERFLOW - a.digits[i]; //create negative value (as it has a larger range by 1)
if (a.exponent > 0) and (result <> 0) then
for i := 1 to a.exponent do
result := result * ELEMENT_OVERFLOW;
if not a.signed then result := -result;
end;
function BigDecimalToInt64(const a: BigDecimal): Int64;
var
i: Integer;
begin
result := 0;
for i := high(a.digits) downto max(0, - a.exponent) do
result := result * ELEMENT_OVERFLOW - a.digits[i]; //create negative value (as it has a larger range by 1)
if (a.exponent > 0) and (result <> 0) then
for i := 1 to a.exponent do
result := result * ELEMENT_OVERFLOW;
if not a.signed then result := -result;
end;
{$ifdef FPC_HAS_TYPE_Extended}
function BigDecimalToExtended(const a: BigDecimal): Extended;
begin
result := a.toExtended;
end;
{$endif}
function roundInRange(mi, exact, ma: BigDecimal): BigDecimal;
function safebin(const bd: BigDecimal; bini: Integer): BigDecimalBin; inline;
begin
bini -= bd.exponent;
if (bini < 0) or (bini > high(bd.digits)) then result := 0
else result := bd.digits[bini];
end;
procedure findDifferentDigit(const a, b: BigDecimal; highExp, lowExp: integer; out bin, digit: integer; out digitA, digitB: BigDecimalBin);
var aBin, bbin: BigDecimalBin;
atemp, btemp: BigDecimalBin;
bini: Integer;
digiti: Integer;
begin
for bini := highExp downto lowExp do begin
aBin := safebin(a, bini);
bBin := safebin(b, bini);
if aBin = bbin then continue; //exact must be in-between
for digiti := DIGITS_PER_ELEMENT - 1 downto 0 do begin
atemp := aBin div powersOf10[digiti];
btemp := bbin div powersOf10[digiti];
if atemp <> btemp then begin
digitA := atemp mod 10;
digitB := btemp mod 10;
bin := bini;
digit := digiti;
exit;
end;
end;
end;
//the numbers should be different
assert(false);
digit := -1;
end;
function cutoff(bin, digit: integer; digitDelta: BigDecimalBin = 0): BigDecimal;
begin
result := round(exact, bin * DIGITS_PER_ELEMENT + digit, bfrmTrunc);
if digitDelta <> 0 then begin
bin -= result.exponent;
assert(bin >= 0);
if (bin > high(result.digits)) then
SetLength(result.digits, bin + 1);
result.digits[bin] := result.digits[bin] + powersOf10[digit] * digitDelta; //there should be no overflow, since the function is only called if there is a digit to inc/decrement
end;
end;
var
highestExp, highestBin: integer;
digit, bin: Integer;
midigit,exdigit,madigit: BigDecimalBin;
canUseMaDigit: Boolean;
nextdigit: BigDecimalBin;
i: Integer;
mitemp: BigDecimalBin;
extemp: BigDecimalBin;
bin2: Integer;
mibin: BigDecimalBin;
exbin: BigDecimalBin;
begin
if (mi.signed and not mi.isZero()) or ma.signed then begin //if mi is signed 0, treat it as unsigned 0
//if ma is signed, mi must be signed and non zero
if not ma.signed and not ma.isZero() then result.setZero() //0 has the minimal number of non-zero digits, but can only be returned if 0 < ma. Then we know mi < 0 or we would not be here
else begin
//possible cases mi signed and ma not signed => ma = 0
// mi signed and ma signed
assert(exact.signed);
mi.signed := false; exact.signed := false; ma.signed := false;
result := roundInRange(ma, exact, mi);
result.signed := not result.isZero();
end;
exit;
end;
//find the first digit pos di such that mi[1..di] < exact[1..di] < ma[1..di]
//if exact[di+1] >= 5, round up, otherwise truncate
//if that would take it outside the bound, search a later digit
result := exact;
highestExp := min(min(mi.exponent, exact.exponent), ma.exponent);
highestBin := max(max(mi.exponent + high(mi.digits),
exact.exponent + high(exact.digits)),
ma.exponent + high(ma.digits));
findDifferentDigit(mi,ma,highestBin, highestExp, bin, digit, midigit, madigit);
if digit = -1 then exit;
exdigit := (safebin(exact, bin) div powersOf10[digit]) mod 10;
if digit = 0 then nextdigit := safebin(exact, bin-1) div powersOf10[DIGITS_PER_ELEMENT-1]
else nextdigit := (safebin(exact, bin) div powersOf10[digit - 1]) mod 10;
//we know midigit <= exdigit <= madigit and midigit < madigit
if (nextdigit < 5) and (midigit < exdigit) and (exdigit < madigit) then
exit(cutoff(bin, digit)); //if we do not want to round up and the current digit is in range
if exdigit + 1 < madigit then exit(cutoff(bin, digit, + 1)); //if we do want to round up. Or midigit = exdigit so we must round up (we know exdigit < madigit since exdigit + 1 < madigit)
//so exdigit in \{ madigit - 1, madigit \} because exdigit + 1 >= madigit, but exdigit + 1 <= madigit + 1
//check if there is a further non zero in ma. If ma is ...matemp00000000+ then we cannot use madigt
canUseMaDigit := safebin(ma, bin) mod powersOf10[digit] <> 0;
if not canUseMaDigit then
for i := bin - 1 downto ma.exponent do
if safebin(ma, i) <> 0 then begin
canUseMaDigit := true;
break;
end;
if canUseMaDigit then begin
//try exdigit and exdigit + 1, order depending, if we want to round up or not
if nextdigit >= 5 then begin
if exdigit + 1 <= madigit {we know midigit < exdigit +1 } then exit(cutoff(bin, digit, +1));
if midigit < exdigit {we know exdigit <= madigit} then exit(cutoff(bin, digit));
end else begin
if midigit < exdigit {we know exdigit <= madigit} then exit(cutoff(bin, digit));
if exdigit + 1 <= madigit {we know midigit < exdigit +1 } then exit(cutoff(bin, digit, +1));
end;
//now we know exdigit < madigit. Otherwise exdigit = madigit and midigit >= exdigit = madigit, which must not be
//Therefore exdigit + 1 <= madigit.
Assert(false);
exit(exact);
end;
//we cannot round up.. if we could, we would have done so above
if (midigit < exdigit) and (exdigit < madigit) then exit(cutoff(bin, digit));
//Try to round down
if (midigit < exdigit - 1) and (exdigit - 1 < madigit) then exit(cutoff(bin, digit, -1));
//now either (midigit = exdigit) or (exdigit = madigit)
//we know midigit < madigit
//if midigit + 1 < madigit, i.e. madigit >= midigit + 2, then
// in the case exdigit = midigit, exdigit + 1 = midigit + 1 < madigit would have been true above
// in the case exdigit = madigit, exdigit - 1 = madigit - 1 > midigit would have been true in the previous if
assert(midigit + 1 = madigit);
//So the situation is like (123, 1235, 124) and we need additional digits
for bin2 := bin downto highestExp do begin
mibin := safebin(mi, bin2);
exbin := safebin(exact, bin2);
if bin2 = bin then begin
//we already know digit is useless, so fill everrything left from it with 9
mibin := mibin mod powersOf10[digit] + ((ELEMENT_OVERFLOW - 1) - (ELEMENT_OVERFLOW - 1) mod powersOf10[digit]);
exbin := exbin mod powersOf10[digit] + ((ELEMENT_OVERFLOW - 1) - (ELEMENT_OVERFLOW - 1) mod powersOf10[digit]);
end;
if (mibin = exbin) and (mibin = ELEMENT_OVERFLOW - 1) then continue; //exact must be in-between
for digit := DIGITS_PER_ELEMENT - 1 downto 0 do begin
mitemp := mibin div powersOf10[digit];
extemp := exbin div powersOf10[digit];
if mitemp <> extemp then exit(cutoff(bin2, digit)) //different digits, use exact
else if extemp mod 10 <> 9 then exit(cutoff(bin2, digit, +1));
end;
end;
assert(false);
end;
function getDigit(const v: BigDecimal; digit: integer): BigDecimalBin;
begin
result := v.getDigit(digit);
end;
procedure setZero(out r: BigDecimal);
begin
r.setZero();
end;
procedure setOne(out r: BigDecimal);
begin
r.setOne();
end;
{$ifdef FPC_HAS_TYPE_Single}
procedure splitFloat(const v: single; out sign: boolean; out exponent: integer; out mantissa: QWord); overload;
begin
sign := (PDWord(@v)^ shr 31) <> 0;
exponent := (PDWord(@v)^ shr 23) and $FF;
mantissa := PDWord(@v)^ and DWord($7FFFFF);
end;
{$endif}
{$ifdef FPC_HAS_TYPE_Double}
procedure splitFloat(const v: double; out sign: boolean; out exponent: integer; out mantissa: QWord); overload;
begin
sign := (PQWord(@v)^ shr 63) <> 0;
exponent := (PQWord(@v)^ shr 52) and $7FF;
mantissa := PQWord(@v)^ and QWord($000FFFFFFFFFFFFF);
end;
{$endif}
{$ifdef FPC_HAS_TYPE_Extended}