1 ------------------------------------------------------------------------------
3 -- GNAT COMPILER COMPONENTS --
5 -- S Y S T E M . G E N E R I C _ B I G N U M S --
9 -- Copyright (C) 2012-2024, Free Software Foundation, Inc. --
11 -- GNAT is free software; you can redistribute it and/or modify it under --
12 -- terms of the GNU General Public License as published by the Free Soft- --
13 -- ware Foundation; either version 3, or (at your option) any later ver- --
14 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
15 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
16 -- or FITNESS FOR A PARTICULAR PURPOSE. --
18 -- As a special exception under Section 7 of GPL version 3, you are granted --
19 -- additional permissions described in the GCC Runtime Library Exception, --
20 -- version 3.1, as published by the Free Software Foundation. --
22 -- You should have received a copy of the GNU General Public License and --
23 -- a copy of the GCC Runtime Library Exception along with this program; --
24 -- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
25 -- <http://www.gnu.org/licenses/>. --
27 -- GNAT was originally developed by the GNAT team at New York University. --
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
30 ------------------------------------------------------------------------------
32 -- This package provides arbitrary precision signed integer arithmetic.
34 package body System
.Generic_Bignums
is
37 -- So that operations on Unsigned_32/Unsigned_64 are available
41 type DD
is mod Base
** 2;
42 -- Double length digit used for intermediate computations
44 function MSD
(X
: DD
) return SD
is (SD
(X
/ Base
));
45 function LSD
(X
: DD
) return SD
is (SD
(X
mod Base
));
46 -- Most significant and least significant digit of double digit value
48 function "&" (X
, Y
: SD
) return DD
is (DD
(X
) * Base
+ DD
(Y
));
49 -- Compose double digit value from two single digit values
51 subtype LLI
is Long_Long_Integer;
52 subtype LLLI
is Long_Long_Long_Integer
;
54 LLLI_Is_128
: constant Boolean := Long_Long_Long_Integer
'Size = 128;
55 -- True if Long_Long_Long_Integer is 128-bit large
57 One_Data
: constant Digit_Vector
(1 .. 1) := [1];
60 Zero_Data
: constant Digit_Vector
(1 .. 0) := [];
63 -----------------------
64 -- Local Subprograms --
65 -----------------------
70 Y_Neg
: Boolean) return Big_Integer
72 Pre
=> X
'First = 1 and then Y
'First = 1;
73 -- This procedure adds two signed numbers returning the Sum, it is used
74 -- for both addition and subtraction. The value computed is X + Y, with
75 -- X_Neg and Y_Neg giving the signs of the operands.
77 type Compare_Result
is (LT
, EQ
, GT
);
78 -- Indicates result of comparison in following call
82 X_Neg
, Y_Neg
: Boolean) return Compare_Result
84 Pre
=> X
'First = 1 and then Y
'First = 1;
85 -- Compare (X with sign X_Neg) with (Y with sign Y_Neg), and return the
86 -- result of the signed comparison.
90 Quotient
: out Big_Integer
;
91 Remainder
: out Big_Integer
;
92 Discard_Quotient
: Boolean := False;
93 Discard_Remainder
: Boolean := False);
94 -- Returns the Quotient and Remainder from dividing abs (X) by abs (Y). The
95 -- values of X and Y are not modified. If Discard_Quotient is True, then
96 -- Quotient is undefined on return, and if Discard_Remainder is True, then
97 -- Remainder is undefined on return. Service routine for Big_Div/Rem/Mod.
101 Neg
: Boolean := False) return Big_Integer
;
102 -- Given a digit vector and sign, allocate and construct a big integer
103 -- value. Note that X may have leading zeroes which must be removed, and if
104 -- the result is zero, the sign is forced positive.
105 -- If X is too big, Storage_Error is raised.
107 function "**" (X
: Bignum
; Y
: SD
) return Big_Integer
;
108 -- Exponentiation routine where we know right operand is one word
115 (X
, Y
: Digit_Vector
;
117 Y_Neg
: Boolean) return Big_Integer
120 -- If signs are the same, we are doing an addition, it is convenient to
121 -- ensure that the first operand is the longer of the two.
123 if X_Neg
= Y_Neg
then
124 if X
'Last < Y
'Last then
125 return Add
(X
=> Y
, Y
=> X
, X_Neg
=> Y_Neg
, Y_Neg
=> X_Neg
);
127 -- Here signs are the same, and the first operand is the longer
130 pragma Assert
(X_Neg
= Y_Neg
and then X
'Last >= Y
'Last);
132 -- Do addition, putting result in Sum (allowing for carry)
135 Sum
: Digit_Vector
(0 .. X
'Last);
140 for J
in reverse 1 .. X
'Last loop
141 RD
:= RD
+ DD
(X
(J
));
143 if J
>= 1 + (X
'Last - Y
'Last) then
144 RD
:= RD
+ DD
(Y
(J
- (X
'Last - Y
'Last)));
152 return Normalize
(Sum
, X_Neg
);
156 -- Signs are different so really this is a subtraction, we want to make
157 -- sure that the largest magnitude operand is the first one, and then
158 -- the result will have the sign of the first operand.
162 CR
: constant Compare_Result
:= Compare
(X
, Y
, False, False);
166 return Normalize
(Zero_Data
);
169 return Add
(X
=> Y
, Y
=> X
, X_Neg
=> Y_Neg
, Y_Neg
=> X_Neg
);
172 pragma Assert
(X_Neg
/= Y_Neg
and then CR
= GT
);
174 -- Do subtraction, putting result in Diff
177 Diff
: Digit_Vector
(1 .. X
'Length);
182 for J
in reverse 1 .. X
'Last loop
183 RD
:= RD
+ DD
(X
(J
));
185 if J
>= 1 + (X
'Last - Y
'Last) then
186 RD
:= RD
- DD
(Y
(J
- (X
'Last - Y
'Last)));
189 Diff
(J
) := LSD
(RD
);
190 RD
:= (if RD
< Base
then 0 else -1);
193 return Normalize
(Diff
, X_Neg
);
204 function Big_Abs
(X
: Bignum
) return Big_Integer
is
206 return Normalize
(X
.D
);
213 function Big_Add
(X
, Y
: Bignum
) return Big_Integer
is
215 return Add
(X
.D
, Y
.D
, X
.Neg
, Y
.Neg
);
222 -- This table is excerpted from RM 4.5.5(28-30) and shows how the result
223 -- varies with the signs of the operands.
241 function Big_Div
(X
, Y
: Bignum
) return Big_Integer
is
242 Q
, R
: aliased Big_Integer
;
244 Div_Rem
(X
, Y
, Q
, R
, Discard_Remainder
=> True);
245 To_Bignum
(Q
).Neg
:= To_Bignum
(Q
).Len
> 0 and then (X
.Neg
xor Y
.Neg
);
253 function "**" (X
: Bignum
; Y
: SD
) return Big_Integer
is
260 return Normalize
(One_Data
);
265 return Normalize
(X
.D
);
270 return Big_Mul
(X
, X
);
272 -- For X greater than 2, use the recursion
274 -- X even, X ** Y = (X ** (Y/2)) ** 2;
275 -- X odd, X ** Y = (X ** (Y/2)) ** 2 * X;
279 XY2
: aliased Big_Integer
:= X
** (Y
/ 2);
280 XY2S
: aliased Big_Integer
:=
281 Big_Mul
(To_Bignum
(XY2
), To_Bignum
(XY2
));
284 Free_Big_Integer
(XY2
);
286 if (Y
and 1) = 0 then
289 return Res
: constant Big_Integer
:=
290 Big_Mul
(To_Bignum
(XY2S
), X
)
292 Free_Big_Integer
(XY2S
);
303 function Big_Exp
(X
, Y
: Bignum
) return Big_Integer
is
305 -- Error if right operand negative
308 raise Constraint_Error
with "exponentiation to negative power";
310 -- X ** 0 is always 1 (including 0 ** 0, so do this test first)
313 return Normalize
(One_Data
);
315 -- 0 ** X is always 0 (for X non-zero)
318 return Normalize
(Zero_Data
);
321 -- (-1) ** Y = +/-1 depending on whether Y is even or odd
323 elsif X
.Len
= 1 and then X
.D
(1) = 1 then
325 (X
.D
, Neg
=> X
.Neg
and then (Y
.D
(Y
.Len
) and 1) = 1);
327 -- If the absolute value of the base is greater than 1, then the
328 -- exponent must not be bigger than one word, otherwise the result
329 -- is ludicrously large, and we just signal Storage_Error right away.
332 raise Storage_Error
with "exponentiation result is too large";
334 -- Special case (+/-)2 ** K, where K is 1 .. 31 using a shift
336 elsif X
.Len
= 1 and then X
.D
(1) = 2 and then Y
.D
(1) < 32 then
338 D
: constant Digit_Vector
(1 .. 1) :=
339 [Shift_Left
(SD
'(1), Natural (Y.D (1)))];
341 return Normalize (D, X.Neg);
344 -- Remaining cases have right operand of one word
355 function Big_And (X, Y : Bignum) return Big_Integer is
357 if X.Len > Y.Len then
358 return Big_And (X => Y, Y => X);
361 -- X is the smallest integer
364 Result : Digit_Vector (1 .. X.Len);
365 Diff : constant Length := Y.Len - X.Len;
367 for J in 1 .. X.Len loop
368 Result (J) := X.D (J) and Y.D (J + Diff);
371 return Normalize (Result, X.Neg and Y.Neg);
379 function Big_Or (X, Y : Bignum) return Big_Integer is
381 if X.Len < Y.Len then
382 return Big_Or (X => Y, Y => X);
385 -- X is the largest integer
388 Result : Digit_Vector (1 .. X.Len);
390 Diff : constant Length := X.Len - Y.Len;
395 while Index <= Diff loop
396 Result (Index) := X.D (Index);
400 for J in 1 .. Y.Len loop
401 Result (Index) := X.D (Index) or Y.D (J);
405 return Normalize (Result, X.Neg or Y.Neg);
413 function Big_Shift_Left (X : Bignum; Amount : Natural) return Big_Integer is
416 raise Constraint_Error;
417 elsif Amount = 0 then
418 return Allocate_Big_Integer (X.D, False);
422 Shift : constant Natural := Amount rem SD'Size;
423 Result : Digit_Vector (0 .. X.Len + Amount / SD'Size);
427 for J in X.Len + 1 .. Result'Last loop
431 for J in reverse 1 .. X.Len loop
432 Result (J) := Shift_Left (X.D (J), Shift) or Carry;
433 Carry := Shift_Right (X.D (J), SD'Size - Shift);
437 return Normalize (Result, False);
441 ---------------------
442 -- Big_Shift_Right --
443 ---------------------
445 function Big_Shift_Right
446 (X : Bignum; Amount : Natural) return Big_Integer is
449 raise Constraint_Error;
450 elsif Amount = 0 then
451 return Allocate_Big_Integer (X.D, False);
455 Shift : constant Natural := Amount rem SD'Size;
456 Result : Digit_Vector (1 .. X.Len - Amount / SD'Size);
460 for J in 1 .. Result'Last - 1 loop
461 Result (J) := Shift_Right (X.D (J), Shift) or Carry;
462 Carry := Shift_Left (X.D (J), SD'Size - Shift);
465 Result (Result'Last) :=
466 Shift_Right (X.D (Result'Last), Shift) or Carry;
468 return Normalize (Result, False);
476 function Big_EQ (X, Y : Bignum) return Boolean is
478 return Compare (X.D, Y.D, X.Neg, Y.Neg) = EQ;
485 function Big_GE (X, Y : Bignum) return Boolean is
487 return Compare (X.D, Y.D, X.Neg, Y.Neg) /= LT;
494 function Big_GT (X, Y : Bignum) return Boolean is
496 return Compare (X.D, Y.D, X.Neg, Y.Neg) = GT;
503 function Big_LE (X, Y : Bignum) return Boolean is
505 return Compare (X.D, Y.D, X.Neg, Y.Neg) /= GT;
512 function Big_LT (X, Y : Bignum) return Boolean is
514 return Compare (X.D, Y.D, X.Neg, Y.Neg) = LT;
521 -- This table is excerpted from RM 4.5.5(28-30) and shows how the result
522 -- of Rem and Mod vary with the signs of the operands.
524 -- A B A mod B A rem B A B A mod B A rem B
526 -- 10 5 0 0 -10 5 0 0
527 -- 11 5 1 1 -11 5 4 -1
528 -- 12 5 2 2 -12 5 3 -2
529 -- 13 5 3 3 -13 5 2 -3
530 -- 14 5 4 4 -14 5 1 -4
532 -- A B A mod B A rem B A B A mod B A rem B
534 -- 10 -5 0 0 -10 -5 0 0
535 -- 11 -5 -4 1 -11 -5 -1 -1
536 -- 12 -5 -3 2 -12 -5 -2 -2
537 -- 13 -5 -2 3 -13 -5 -3 -3
538 -- 14 -5 -1 4 -14 -5 -4 -4
540 function Big_Mod (X, Y : Bignum) return Big_Integer is
541 Q, R : aliased Big_Integer;
544 -- If signs are same, result is same as Rem
546 if X.Neg = Y.Neg then
547 return Big_Rem (X, Y);
549 -- Case where Mod is different
554 Div_Rem (X, Y, Q, R, Discard_Quotient => True);
556 -- Zero result is unchanged
558 if To_Bignum (R).Len = 0 then
561 -- Otherwise adjust result
565 T1 : aliased Big_Integer := Big_Sub (Y, To_Bignum (R));
567 To_Bignum (T1).Neg := Y.Neg;
568 Free_Big_Integer (R);
579 function Big_Mul (X, Y : Bignum) return Big_Integer is
580 Result : Digit_Vector (1 .. X.Len + Y.Len) := [others => 0];
581 -- Accumulate result (max length of result is sum of operand lengths)
584 -- Current result digit
590 for J in 1 .. X.Len loop
591 for K in 1 .. Y.Len loop
592 L := Result'Last - (X.Len - J) - (Y.Len - K);
593 D := DD (X.D (J)) * DD (Y.D (K)) + DD (Result (L));
594 Result (L) := LSD (D);
597 -- D is carry which must be propagated
599 while D /= 0 and then L >= 1 loop
601 D := D + DD (Result (L));
602 Result (L) := LSD (D);
606 -- Must not have a carry trying to extend max length
608 pragma Assert (D = 0);
614 return Normalize (Result, X.Neg xor Y.Neg);
621 function Big_NE (X, Y : Bignum) return Boolean is
623 return Compare (X.D, Y.D, X.Neg, Y.Neg) /= EQ;
630 function Big_Neg (X : Bignum) return Big_Integer is
632 return Normalize (X.D, not X.Neg);
639 -- This table is excerpted from RM 4.5.5(28-30) and shows how the result
640 -- varies with the signs of the operands.
642 -- A B A rem B A B A rem B
650 -- A B A rem B A B A rem B
658 function Big_Rem (X, Y : Bignum) return Big_Integer is
659 Q, R : aliased Big_Integer;
661 Div_Rem (X, Y, Q, R, Discard_Quotient => True);
662 To_Bignum (R).Neg := To_Bignum (R).Len > 0 and then X.Neg;
670 function Big_Sub (X, Y : Bignum) return Big_Integer is
672 -- If right operand zero, return left operand (avoiding sharing)
675 return Normalize (X.D, X.Neg);
677 -- Otherwise add negative of right operand
680 return Add (X.D, Y.D, X.Neg, not Y.Neg);
689 (X, Y : Digit_Vector;
690 X_Neg, Y_Neg : Boolean) return Compare_Result
693 -- Signs are different, that's decisive, since 0 is always plus
695 if X_Neg /= Y_Neg then
696 return (if X_Neg then LT else GT);
698 -- Lengths are different, that's decisive since no leading zeroes
700 elsif X'Last /= Y'Last then
701 return (if X'Last > Y'Last xor X_Neg then GT else LT);
703 -- Need to compare data
706 for J in X'Range loop
707 if X (J) /= Y (J) then
708 return (if X (J) > Y (J) xor X_Neg then GT else LT);
722 Quotient : out Big_Integer;
723 Remainder : out Big_Integer;
724 Discard_Quotient : Boolean := False;
725 Discard_Remainder : Boolean := False) is
727 -- Error if division by zero
730 raise Constraint_Error with "division by zero";
733 -- Handle simple cases with special tests
735 -- If X < Y then quotient is zero and remainder is X
737 if Compare (X.D, Y.D, False, False) = LT then
738 if not Discard_Quotient then
739 Quotient := Normalize (Zero_Data);
742 if not Discard_Remainder then
743 Remainder := Normalize (X.D);
748 -- If both X and Y are less than 2**63-1, we can use Long_Long_Integer
749 -- arithmetic. Note it is good not to do an accurate range check against
750 -- Long_Long_Integer since -2**63 / -1 overflows.
752 elsif (X.Len <= 1 or else (X.Len = 2 and then X.D (1) < 2**31))
754 (Y.Len <= 1 or else (Y.Len = 2 and then Y.D (1) < 2**31))
757 A : constant LLI := abs (From_Bignum (X));
758 B : constant LLI := abs (From_Bignum (Y));
760 if not Discard_Quotient then
761 Quotient := To_Bignum (A / B);
764 if not Discard_Remainder then
765 Remainder := To_Bignum (A rem B);
771 -- Easy case if divisor is one digit
776 Div : constant DD := DD (Y.D (1));
778 Result : Digit_Vector (1 .. X.Len);
779 Remdr : Digit_Vector (1 .. 1);
783 for J in 1 .. X.Len loop
784 ND := Base * ND + DD (X.D (J));
785 pragma Assert (Div /= 0);
786 Result (J) := SD (ND / Div);
790 if not Discard_Quotient then
791 Quotient := Normalize (Result);
794 if not Discard_Remainder then
795 Remdr (1) := SD (ND);
796 Remainder := Normalize (Remdr);
803 -- The complex full multi-precision case. We will employ algorithm
804 -- D defined in the section "The Classical Algorithms" (sec. 4.3.1)
805 -- of Donald Knuth's "The Art of Computer Programming", Vol. 2, 2nd
806 -- edition. The terminology is adjusted for this section to match that
809 -- We are dividing X.Len digits of X (called u here) by Y.Len digits
810 -- of Y (called v here), developing the quotient and remainder. The
811 -- numbers are represented using Base, which was chosen so that we have
812 -- the operations of multiplying to single digits (SD) to form a double
813 -- digit (DD), and dividing a double digit (DD) by a single digit (SD)
814 -- to give a single digit quotient and a single digit remainder.
816 -- Algorithm D from Knuth
818 -- Comments here with square brackets are directly from Knuth
820 Algorithm_D : declare
822 -- The following lower case variables correspond exactly to the
823 -- terminology used in algorithm D.
825 m : constant Length := X.Len - Y.Len;
826 n : constant Length := Y.Len;
827 b : constant DD := Base;
829 u : Digit_Vector (0 .. m + n);
830 v : Digit_Vector (1 .. n);
831 q : Digit_Vector (0 .. m);
832 r : Digit_Vector (1 .. n);
834 u0 : SD renames u (0);
835 v1 : SD renames v (1);
836 v2 : SD renames v (2);
845 -- Initialize data of left and right operands
847 for J in 1 .. m + n loop
855 -- [Division of nonnegative integers.] Given nonnegative integers u
856 -- = (ul,u2..um+n) and v = (v1,v2..vn), where v1 /= 0 and n > 1, we
857 -- form the quotient u / v = (q0,ql..qm) and the remainder u mod v =
860 pragma Assert (v1 /= 0);
861 pragma Assert (n > 1);
863 -- Dl. [Normalize.] Set d = b/(vl + 1). Then set (u0,u1,u2..um+n)
864 -- equal to (u1,u2..um+n) times d, and set (v1,v2..vn) equal to
865 -- (v1,v2..vn) times d. Note the introduction of a new digit position
866 -- u0 at the left of u1; if d = 1 all we need to do in this step is
869 d := b / (DD (v1) + 1);
880 -- Multiply Dividend (u) by d
883 for J in reverse 1 .. m + n loop
884 Tmp := DD (u (J)) * d + Carry;
891 -- Multiply Divisor (v) by d
894 for J in reverse 1 .. n loop
895 Tmp := DD (v (J)) * d + Carry;
900 pragma Assert (Carry = 0);
904 -- D2. [Initialize j.] Set j = 0. The loop on j, steps D2 through D7,
905 -- will be essentially a division of (uj, uj+1..uj+n) by (v1,v2..vn)
906 -- to get a single quotient digit qj.
910 -- Loop through digits
913 -- Note: In the original printing, step D3 was as follows:
915 -- D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise
916 -- set qhat to (uj,uj+1)/v1. Now test if v2 * qhat is greater than
917 -- (uj*b + uj+1 - qhat*v1)*b + uj+2. If so, decrease qhat by 1 and
920 -- This had a bug not discovered till 1995, see Vol 2 errata:
921 -- http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz. Under
922 -- rare circumstances the expression in the test could overflow.
923 -- This version was further corrected in 2005, see Vol 2 errata:
924 -- http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz.
925 -- The code below is the fixed version of this step.
927 -- D3. [Calculate qhat.] Set qhat to (uj,uj+1)/v1 and rhat to
928 -- to (uj,uj+1) mod v1.
930 temp := u (j) & u (j + 1);
931 qhat := temp / DD (v1);
932 rhat := temp mod DD (v1);
934 -- D3 (continued). Now test if qhat >= b or v2*qhat > (rhat,uj+2):
935 -- if so, decrease qhat by 1, increase rhat by v1, and repeat this
936 -- test if rhat < b. [The test on v2 determines at high speed
937 -- most of the cases in which the trial value qhat is one too
938 -- large, and eliminates all cases where qhat is two too large.]
941 or else DD (v2) * qhat > LSD (rhat) & u (j + 2)
944 rhat := rhat + DD (v1);
948 -- D4. [Multiply and subtract.] Replace (uj,uj+1..uj+n) by
949 -- (uj,uj+1..uj+n) minus qhat times (v1,v2..vn). This step
950 -- consists of a simple multiplication by a one-place number,
951 -- combined with a subtraction.
953 -- The digits (uj,uj+1..uj+n) are always kept positive; if the
954 -- result of this step is actually negative then (uj,uj+1..uj+n)
955 -- is left as the true value plus b**(n+1), i.e. as the b's
956 -- complement of the true value, and a "borrow" to the left is
965 -- Records if subtraction causes a negative result, requiring
966 -- an add back (case where qhat turned out to be 1 too large).
970 for K in reverse 1 .. n loop
971 Temp := qhat * DD (v (K)) + DD (Borrow);
972 Borrow := MSD (Temp);
974 if LSD (Temp) > u (j + K) then
975 Borrow := Borrow + 1;
978 u (j + K) := u (j + K) - LSD (Temp);
981 Negative := u (j) < Borrow;
982 u (j) := u (j) - Borrow;
984 -- D5. [Test remainder.] Set qj = qhat. If the result of step
985 -- D4 was negative, we will do the add back step (step D6).
991 -- D6. [Add back.] Decrease qj by 1, and add (0,v1,v2..vn)
992 -- to (uj,uj+1,uj+2..uj+n). (A carry will occur to the left
993 -- of uj, and it is be ignored since it cancels with the
994 -- borrow that occurred in D4.)
999 for K in reverse 1 .. n loop
1000 Temp := DD (v (K)) + DD (u (j + K)) + Carry;
1001 u (j + K) := LSD (Temp);
1002 Carry := Temp / Base;
1005 u (j) := u (j) + SD (Carry);
1009 -- D7. [Loop on j.] Increase j by one. Now if j <= m, go back to
1010 -- D3 (the start of the loop on j).
1013 exit when not (j <= m);
1016 -- D8. [Unnormalize.] Now (qo,ql..qm) is the desired quotient, and
1017 -- the desired remainder may be obtained by dividing (um+1..um+n)
1020 if not Discard_Quotient then
1021 Quotient := Normalize (q);
1024 if not Discard_Remainder then
1030 for K in 1 .. n loop
1031 Remdr := Base * Remdr + DD (u (m + K));
1032 r (K) := SD (Remdr / d);
1033 Remdr := Remdr rem d;
1036 pragma Assert (Remdr = 0);
1039 Remainder := Normalize (r);
1048 function From_Bignum (X : Bignum) return Long_Long_Long_Integer is
1053 elsif X.Len = 1 then
1054 return (if X.Neg then -LLLI (X.D (1)) else LLLI (X.D (1)));
1056 elsif X.Len = 2 then
1058 Mag : constant DD := X.D (1) & X.D (2);
1060 if X.Neg and then (Mag <= 2 ** 63 or else LLLI_Is_128) then
1062 elsif Mag < 2 ** 63 or else LLLI_Is_128 then
1067 elsif X.Len = 3 and then LLLI_Is_128 then
1069 Hi : constant SD := X.D (1);
1070 Lo : constant DD := X.D (2) & X.D (3);
1071 Mag : constant Unsigned_128 :=
1072 Shift_Left (Unsigned_128 (Hi), 64) + Unsigned_128 (Lo);
1074 return (if X.Neg then -LLLI (Mag) else LLLI (Mag));
1077 elsif X.Len = 4 and then LLLI_Is_128 then
1079 Hi : constant DD := X.D (1) & X.D (2);
1080 Lo : constant DD := X.D (3) & X.D (4);
1081 Mag : constant Unsigned_128 :=
1082 Shift_Left (Unsigned_128 (Hi), 64) + Unsigned_128 (Lo);
1085 and then (Hi < 2 ** 63 or else (Hi = 2 ** 63 and then Lo = 0))
1088 elsif Hi < 2 ** 63 then
1094 raise Constraint_Error with "expression value out of range";
1097 function From_Bignum (X : Bignum) return Long_Long_Integer is
1099 return Long_Long_Integer (Long_Long_Long_Integer'(From_Bignum
(X
)));
1102 function From_Bignum
(X
: Bignum
) return Unsigned_128
is
1107 elsif X
.Len
= 0 then
1110 elsif X
.Len
= 1 then
1111 return Unsigned_128
(X
.D
(1));
1113 elsif X
.Len
= 2 then
1114 return Unsigned_128
(DD
'(X.D (1) & X.D (2)));
1116 elsif X.Len = 3 and then LLLI_Is_128 then
1118 Shift_Left (Unsigned_128 (X.D (1)), 64) +
1119 Unsigned_128 (DD'(X
.D
(2) & X
.D
(3)));
1121 elsif X
.Len
= 4 and then LLLI_Is_128
then
1123 Shift_Left
(Unsigned_128
(DD
'(X.D (1) & X.D (2))), 64) +
1124 Unsigned_128 (DD'(X
.D
(3) & X
.D
(4)));
1127 raise Constraint_Error
with "expression value out of range";
1130 function From_Bignum
(X
: Bignum
) return Unsigned_64
is
1132 return Unsigned_64
(Unsigned_128
'(From_Bignum (X)));
1135 -------------------------
1136 -- Bignum_In_LLI_Range --
1137 -------------------------
1139 function Bignum_In_LLI_Range (X : Bignum) return Boolean is
1141 -- If length is 0 or 1, definitely fits
1146 -- If length is greater than 2, definitely does not fit
1148 elsif X.Len > 2 then
1151 -- Length is 2, more tests needed
1155 Mag : constant DD := X.D (1) & X.D (2);
1157 return Mag < 2 ** 63 or else (X.Neg and then Mag = 2 ** 63);
1160 end Bignum_In_LLI_Range;
1166 Bignum_Limit : constant := 200;
1170 Neg : Boolean := False) return Big_Integer
1176 while J <= X'Last and then X (J) = 0 loop
1180 if X'Last - J > Bignum_Limit then
1181 raise Storage_Error with "big integer limit exceeded";
1184 return Allocate_Big_Integer (X (J .. X'Last), J <= X'Last and then Neg);
1191 function To_Bignum (X : Long_Long_Long_Integer) return Big_Integer is
1193 function Convert_128
1194 (X : Long_Long_Long_Integer; Neg : Boolean) return Big_Integer;
1195 -- Convert a 128 bits natural integer to a Big_Integer
1201 function Convert_128
1202 (X : Long_Long_Long_Integer; Neg : Boolean) return Big_Integer
1204 Vector : Digit_Vector (1 .. 4);
1205 High : constant Unsigned_64 :=
1206 Unsigned_64 (Shift_Right (Unsigned_128 (X), 64));
1207 Low : constant Unsigned_64 :=
1208 Unsigned_64 (Unsigned_128 (X) and 16#FFFF_FFFF_FFFF_FFFF#);
1211 Vector (1) := SD (High / Base);
1212 Vector (2) := SD (High mod Base);
1213 Vector (3) := SD (Low / Base);
1214 Vector (4) := SD (Low mod Base);
1215 return Normalize (Vector, Neg);
1220 return Allocate_Big_Integer ([], False);
1224 elsif X in -(2 ** 32 - 1) .. +(2 ** 32 - 1) then
1225 return Allocate_Big_Integer ([SD (abs X)], X < 0);
1227 -- Large negative number annoyance
1229 elsif X = -2 ** 63 then
1230 return Allocate_Big_Integer ([2 ** 31, 0], True);
1232 elsif LLLI_Is_128 and then X = Long_Long_Long_Integer'First then
1233 return Allocate_Big_Integer ([2 ** 31, 0, 0, 0], True);
1235 -- Other negative numbers
1239 return Convert_128 (-X, True);
1241 return Allocate_Big_Integer
1242 ((SD ((-X) / Base), SD ((-X) mod Base)), True);
1249 return Convert_128 (X, False);
1251 return Allocate_Big_Integer
1252 ((SD (X / Base), SD (X mod Base)), False);
1257 function To_Bignum (X : Long_Long_Integer) return Big_Integer is
1259 return To_Bignum (Long_Long_Long_Integer (X));
1262 function To_Bignum (X : Unsigned_128) return Big_Integer is
1265 return Allocate_Big_Integer ([], False);
1269 elsif X < 2 ** 32 then
1270 return Allocate_Big_Integer ([SD (X)], False);
1274 elsif Shift_Right (X, 32) < 2 ** 32 then
1275 return Allocate_Big_Integer ([SD (X / Base), SD (X mod Base)], False);
1277 -- Three or four word result
1281 Vector : Digit_Vector (1 .. 4);
1282 High : constant Unsigned_64 := Unsigned_64 (Shift_Right (X, 64));
1283 Low : constant Unsigned_64 :=
1284 Unsigned_64 (X and 16#FFFF_FFFF_FFFF_FFFF#);
1287 Vector (1) := SD (High / Base);
1288 Vector (2) := SD (High mod Base);
1289 Vector (3) := SD (Low / Base);
1290 Vector (4) := SD (Low mod Base);
1291 return Normalize (Vector, False);
1296 function To_Bignum (X : Unsigned_64) return Big_Integer is
1298 return To_Bignum (Unsigned_128 (X));
1305 Hex_Chars : constant array (0 .. 15) of Character := "0123456789ABCDEF";
1308 (X : Bignum; Width : Natural := 0; Base : Positive := 10) return String
1310 Big_Base : aliased Bignum_Data := (1, False, [SD (Base)]);
1312 function Add_Base (S : String) return String;
1313 -- Add base information if Base /= 10
1315 function Leading_Padding
1317 Min_Length : Natural;
1318 Char : Character := ' ') return String;
1319 -- Return padding of Char concatenated with Str so that the resulting
1320 -- string is at least Min_Length long.
1322 function Image (Arg : Bignum) return String;
1323 -- Return image of Arg, assuming Arg is positive.
1325 function Image (N : Natural) return String;
1326 -- Return image of N, with no leading space.
1332 function Add_Base (S : String) return String is
1337 return Image (Base) & "#" & S & "#";
1345 function Image (N : Natural) return String is
1346 S : constant String := Natural'Image (N);
1348 return S (2 .. S'Last);
1351 function Image (Arg : Bignum) return String is
1353 if Big_LT (Arg, Big_Base'Unchecked_Access) then
1354 return [Hex_Chars (Natural (LLI'(From_Bignum
(Arg
))))];
1357 Div
: aliased Big_Integer
;
1358 Remain
: aliased Big_Integer
;
1362 Div_Rem
(Arg
, Big_Base
'Unchecked_Access, Div
, Remain
);
1363 R
:= Natural (LLI
'(From_Bignum (To_Bignum (Remain))));
1364 Free_Big_Integer (Remain);
1366 return S : constant String :=
1367 Image (To_Bignum (Div)) & Hex_Chars (R)
1369 Free_Big_Integer (Div);
1375 ---------------------
1376 -- Leading_Padding --
1377 ---------------------
1379 function Leading_Padding
1381 Min_Length : Natural;
1382 Char : Character := ' ') return String is
1384 return [1 .. Integer'Max (Integer (Min_Length) - Str'Length, 0)
1386 end Leading_Padding;
1388 Zero : aliased Bignum_Data := (0, False, D => Zero_Data);
1391 if Big_LT (X, Zero'Unchecked_Access) then
1393 X_Pos : aliased Bignum_Data := (X.Len, not X.Neg, X.D);
1395 return Leading_Padding
1396 ("-" & Add_Base (Image (X_Pos'Unchecked_Access)), Width);
1399 return Leading_Padding (" " & Add_Base (Image (X)), Width);
1407 function Is_Zero (X : Bignum) return Boolean is
1408 (X /= null and then X.D = Zero_Data);
1410 end System.Generic_Bignums;