Mercurial > hg > CbC > CbC_gcc
view gcc/ada/libgnat/s-genbig.adb @ 145:1830386684a0
gcc-9.2.0
| author | anatofuz |
|---|---|
| date | Thu, 13 Feb 2020 11:34:05 +0900 |
| parents | |
| children |
line wrap: on
line source
------------------------------------------------------------------------------ -- -- -- GNAT COMPILER COMPONENTS -- -- -- -- S Y S T E M . G E N E R I C _ B I G N U M S -- -- -- -- B o d y -- -- -- -- Copyright (C) 2012-2019, Free Software Foundation, Inc. -- -- -- -- GNAT is free software; you can redistribute it and/or modify it under -- -- terms of the GNU General Public License as published by the Free Soft- -- -- ware Foundation; either version 3, or (at your option) any later ver- -- -- sion. GNAT is distributed in the hope that it will be useful, but WITH- -- -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -- -- or FITNESS FOR A PARTICULAR PURPOSE. -- -- -- -- As a special exception under Section 7 of GPL version 3, you are granted -- -- additional permissions described in the GCC Runtime Library Exception, -- -- version 3.1, as published by the Free Software Foundation. -- -- -- -- You should have received a copy of the GNU General Public License and -- -- a copy of the GCC Runtime Library Exception along with this program; -- -- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see -- -- <http://www.gnu.org/licenses/>. -- -- -- -- GNAT was originally developed by the GNAT team at New York University. -- -- Extensive contributions were provided by Ada Core Technologies Inc. -- -- -- ------------------------------------------------------------------------------ -- This package provides arbitrary precision signed integer arithmetic. with System; use System; with System.Secondary_Stack; use System.Secondary_Stack; with System.Storage_Elements; use System.Storage_Elements; package body System.Generic_Bignums is use Interfaces; -- So that operations on Unsigned_32/Unsigned_64 are available type DD is mod Base ** 2; -- Double length digit used for intermediate computations function MSD (X : DD) return SD is (SD (X / Base)); function LSD (X : DD) return SD is (SD (X mod Base)); -- Most significant and least significant digit of double digit value function "&" (X, Y : SD) return DD is (DD (X) * Base + DD (Y)); -- Compose double digit value from two single digit values subtype LLI is Long_Long_Integer; One_Data : constant Digit_Vector (1 .. 1) := (1 => 1); -- Constant one Zero_Data : constant Digit_Vector (1 .. 0) := (1 .. 0 => 0); -- Constant zero ----------------------- -- Local Subprograms -- ----------------------- function Add (X, Y : Digit_Vector; X_Neg : Boolean; Y_Neg : Boolean) return Bignum with Pre => X'First = 1 and then Y'First = 1; -- This procedure adds two signed numbers returning the Sum, it is used -- for both addition and subtraction. The value computed is X + Y, with -- X_Neg and Y_Neg giving the signs of the operands. function Allocate_Bignum (Len : Length) return Bignum with Post => Allocate_Bignum'Result.Len = Len; -- Allocate Bignum value of indicated length on secondary stack. On return -- the Neg and D fields are left uninitialized. type Compare_Result is (LT, EQ, GT); -- Indicates result of comparison in following call function Compare (X, Y : Digit_Vector; X_Neg, Y_Neg : Boolean) return Compare_Result with Pre => X'First = 1 and then Y'First = 1; -- Compare (X with sign X_Neg) with (Y with sign Y_Neg), and return the -- result of the signed comparison. procedure Div_Rem (X, Y : Bignum; Quotient : out Bignum; Remainder : out Bignum; Discard_Quotient : Boolean := False; Discard_Remainder : Boolean := False); -- Returns the Quotient and Remainder from dividing abs (X) by abs (Y). The -- values of X and Y are not modified. If Discard_Quotient is True, then -- Quotient is undefined on return, and if Discard_Remainder is True, then -- Remainder is undefined on return. Service routine for Big_Div/Rem/Mod. procedure Free_Bignum (X : Bignum) is null; -- Called to free a Bignum value used in intermediate computations. In -- this implementation using the secondary stack, it does nothing at all, -- because we rely on Mark/Release, but it may be of use for some -- alternative implementation. function Normalize (X : Digit_Vector; Neg : Boolean := False) return Bignum; -- Given a digit vector and sign, allocate and construct a Bignum value. -- Note that X may have leading zeroes which must be removed, and if the -- result is zero, the sign is forced positive. --------- -- Add -- --------- function Add (X, Y : Digit_Vector; X_Neg : Boolean; Y_Neg : Boolean) return Bignum is begin -- If signs are the same, we are doing an addition, it is convenient to -- ensure that the first operand is the longer of the two. if X_Neg = Y_Neg then if X'Last < Y'Last then return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg); -- Here signs are the same, and the first operand is the longer else pragma Assert (X_Neg = Y_Neg and then X'Last >= Y'Last); -- Do addition, putting result in Sum (allowing for carry) declare Sum : Digit_Vector (0 .. X'Last); RD : DD; begin RD := 0; for J in reverse 1 .. X'Last loop RD := RD + DD (X (J)); if J >= 1 + (X'Last - Y'Last) then RD := RD + DD (Y (J - (X'Last - Y'Last))); end if; Sum (J) := LSD (RD); RD := RD / Base; end loop; Sum (0) := SD (RD); return Normalize (Sum, X_Neg); end; end if; -- Signs are different so really this is a subtraction, we want to make -- sure that the largest magnitude operand is the first one, and then -- the result will have the sign of the first operand. else declare CR : constant Compare_Result := Compare (X, Y, False, False); begin if CR = EQ then return Normalize (Zero_Data); elsif CR = LT then return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg); else pragma Assert (X_Neg /= Y_Neg and then CR = GT); -- Do subtraction, putting result in Diff declare Diff : Digit_Vector (1 .. X'Length); RD : DD; begin RD := 0; for J in reverse 1 .. X'Last loop RD := RD + DD (X (J)); if J >= 1 + (X'Last - Y'Last) then RD := RD - DD (Y (J - (X'Last - Y'Last))); end if; Diff (J) := LSD (RD); RD := (if RD < Base then 0 else -1); end loop; return Normalize (Diff, X_Neg); end; end if; end; end if; end Add; --------------------- -- Allocate_Bignum -- --------------------- function Allocate_Bignum (Len : Length) return Bignum is Addr : Address; begin -- Allocation on the heap if not Use_Secondary_Stack then declare B : Bignum; begin B := new Bignum_Data'(Len, False, (others => 0)); return B; end; -- Allocation on the secondary stack else -- Note: The approach used here is designed to avoid strict aliasing -- warnings that appeared previously using unchecked conversion. SS_Allocate (Addr, Storage_Offset (4 + 4 * Len)); declare B : Bignum; for B'Address use Addr'Address; pragma Import (Ada, B); BD : Bignum_Data (Len); for BD'Address use Addr; pragma Import (Ada, BD); -- Expose a writable view of discriminant BD.Len so that we can -- initialize it. We need to use the exact layout of the record -- to ensure that the Length field has 24 bits as expected. type Bignum_Data_Header is record Len : Length; Neg : Boolean; end record; for Bignum_Data_Header use record Len at 0 range 0 .. 23; Neg at 3 range 0 .. 7; end record; BDH : Bignum_Data_Header; for BDH'Address use BD'Address; pragma Import (Ada, BDH); pragma Assert (BDH.Len'Size = BD.Len'Size); begin BDH.Len := Len; return B; end; end if; end Allocate_Bignum; ------------- -- Big_Abs -- ------------- function Big_Abs (X : Bignum) return Bignum is begin return Normalize (X.D); end Big_Abs; ------------- -- Big_Add -- ------------- function Big_Add (X, Y : Bignum) return Bignum is begin return Add (X.D, Y.D, X.Neg, Y.Neg); end Big_Add; ------------- -- Big_Div -- ------------- -- This table is excerpted from RM 4.5.5(28-30) and shows how the result -- varies with the signs of the operands. -- A B A/B A B A/B -- -- 10 5 2 -10 5 -2 -- 11 5 2 -11 5 -2 -- 12 5 2 -12 5 -2 -- 13 5 2 -13 5 -2 -- 14 5 2 -14 5 -2 -- -- A B A/B A B A/B -- -- 10 -5 -2 -10 -5 2 -- 11 -5 -2 -11 -5 2 -- 12 -5 -2 -12 -5 2 -- 13 -5 -2 -13 -5 2 -- 14 -5 -2 -14 -5 2 function Big_Div (X, Y : Bignum) return Bignum is Q, R : Bignum; begin Div_Rem (X, Y, Q, R, Discard_Remainder => True); Q.Neg := Q.Len > 0 and then (X.Neg xor Y.Neg); return Q; end Big_Div; ------------- -- Big_Exp -- ------------- function Big_Exp (X, Y : Bignum) return Bignum is function "**" (X : Bignum; Y : SD) return Bignum; -- Internal routine where we know right operand is one word ---------- -- "**" -- ---------- function "**" (X : Bignum; Y : SD) return Bignum is begin case Y is -- X ** 0 is 1 when 0 => return Normalize (One_Data); -- X ** 1 is X when 1 => return Normalize (X.D); -- X ** 2 is X * X when 2 => return Big_Mul (X, X); -- For X greater than 2, use the recursion -- X even, X ** Y = (X ** (Y/2)) ** 2; -- X odd, X ** Y = (X ** (Y/2)) ** 2 * X; when others => declare XY2 : constant Bignum := X ** (Y / 2); XY2S : constant Bignum := Big_Mul (XY2, XY2); Res : Bignum; begin Free_Bignum (XY2); -- Raise storage error if intermediate value is getting too -- large, which we arbitrarily define as 200 words for now. if XY2S.Len > 200 then Free_Bignum (XY2S); raise Storage_Error with "exponentiation result is too large"; end if; -- Otherwise take care of even/odd cases if (Y and 1) = 0 then return XY2S; else Res := Big_Mul (XY2S, X); Free_Bignum (XY2S); return Res; end if; end; end case; end "**"; -- Start of processing for Big_Exp begin -- Error if right operand negative if Y.Neg then raise Constraint_Error with "exponentiation to negative power"; -- X ** 0 is always 1 (including 0 ** 0, so do this test first) elsif Y.Len = 0 then return Normalize (One_Data); -- 0 ** X is always 0 (for X non-zero) elsif X.Len = 0 then return Normalize (Zero_Data); -- (+1) ** Y = 1 -- (-1) ** Y = +/-1 depending on whether Y is even or odd elsif X.Len = 1 and then X.D (1) = 1 then return Normalize (X.D, Neg => X.Neg and then ((Y.D (Y.Len) and 1) = 1)); -- If the absolute value of the base is greater than 1, then the -- exponent must not be bigger than one word, otherwise the result -- is ludicrously large, and we just signal Storage_Error right away. elsif Y.Len > 1 then raise Storage_Error with "exponentiation result is too large"; -- Special case (+/-)2 ** K, where K is 1 .. 31 using a shift elsif X.Len = 1 and then X.D (1) = 2 and then Y.D (1) < 32 then declare D : constant Digit_Vector (1 .. 1) := (1 => Shift_Left (SD'(1), Natural (Y.D (1)))); begin return Normalize (D, X.Neg); end; -- Remaining cases have right operand of one word else return X ** Y.D (1); end if; end Big_Exp; ------------ -- Big_EQ -- ------------ function Big_EQ (X, Y : Bignum) return Boolean is begin return Compare (X.D, Y.D, X.Neg, Y.Neg) = EQ; end Big_EQ; ------------ -- Big_GE -- ------------ function Big_GE (X, Y : Bignum) return Boolean is begin return Compare (X.D, Y.D, X.Neg, Y.Neg) /= LT; end Big_GE; ------------ -- Big_GT -- ------------ function Big_GT (X, Y : Bignum) return Boolean is begin return Compare (X.D, Y.D, X.Neg, Y.Neg) = GT; end Big_GT; ------------ -- Big_LE -- ------------ function Big_LE (X, Y : Bignum) return Boolean is begin return Compare (X.D, Y.D, X.Neg, Y.Neg) /= GT; end Big_LE; ------------ -- Big_LT -- ------------ function Big_LT (X, Y : Bignum) return Boolean is begin return Compare (X.D, Y.D, X.Neg, Y.Neg) = LT; end Big_LT; ------------- -- Big_Mod -- ------------- -- This table is excerpted from RM 4.5.5(28-30) and shows how the result -- of Rem and Mod vary with the signs of the operands. -- A B A mod B A rem B A B A mod B A rem B -- 10 5 0 0 -10 5 0 0 -- 11 5 1 1 -11 5 4 -1 -- 12 5 2 2 -12 5 3 -2 -- 13 5 3 3 -13 5 2 -3 -- 14 5 4 4 -14 5 1 -4 -- A B A mod B A rem B A B A mod B A rem B -- 10 -5 0 0 -10 -5 0 0 -- 11 -5 -4 1 -11 -5 -1 -1 -- 12 -5 -3 2 -12 -5 -2 -2 -- 13 -5 -2 3 -13 -5 -3 -3 -- 14 -5 -1 4 -14 -5 -4 -4 function Big_Mod (X, Y : Bignum) return Bignum is Q, R : Bignum; begin -- If signs are same, result is same as Rem if X.Neg = Y.Neg then return Big_Rem (X, Y); -- Case where Mod is different else -- Do division Div_Rem (X, Y, Q, R, Discard_Quotient => True); -- Zero result is unchanged if R.Len = 0 then return R; -- Otherwise adjust result else declare T1 : constant Bignum := Big_Sub (Y, R); begin T1.Neg := Y.Neg; Free_Bignum (R); return T1; end; end if; end if; end Big_Mod; ------------- -- Big_Mul -- ------------- function Big_Mul (X, Y : Bignum) return Bignum is Result : Digit_Vector (1 .. X.Len + Y.Len) := (others => 0); -- Accumulate result (max length of result is sum of operand lengths) L : Length; -- Current result digit D : DD; -- Result digit begin for J in 1 .. X.Len loop for K in 1 .. Y.Len loop L := Result'Last - (X.Len - J) - (Y.Len - K); D := DD (X.D (J)) * DD (Y.D (K)) + DD (Result (L)); Result (L) := LSD (D); D := D / Base; -- D is carry which must be propagated while D /= 0 and then L >= 1 loop L := L - 1; D := D + DD (Result (L)); Result (L) := LSD (D); D := D / Base; end loop; -- Must not have a carry trying to extend max length pragma Assert (D = 0); end loop; end loop; -- Return result return Normalize (Result, X.Neg xor Y.Neg); end Big_Mul; ------------ -- Big_NE -- ------------ function Big_NE (X, Y : Bignum) return Boolean is begin return Compare (X.D, Y.D, X.Neg, Y.Neg) /= EQ; end Big_NE; ------------- -- Big_Neg -- ------------- function Big_Neg (X : Bignum) return Bignum is begin return Normalize (X.D, not X.Neg); end Big_Neg; ------------- -- Big_Rem -- ------------- -- This table is excerpted from RM 4.5.5(28-30) and shows how the result -- varies with the signs of the operands. -- A B A rem B A B A rem B -- 10 5 0 -10 5 0 -- 11 5 1 -11 5 -1 -- 12 5 2 -12 5 -2 -- 13 5 3 -13 5 -3 -- 14 5 4 -14 5 -4 -- A B A rem B A B A rem B -- 10 -5 0 -10 -5 0 -- 11 -5 1 -11 -5 -1 -- 12 -5 2 -12 -5 -2 -- 13 -5 3 -13 -5 -3 -- 14 -5 4 -14 -5 -4 function Big_Rem (X, Y : Bignum) return Bignum is Q, R : Bignum; begin Div_Rem (X, Y, Q, R, Discard_Quotient => True); R.Neg := R.Len > 0 and then X.Neg; return R; end Big_Rem; ------------- -- Big_Sub -- ------------- function Big_Sub (X, Y : Bignum) return Bignum is begin -- If right operand zero, return left operand (avoiding sharing) if Y.Len = 0 then return Normalize (X.D, X.Neg); -- Otherwise add negative of right operand else return Add (X.D, Y.D, X.Neg, not Y.Neg); end if; end Big_Sub; ------------- -- Compare -- ------------- function Compare (X, Y : Digit_Vector; X_Neg, Y_Neg : Boolean) return Compare_Result is begin -- Signs are different, that's decisive, since 0 is always plus if X_Neg /= Y_Neg then return (if X_Neg then LT else GT); -- Lengths are different, that's decisive since no leading zeroes elsif X'Last /= Y'Last then return (if (X'Last > Y'Last) xor X_Neg then GT else LT); -- Need to compare data else for J in X'Range loop if X (J) /= Y (J) then return (if (X (J) > Y (J)) xor X_Neg then GT else LT); end if; end loop; return EQ; end if; end Compare; ------------- -- Div_Rem -- ------------- procedure Div_Rem (X, Y : Bignum; Quotient : out Bignum; Remainder : out Bignum; Discard_Quotient : Boolean := False; Discard_Remainder : Boolean := False) is begin -- Error if division by zero if Y.Len = 0 then raise Constraint_Error with "division by zero"; end if; -- Handle simple cases with special tests -- If X < Y then quotient is zero and remainder is X if Compare (X.D, Y.D, False, False) = LT then Remainder := Normalize (X.D); Quotient := Normalize (Zero_Data); return; -- If both X and Y are less than 2**63-1, we can use Long_Long_Integer -- arithmetic. Note it is good not to do an accurate range check against -- Long_Long_Integer since -2**63 / -1 overflows. elsif (X.Len <= 1 or else (X.Len = 2 and then X.D (1) < 2**31)) and then (Y.Len <= 1 or else (Y.Len = 2 and then Y.D (1) < 2**31)) then declare A : constant LLI := abs (From_Bignum (X)); B : constant LLI := abs (From_Bignum (Y)); begin Quotient := To_Bignum (A / B); Remainder := To_Bignum (A rem B); return; end; -- Easy case if divisor is one digit elsif Y.Len = 1 then declare ND : DD; Div : constant DD := DD (Y.D (1)); Result : Digit_Vector (1 .. X.Len); Remdr : Digit_Vector (1 .. 1); begin ND := 0; for J in 1 .. X.Len loop ND := Base * ND + DD (X.D (J)); Result (J) := SD (ND / Div); ND := ND rem Div; end loop; Quotient := Normalize (Result); Remdr (1) := SD (ND); Remainder := Normalize (Remdr); return; end; end if; -- The complex full multi-precision case. We will employ algorithm -- D defined in the section "The Classical Algorithms" (sec. 4.3.1) -- of Donald Knuth's "The Art of Computer Programming", Vol. 2, 2nd -- edition. The terminology is adjusted for this section to match that -- reference. -- We are dividing X.Len digits of X (called u here) by Y.Len digits -- of Y (called v here), developing the quotient and remainder. The -- numbers are represented using Base, which was chosen so that we have -- the operations of multiplying to single digits (SD) to form a double -- digit (DD), and dividing a double digit (DD) by a single digit (SD) -- to give a single digit quotient and a single digit remainder. -- Algorithm D from Knuth -- Comments here with square brackets are directly from Knuth Algorithm_D : declare -- The following lower case variables correspond exactly to the -- terminology used in algorithm D. m : constant Length := X.Len - Y.Len; n : constant Length := Y.Len; b : constant DD := Base; u : Digit_Vector (0 .. m + n); v : Digit_Vector (1 .. n); q : Digit_Vector (0 .. m); r : Digit_Vector (1 .. n); u0 : SD renames u (0); v1 : SD renames v (1); v2 : SD renames v (2); d : DD; j : Length; qhat : DD; rhat : DD; temp : DD; begin -- Initialize data of left and right operands for J in 1 .. m + n loop u (J) := X.D (J); end loop; for J in 1 .. n loop v (J) := Y.D (J); end loop; -- [Division of nonnegative integers.] Given nonnegative integers u -- = (ul,u2..um+n) and v = (v1,v2..vn), where v1 /= 0 and n > 1, we -- form the quotient u / v = (q0,ql..qm) and the remainder u mod v = -- (r1,r2..rn). pragma Assert (v1 /= 0); pragma Assert (n > 1); -- Dl. [Normalize.] Set d = b/(vl + 1). Then set (u0,u1,u2..um+n) -- equal to (u1,u2..um+n) times d, and set (v1,v2..vn) equal to -- (v1,v2..vn) times d. Note the introduction of a new digit position -- u0 at the left of u1; if d = 1 all we need to do in this step is -- to set u0 = 0. d := b / (DD (v1) + 1); if d = 1 then u0 := 0; else declare Carry : DD; Tmp : DD; begin -- Multiply Dividend (u) by d Carry := 0; for J in reverse 1 .. m + n loop Tmp := DD (u (J)) * d + Carry; u (J) := LSD (Tmp); Carry := Tmp / Base; end loop; u0 := SD (Carry); -- Multiply Divisor (v) by d Carry := 0; for J in reverse 1 .. n loop Tmp := DD (v (J)) * d + Carry; v (J) := LSD (Tmp); Carry := Tmp / Base; end loop; pragma Assert (Carry = 0); end; end if; -- D2. [Initialize j.] Set j = 0. The loop on j, steps D2 through D7, -- will be essentially a division of (uj, uj+1..uj+n) by (v1,v2..vn) -- to get a single quotient digit qj. j := 0; -- Loop through digits loop -- Note: In the original printing, step D3 was as follows: -- D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise -- set qhat to (uj,uj+1)/v1. Now test if v2 * qhat is greater than -- (uj*b + uj+1 - qhat*v1)*b + uj+2. If so, decrease qhat by 1 and -- repeat this test -- This had a bug not discovered till 1995, see Vol 2 errata: -- http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz. Under -- rare circumstances the expression in the test could overflow. -- This version was further corrected in 2005, see Vol 2 errata: -- http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz. -- The code below is the fixed version of this step. -- D3. [Calculate qhat.] Set qhat to (uj,uj+1)/v1 and rhat to -- to (uj,uj+1) mod v1. temp := u (j) & u (j + 1); qhat := temp / DD (v1); rhat := temp mod DD (v1); -- D3 (continued). Now test if qhat >= b or v2*qhat > (rhat,uj+2): -- if so, decrease qhat by 1, increase rhat by v1, and repeat this -- test if rhat < b. [The test on v2 determines at high speed -- most of the cases in which the trial value qhat is one too -- large, and eliminates all cases where qhat is two too large.] while qhat >= b or else DD (v2) * qhat > LSD (rhat) & u (j + 2) loop qhat := qhat - 1; rhat := rhat + DD (v1); exit when rhat >= b; end loop; -- D4. [Multiply and subtract.] Replace (uj,uj+1..uj+n) by -- (uj,uj+1..uj+n) minus qhat times (v1,v2..vn). This step -- consists of a simple multiplication by a one-place number, -- combined with a subtraction. -- The digits (uj,uj+1..uj+n) are always kept positive; if the -- result of this step is actually negative then (uj,uj+1..uj+n) -- is left as the true value plus b**(n+1), i.e. as the b's -- complement of the true value, and a "borrow" to the left is -- remembered. declare Borrow : SD; Carry : DD; Temp : DD; Negative : Boolean; -- Records if subtraction causes a negative result, requiring -- an add back (case where qhat turned out to be 1 too large). begin Borrow := 0; for K in reverse 1 .. n loop Temp := qhat * DD (v (K)) + DD (Borrow); Borrow := MSD (Temp); if LSD (Temp) > u (j + K) then Borrow := Borrow + 1; end if; u (j + K) := u (j + K) - LSD (Temp); end loop; Negative := u (j) < Borrow; u (j) := u (j) - Borrow; -- D5. [Test remainder.] Set qj = qhat. If the result of step -- D4 was negative, we will do the add back step (step D6). q (j) := LSD (qhat); if Negative then -- D6. [Add back.] Decrease qj by 1, and add (0,v1,v2..vn) -- to (uj,uj+1,uj+2..uj+n). (A carry will occur to the left -- of uj, and it is be ignored since it cancels with the -- borrow that occurred in D4.) q (j) := q (j) - 1; Carry := 0; for K in reverse 1 .. n loop Temp := DD (v (K)) + DD (u (j + K)) + Carry; u (j + K) := LSD (Temp); Carry := Temp / Base; end loop; u (j) := u (j) + SD (Carry); end if; end; -- D7. [Loop on j.] Increase j by one. Now if j <= m, go back to -- D3 (the start of the loop on j). j := j + 1; exit when not (j <= m); end loop; -- D8. [Unnormalize.] Now (qo,ql..qm) is the desired quotient, and -- the desired remainder may be obtained by dividing (um+1..um+n) -- by d. if not Discard_Quotient then Quotient := Normalize (q); end if; if not Discard_Remainder then declare Remdr : DD; begin Remdr := 0; for K in 1 .. n loop Remdr := Base * Remdr + DD (u (m + K)); r (K) := SD (Remdr / d); Remdr := Remdr rem d; end loop; pragma Assert (Remdr = 0); end; Remainder := Normalize (r); end if; end Algorithm_D; end Div_Rem; ----------------- -- From_Bignum -- ----------------- function From_Bignum (X : Bignum) return Long_Long_Integer is begin if X.Len = 0 then return 0; elsif X.Len = 1 then return (if X.Neg then -LLI (X.D (1)) else LLI (X.D (1))); elsif X.Len = 2 then declare Mag : constant DD := X.D (1) & X.D (2); begin if X.Neg and then Mag <= 2 ** 63 then return -LLI (Mag); elsif Mag < 2 ** 63 then return LLI (Mag); end if; end; end if; raise Constraint_Error with "expression value out of range"; end From_Bignum; ------------------------- -- Bignum_In_LLI_Range -- ------------------------- function Bignum_In_LLI_Range (X : Bignum) return Boolean is begin -- If length is 0 or 1, definitely fits if X.Len <= 1 then return True; -- If length is greater than 2, definitely does not fit elsif X.Len > 2 then return False; -- Length is 2, more tests needed else declare Mag : constant DD := X.D (1) & X.D (2); begin return Mag < 2 ** 63 or else (X.Neg and then Mag = 2 ** 63); end; end if; end Bignum_In_LLI_Range; --------------- -- Normalize -- --------------- function Normalize (X : Digit_Vector; Neg : Boolean := False) return Bignum is B : Bignum; J : Length; begin J := X'First; while J <= X'Last and then X (J) = 0 loop J := J + 1; end loop; B := Allocate_Bignum (X'Last - J + 1); B.Neg := B.Len > 0 and then Neg; B.D := X (J .. X'Last); return B; end Normalize; --------------- -- To_Bignum -- --------------- function To_Bignum (X : Long_Long_Integer) return Bignum is R : Bignum; begin if X = 0 then R := Allocate_Bignum (0); -- One word result elsif X in -(2 ** 32 - 1) .. +(2 ** 32 - 1) then R := Allocate_Bignum (1); R.D (1) := SD (abs (X)); -- Largest negative number annoyance elsif X = Long_Long_Integer'First then R := Allocate_Bignum (2); R.D (1) := 2 ** 31; R.D (2) := 0; -- Normal two word case else R := Allocate_Bignum (2); R.D (2) := SD (abs (X) mod Base); R.D (1) := SD (abs (X) / Base); end if; R.Neg := X < 0; return R; end To_Bignum; function To_Bignum (X : Unsigned_64) return Bignum is R : Bignum; begin if X = 0 then R := Allocate_Bignum (0); -- One word result elsif X < 2 ** 32 then R := Allocate_Bignum (1); R.D (1) := SD (X); -- Two word result else R := Allocate_Bignum (2); R.D (2) := SD (X mod Base); R.D (1) := SD (X / Base); end if; R.Neg := False; return R; end To_Bignum; ------------- -- Is_Zero -- ------------- function Is_Zero (X : Bignum) return Boolean is (X /= null and then X.D = Zero_Data); end System.Generic_Bignums;