CbC/CbC_gcc: gcc/ada/libgnat/s-genbig.adb comparison

comparison 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

comparison

equal deleted inserted replaced

-:84e7813d76e9
+:1830386684a0
+------------------------------------------------------------------------------
+--                                                                          --
+--                         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;

Mercurial > hg > CbC > CbC_gcc

comparison gcc/ada/libgnat/s-genbig.adb @ 145:1830386684a0