aboutsummaryrefslogtreecommitdiffstats
path: root/gcc-4.9/gcc/ada/s-bignum.adb
diff options
context:
space:
mode:
Diffstat (limited to 'gcc-4.9/gcc/ada/s-bignum.adb')
-rw-r--r--gcc-4.9/gcc/ada/s-bignum.adb1105
1 files changed, 1105 insertions, 0 deletions
diff --git a/gcc-4.9/gcc/ada/s-bignum.adb b/gcc-4.9/gcc/ada/s-bignum.adb
new file mode 100644
index 000000000..d12a9e3cd
--- /dev/null
+++ b/gcc-4.9/gcc/ada/s-bignum.adb
@@ -0,0 +1,1105 @@
+------------------------------------------------------------------------------
+-- --
+-- GNAT COMPILER COMPONENTS --
+-- --
+-- S Y S T E M . B I G N U M S --
+-- --
+-- B o d y --
+-- --
+-- Copyright (C) 2012-2013, 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 for
+-- use in computing intermediate values in expressions for the case where
+-- pragma Overflow_Check (Eliminate) is in effect.
+
+with System; use System;
+with System.Secondary_Stack; use System.Secondary_Stack;
+with System.Storage_Elements; use System.Storage_Elements;
+
+package body System.Bignums is
+
+ use Interfaces;
+ -- So that operations on Unsigned_32 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
+ -- Change the if False here to if True to get allocation on the heap
+ -- instead of the secondary stack, which is convenient for debugging
+ -- System.Bignum itself.
+
+ if False then
+ declare
+ B : Bignum;
+ begin
+ B := new Bignum_Data'(Len, False, (others => 0));
+ return B;
+ end;
+
+ -- Normal case of 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 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;
+
+end System.Bignums;