gcc/ada/s-bignum.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT COMPILER COMPONENTS                         --
   4 --                                                                          --
   5 --                       S Y S T E M . B I G N U M S                        --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --          Copyright (C) 2012-2015, Free Software Foundation, Inc.         --
  10 --                                                                          --
  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.                                     --
  17 --                                                                          --
  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.               --
  21 --                                                                          --
  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/>.                                          --
  26 --                                                                          --
  27 -- GNAT was originally developed  by the GNAT team at  New York University. --
  28 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  29 --                                                                          --
  30 ------------------------------------------------------------------------------
  31
  32 --  This package provides arbitrary precision signed integer arithmetic for
  33 --  use in computing intermediate values in expressions for the case where
  34 --  pragma Overflow_Check (Eliminate) is in effect.
  35
  36 with System;                  use System;
  37 with System.Secondary_Stack;  use System.Secondary_Stack;
  38 with System.Storage_Elements; use System.Storage_Elements;
  39
  40 package body System.Bignums is
  41
  42    use Interfaces;
  43    --  So that operations on Unsigned_32 are available
  44
  45    type DD is mod Base ** 2;
  46    --  Double length digit used for intermediate computations
  47
  48    function MSD (X : DD) return SD is (SD (X / Base));
  49    function LSD (X : DD) return SD is (SD (X mod Base));
  50    --  Most significant and least significant digit of double digit value
  51
  52    function "&" (X, Y : SD) return DD is (DD (X) * Base + DD (Y));
  53    --  Compose double digit value from two single digit values
  54
  55    subtype LLI is Long_Long_Integer;
  56
  57    One_Data : constant Digit_Vector (1 .. 1) := (1 => 1);
  58    --  Constant one
  59
  60    Zero_Data : constant Digit_Vector (1 .. 0) := (1 .. 0 => 0);
  61    --  Constant zero
  62
  63    -----------------------
  64    -- Local Subprograms --
  65    -----------------------
  66
  67    function Add
  68      (X, Y  : Digit_Vector;
  69       X_Neg : Boolean;
  70       Y_Neg : Boolean) return Bignum
  71    with
  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.
  76
  77    function Allocate_Bignum (Len : Length) return Bignum with
  78      Post => Allocate_Bignum'Result.Len = Len;
  79    --  Allocate Bignum value of indicated length on secondary stack. On return
  80    --  the Neg and D fields are left uninitialized.
  81
  82    type Compare_Result is (LT, EQ, GT);
  83    --  Indicates result of comparison in following call
  84
  85    function Compare
  86      (X, Y         : Digit_Vector;
  87       X_Neg, Y_Neg : Boolean) return Compare_Result
  88    with
  89      Pre => X'First = 1 and then Y'First = 1;
  90    --  Compare (X with sign X_Neg) with (Y with sign Y_Neg), and return the
  91    --  result of the signed comparison.
  92
  93    procedure Div_Rem
  94      (X, Y              : Bignum;
  95       Quotient          : out Bignum;
  96       Remainder         : out Bignum;
  97       Discard_Quotient  : Boolean := False;
  98       Discard_Remainder : Boolean := False);
  99    --  Returns the Quotient and Remainder from dividing abs (X) by abs (Y). The
 100    --  values of X and Y are not modified. If Discard_Quotient is True, then
 101    --  Quotient is undefined on return, and if Discard_Remainder is True, then
 102    --  Remainder is undefined on return. Service routine for Big_Div/Rem/Mod.
 103
 104    procedure Free_Bignum (X : Bignum) is null;
 105    --  Called to free a Bignum value used in intermediate computations. In
 106    --  this implementation using the secondary stack, it does nothing at all,
 107    --  because we rely on Mark/Release, but it may be of use for some
 108    --  alternative implementation.
 109
 110    function Normalize
 111      (X   : Digit_Vector;
 112       Neg : Boolean := False) return Bignum;
 113    --  Given a digit vector and sign, allocate and construct a Bignum value.
 114    --  Note that X may have leading zeroes which must be removed, and if the
 115    --  result is zero, the sign is forced positive.
 116
 117    ---------
 118    -- Add --
 119    ---------
 120
 121    function Add
 122      (X, Y  : Digit_Vector;
 123       X_Neg : Boolean;
 124       Y_Neg : Boolean) return Bignum
 125    is
 126    begin
 127       --  If signs are the same, we are doing an addition, it is convenient to
 128       --  ensure that the first operand is the longer of the two.
 129
 130       if X_Neg = Y_Neg then
 131          if X'Last < Y'Last then
 132             return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg);
 133
 134          --  Here signs are the same, and the first operand is the longer
 135
 136          else
 137             pragma Assert (X_Neg = Y_Neg and then X'Last >= Y'Last);
 138
 139             --  Do addition, putting result in Sum (allowing for carry)
 140
 141             declare
 142                Sum : Digit_Vector (0 .. X'Last);
 143                RD  : DD;
 144
 145             begin
 146                RD := 0;
 147                for J in reverse 1 .. X'Last loop
 148                   RD := RD + DD (X (J));
 149
 150                   if J >= 1 + (X'Last - Y'Last) then
 151                      RD := RD + DD (Y (J - (X'Last - Y'Last)));
 152                   end if;
 153
 154                   Sum (J) := LSD (RD);
 155                   RD := RD / Base;
 156                end loop;
 157
 158                Sum (0) := SD (RD);
 159                return Normalize (Sum, X_Neg);
 160             end;
 161          end if;
 162
 163       --  Signs are different so really this is a subtraction, we want to make
 164       --  sure that the largest magnitude operand is the first one, and then
 165       --  the result will have the sign of the first operand.
 166
 167       else
 168          declare
 169             CR : constant Compare_Result := Compare (X, Y, False, False);
 170
 171          begin
 172             if CR = EQ then
 173                return Normalize (Zero_Data);
 174
 175             elsif CR = LT then
 176                return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg);
 177
 178             else
 179                pragma Assert (X_Neg /= Y_Neg and then CR = GT);
 180
 181                --  Do subtraction, putting result in Diff
 182
 183                declare
 184                   Diff : Digit_Vector (1 .. X'Length);
 185                   RD   : DD;
 186
 187                begin
 188                   RD := 0;
 189                   for J in reverse 1 .. X'Last loop
 190                      RD := RD + DD (X (J));
 191
 192                      if J >= 1 + (X'Last - Y'Last) then
 193                         RD := RD - DD (Y (J - (X'Last - Y'Last)));
 194                      end if;
 195
 196                      Diff (J) := LSD (RD);
 197                      RD := (if RD < Base then 0 else -1);
 198                   end loop;
 199
 200                   return Normalize (Diff, X_Neg);
 201                end;
 202             end if;
 203          end;
 204       end if;
 205    end Add;
 206
 207    ---------------------
 208    -- Allocate_Bignum --
 209    ---------------------
 210
 211    function Allocate_Bignum (Len : Length) return Bignum is
 212       Addr : Address;
 213
 214    begin
 215       --  Change the if False here to if True to get allocation on the heap
 216       --  instead of the secondary stack, which is convenient for debugging
 217       --  System.Bignum itself.
 218
 219       if False then
 220          declare
 221             B : Bignum;
 222          begin
 223             B := new Bignum_Data'(Len, False, (others => 0));
 224             return B;
 225          end;
 226
 227       --  Normal case of allocation on the secondary stack
 228
 229       else
 230          --  Note: The approach used here is designed to avoid strict aliasing
 231          --  warnings that appeared previously using unchecked conversion.
 232
 233          SS_Allocate (Addr, Storage_Offset (4 + 4 * Len));
 234
 235          declare
 236             B : Bignum;
 237             for B'Address use Addr'Address;
 238             pragma Import (Ada, B);
 239
 240             BD : Bignum_Data (Len);
 241             for BD'Address use Addr;
 242             pragma Import (Ada, BD);
 243
 244             --  Expose a writable view of discriminant BD.Len so that we can
 245             --  initialize it. We need to use the exact layout of the record
 246             --  to ensure that the Length field has 24 bits as expected.
 247
 248             type Bignum_Data_Header is record
 249                Len : Length;
 250                Neg : Boolean;
 251             end record;
 252
 253             for Bignum_Data_Header use record
 254                Len at 0 range 0 .. 23;
 255                Neg at 3 range 0 .. 7;
 256             end record;
 257
 258             BDH : Bignum_Data_Header;
 259             for BDH'Address use BD'Address;
 260             pragma Import (Ada, BDH);
 261
 262             pragma Assert (BDH.Len'Size = BD.Len'Size);
 263
 264          begin
 265             BDH.Len := Len;
 266             return B;
 267          end;
 268       end if;
 269    end Allocate_Bignum;
 270
 271    -------------
 272    -- Big_Abs --
 273    -------------
 274
 275    function Big_Abs (X : Bignum) return Bignum is
 276    begin
 277       return Normalize (X.D);
 278    end Big_Abs;
 279
 280    -------------
 281    -- Big_Add --
 282    -------------
 283
 284    function Big_Add  (X, Y : Bignum) return Bignum is
 285    begin
 286       return Add (X.D, Y.D, X.Neg, Y.Neg);
 287    end Big_Add;
 288
 289    -------------
 290    -- Big_Div --
 291    -------------
 292
 293    --  This table is excerpted from RM 4.5.5(28-30) and shows how the result
 294    --  varies with the signs of the operands.
 295
 296    --   A      B   A/B      A     B    A/B
 297    --
 298    --   10     5    2      -10    5    -2
 299    --   11     5    2      -11    5    -2
 300    --   12     5    2      -12    5    -2
 301    --   13     5    2      -13    5    -2
 302    --   14     5    2      -14    5    -2
 303    --
 304    --   A      B   A/B      A     B    A/B
 305    --
 306    --   10    -5   -2      -10   -5     2
 307    --   11    -5   -2      -11   -5     2
 308    --   12    -5   -2      -12   -5     2
 309    --   13    -5   -2      -13   -5     2
 310    --   14    -5   -2      -14   -5     2
 311
 312    function Big_Div  (X, Y : Bignum) return Bignum is
 313       Q, R : Bignum;
 314    begin
 315       Div_Rem (X, Y, Q, R, Discard_Remainder => True);
 316       Q.Neg := Q.Len > 0 and then (X.Neg xor Y.Neg);
 317       return Q;
 318    end Big_Div;
 319
 320    -------------
 321    -- Big_Exp --
 322    -------------
 323
 324    function Big_Exp  (X, Y : Bignum) return Bignum is
 325
 326       function "**" (X : Bignum; Y : SD) return Bignum;
 327       --  Internal routine where we know right operand is one word
 328
 329       ----------
 330       -- "**" --
 331       ----------
 332
 333       function "**" (X : Bignum; Y : SD) return Bignum is
 334       begin
 335          case Y is
 336
 337             --  X ** 0 is 1
 338
 339             when 0 =>
 340                return Normalize (One_Data);
 341
 342             --  X ** 1 is X
 343
 344             when 1 =>
 345                return Normalize (X.D);
 346
 347             --  X ** 2 is X * X
 348
 349             when 2 =>
 350                return Big_Mul (X, X);
 351
 352             --  For X greater than 2, use the recursion
 353
 354             --  X even, X ** Y = (X ** (Y/2)) ** 2;
 355             --  X odd,  X ** Y = (X ** (Y/2)) ** 2 * X;
 356
 357             when others =>
 358                declare
 359                   XY2  : constant Bignum := X ** (Y / 2);
 360                   XY2S : constant Bignum := Big_Mul (XY2, XY2);
 361                   Res  : Bignum;
 362
 363                begin
 364                   Free_Bignum (XY2);
 365
 366                   --  Raise storage error if intermediate value is getting too
 367                   --  large, which we arbitrarily define as 200 words for now.
 368
 369                   if XY2S.Len > 200 then
 370                      Free_Bignum (XY2S);
 371                      raise Storage_Error with
 372                        "exponentiation result is too large";
 373                   end if;
 374
 375                   --  Otherwise take care of even/odd cases
 376
 377                   if (Y and 1) = 0 then
 378                      return XY2S;
 379
 380                   else
 381                      Res := Big_Mul (XY2S, X);
 382                      Free_Bignum (XY2S);
 383                      return Res;
 384                   end if;
 385                end;
 386          end case;
 387       end "**";
 388
 389    --  Start of processing for Big_Exp
 390
 391    begin
 392       --  Error if right operand negative
 393
 394       if Y.Neg then
 395          raise Constraint_Error with "exponentiation to negative power";
 396
 397       --  X ** 0 is always 1 (including 0 ** 0, so do this test first)
 398
 399       elsif Y.Len = 0 then
 400          return Normalize (One_Data);
 401
 402       --  0 ** X is always 0 (for X non-zero)
 403
 404       elsif X.Len = 0 then
 405          return Normalize (Zero_Data);
 406
 407       --  (+1) ** Y = 1
 408       --  (-1) ** Y = +/-1 depending on whether Y is even or odd
 409
 410       elsif X.Len = 1 and then X.D (1) = 1 then
 411          return Normalize
 412            (X.D, Neg => X.Neg and then ((Y.D (Y.Len) and 1) = 1));
 413
 414       --  If the absolute value of the base is greater than 1, then the
 415       --  exponent must not be bigger than one word, otherwise the result
 416       --  is ludicrously large, and we just signal Storage_Error right away.
 417
 418       elsif Y.Len > 1 then
 419          raise Storage_Error with "exponentiation result is too large";
 420
 421       --  Special case (+/-)2 ** K, where K is 1 .. 31 using a shift
 422
 423       elsif X.Len = 1 and then X.D (1) = 2 and then Y.D (1) < 32 then
 424          declare
 425             D : constant Digit_Vector (1 .. 1) :=
 426                   (1 => Shift_Left (SD'(1), Natural (Y.D (1))));
 427          begin
 428             return Normalize (D, X.Neg);
 429          end;
 430
 431       --  Remaining cases have right operand of one word
 432
 433       else
 434          return X ** Y.D (1);
 435       end if;
 436    end Big_Exp;
 437
 438    ------------
 439    -- Big_EQ --
 440    ------------
 441
 442    function Big_EQ (X, Y : Bignum) return Boolean is
 443    begin
 444       return Compare (X.D, Y.D, X.Neg, Y.Neg) = EQ;
 445    end Big_EQ;
 446
 447    ------------
 448    -- Big_GE --
 449    ------------
 450
 451    function Big_GE (X, Y : Bignum) return Boolean is
 452    begin
 453       return Compare (X.D, Y.D, X.Neg, Y.Neg) /= LT;
 454    end Big_GE;
 455
 456    ------------
 457    -- Big_GT --
 458    ------------
 459
 460    function Big_GT (X, Y : Bignum) return Boolean is
 461    begin
 462       return Compare (X.D, Y.D, X.Neg, Y.Neg) = GT;
 463    end Big_GT;
 464
 465    ------------
 466    -- Big_LE --
 467    ------------
 468
 469    function Big_LE (X, Y : Bignum) return Boolean is
 470    begin
 471       return Compare (X.D, Y.D, X.Neg, Y.Neg) /= GT;
 472    end Big_LE;
 473
 474    ------------
 475    -- Big_LT --
 476    ------------
 477
 478    function Big_LT (X, Y : Bignum) return Boolean is
 479    begin
 480       return Compare (X.D, Y.D, X.Neg, Y.Neg) = LT;
 481    end Big_LT;
 482
 483    -------------
 484    -- Big_Mod --
 485    -------------
 486
 487    --  This table is excerpted from RM 4.5.5(28-30) and shows how the result
 488    --  of Rem and Mod vary with the signs of the operands.
 489
 490    --   A      B    A mod B  A rem B     A     B    A mod B  A rem B
 491
 492    --   10     5       0        0       -10    5       0        0
 493    --   11     5       1        1       -11    5       4       -1
 494    --   12     5       2        2       -12    5       3       -2
 495    --   13     5       3        3       -13    5       2       -3
 496    --   14     5       4        4       -14    5       1       -4
 497
 498    --   A      B    A mod B  A rem B     A     B    A mod B  A rem B
 499
 500    --   10    -5       0        0       -10   -5       0        0
 501    --   11    -5      -4        1       -11   -5      -1       -1
 502    --   12    -5      -3        2       -12   -5      -2       -2
 503    --   13    -5      -2        3       -13   -5      -3       -3
 504    --   14    -5      -1        4       -14   -5      -4       -4
 505
 506    function Big_Mod (X, Y : Bignum) return Bignum is
 507       Q, R : Bignum;
 508
 509    begin
 510       --  If signs are same, result is same as Rem
 511
 512       if X.Neg = Y.Neg then
 513          return Big_Rem (X, Y);
 514
 515       --  Case where Mod is different
 516
 517       else
 518          --  Do division
 519
 520          Div_Rem (X, Y, Q, R, Discard_Quotient => True);
 521
 522          --  Zero result is unchanged
 523
 524          if R.Len = 0 then
 525             return R;
 526
 527          --  Otherwise adjust result
 528
 529          else
 530             declare
 531                T1 : constant Bignum := Big_Sub (Y, R);
 532             begin
 533                T1.Neg := Y.Neg;
 534                Free_Bignum (R);
 535                return T1;
 536             end;
 537          end if;
 538       end if;
 539    end Big_Mod;
 540
 541    -------------
 542    -- Big_Mul --
 543    -------------
 544
 545    function Big_Mul (X, Y : Bignum) return Bignum is
 546       Result : Digit_Vector (1 .. X.Len + Y.Len) := (others => 0);
 547       --  Accumulate result (max length of result is sum of operand lengths)
 548
 549       L : Length;
 550       --  Current result digit
 551
 552       D : DD;
 553       --  Result digit
 554
 555    begin
 556       for J in 1 .. X.Len loop
 557          for K in 1 .. Y.Len loop
 558             L := Result'Last - (X.Len - J) - (Y.Len - K);
 559             D := DD (X.D (J)) * DD (Y.D (K)) + DD (Result (L));
 560             Result (L) := LSD (D);
 561             D := D / Base;
 562
 563             --  D is carry which must be propagated
 564
 565             while D /= 0 and then L >= 1 loop
 566                L := L - 1;
 567                D := D + DD (Result (L));
 568                Result (L) := LSD (D);
 569                D := D / Base;
 570             end loop;
 571
 572             --  Must not have a carry trying to extend max length
 573
 574             pragma Assert (D = 0);
 575          end loop;
 576       end loop;
 577
 578       --  Return result
 579
 580       return Normalize (Result, X.Neg xor Y.Neg);
 581    end Big_Mul;
 582
 583    ------------
 584    -- Big_NE --
 585    ------------
 586
 587    function Big_NE (X, Y : Bignum) return Boolean is
 588    begin
 589       return Compare (X.D, Y.D, X.Neg, Y.Neg) /= EQ;
 590    end Big_NE;
 591
 592    -------------
 593    -- Big_Neg --
 594    -------------
 595
 596    function Big_Neg (X : Bignum) return Bignum is
 597    begin
 598       return Normalize (X.D, not X.Neg);
 599    end Big_Neg;
 600
 601    -------------
 602    -- Big_Rem --
 603    -------------
 604
 605    --  This table is excerpted from RM 4.5.5(28-30) and shows how the result
 606    --  varies with the signs of the operands.
 607
 608    --   A      B   A rem B   A     B   A rem B
 609
 610    --   10     5      0     -10    5      0
 611    --   11     5      1     -11    5     -1
 612    --   12     5      2     -12    5     -2
 613    --   13     5      3     -13    5     -3
 614    --   14     5      4     -14    5     -4
 615
 616    --   A      B  A rem B    A     B   A rem B
 617
 618    --   10    -5     0      -10   -5      0
 619    --   11    -5     1      -11   -5     -1
 620    --   12    -5     2      -12   -5     -2
 621    --   13    -5     3      -13   -5     -3
 622    --   14    -5     4      -14   -5     -4
 623
 624    function Big_Rem (X, Y : Bignum) return Bignum is
 625       Q, R : Bignum;
 626    begin
 627       Div_Rem (X, Y, Q, R, Discard_Quotient => True);
 628       R.Neg := R.Len > 0 and then X.Neg;
 629       return R;
 630    end Big_Rem;
 631
 632    -------------
 633    -- Big_Sub --
 634    -------------
 635
 636    function Big_Sub (X, Y : Bignum) return Bignum is
 637    begin
 638       --  If right operand zero, return left operand (avoiding sharing)
 639
 640       if Y.Len = 0 then
 641          return Normalize (X.D, X.Neg);
 642
 643       --  Otherwise add negative of right operand
 644
 645       else
 646          return Add (X.D, Y.D, X.Neg, not Y.Neg);
 647       end if;
 648    end Big_Sub;
 649
 650    -------------
 651    -- Compare --
 652    -------------
 653
 654    function Compare
 655      (X, Y         : Digit_Vector;
 656       X_Neg, Y_Neg : Boolean) return Compare_Result
 657    is
 658    begin
 659       --  Signs are different, that's decisive, since 0 is always plus
 660
 661       if X_Neg /= Y_Neg then
 662          return (if X_Neg then LT else GT);
 663
 664       --  Lengths are different, that's decisive since no leading zeroes
 665
 666       elsif X'Last /= Y'Last then
 667          return (if (X'Last > Y'Last) xor X_Neg then GT else LT);
 668
 669       --  Need to compare data
 670
 671       else
 672          for J in X'Range loop
 673             if X (J) /= Y (J) then
 674                return (if (X (J) > Y (J)) xor X_Neg then GT else LT);
 675             end if;
 676          end loop;
 677
 678          return EQ;
 679       end if;
 680    end Compare;
 681
 682    -------------
 683    -- Div_Rem --
 684    -------------
 685
 686    procedure Div_Rem
 687      (X, Y              : Bignum;
 688       Quotient          : out Bignum;
 689       Remainder         : out Bignum;
 690       Discard_Quotient  : Boolean := False;
 691       Discard_Remainder : Boolean := False)
 692    is
 693    begin
 694       --  Error if division by zero
 695
 696       if Y.Len = 0 then
 697          raise Constraint_Error with "division by zero";
 698       end if;
 699
 700       --  Handle simple cases with special tests
 701
 702       --  If X < Y then quotient is zero and remainder is X
 703
 704       if Compare (X.D, Y.D, False, False) = LT then
 705          Remainder := Normalize (X.D);
 706          Quotient  := Normalize (Zero_Data);
 707          return;
 708
 709       --  If both X and Y are less than 2**63-1, we can use Long_Long_Integer
 710       --  arithmetic. Note it is good not to do an accurate range check against
 711       --  Long_Long_Integer since -2**63 / -1 overflows.
 712
 713       elsif (X.Len <= 1 or else (X.Len = 2 and then X.D (1) < 2**31))
 714               and then
 715             (Y.Len <= 1 or else (Y.Len = 2 and then Y.D (1) < 2**31))
 716       then
 717          declare
 718             A : constant LLI := abs (From_Bignum (X));
 719             B : constant LLI := abs (From_Bignum (Y));
 720          begin
 721             Quotient  := To_Bignum (A / B);
 722             Remainder := To_Bignum (A rem B);
 723             return;
 724          end;
 725
 726       --  Easy case if divisor is one digit
 727
 728       elsif Y.Len = 1 then
 729          declare
 730             ND  : DD;
 731             Div : constant DD := DD (Y.D (1));
 732
 733             Result : Digit_Vector (1 .. X.Len);
 734             Remdr  : Digit_Vector (1 .. 1);
 735
 736          begin
 737             ND := 0;
 738             for J in 1 .. X.Len loop
 739                ND := Base * ND + DD (X.D (J));
 740                Result (J) := SD (ND / Div);
 741                ND := ND rem Div;
 742             end loop;
 743
 744             Quotient  := Normalize (Result);
 745             Remdr (1) := SD (ND);
 746             Remainder := Normalize (Remdr);
 747             return;
 748          end;
 749       end if;
 750
 751       --  The complex full multi-precision case. We will employ algorithm
 752       --  D defined in the section "The Classical Algorithms" (sec. 4.3.1)
 753       --  of Donald Knuth's "The Art of Computer Programming", Vol. 2, 2nd
 754       --  edition. The terminology is adjusted for this section to match that
 755       --  reference.
 756
 757       --  We are dividing X.Len digits of X (called u here) by Y.Len digits
 758       --  of Y (called v here), developing the quotient and remainder. The
 759       --  numbers are represented using Base, which was chosen so that we have
 760       --  the operations of multiplying to single digits (SD) to form a double
 761       --  digit (DD), and dividing a double digit (DD) by a single digit (SD)
 762       --  to give a single digit quotient and a single digit remainder.
 763
 764       --  Algorithm D from Knuth
 765
 766       --  Comments here with square brackets are directly from Knuth
 767
 768       Algorithm_D : declare
 769
 770          --  The following lower case variables correspond exactly to the
 771          --  terminology used in algorithm D.
 772
 773          m : constant Length := X.Len - Y.Len;
 774          n : constant Length := Y.Len;
 775          b : constant DD     := Base;
 776
 777          u : Digit_Vector (0 .. m + n);
 778          v : Digit_Vector (1 .. n);
 779          q : Digit_Vector (0 .. m);
 780          r : Digit_Vector (1 .. n);
 781
 782          u0 : SD renames u (0);
 783          v1 : SD renames v (1);
 784          v2 : SD renames v (2);
 785
 786          d    : DD;
 787          j    : Length;
 788          qhat : DD;
 789          rhat : DD;
 790          temp : DD;
 791
 792       begin
 793          --  Initialize data of left and right operands
 794
 795          for J in 1 .. m + n loop
 796             u (J) := X.D (J);
 797          end loop;
 798
 799          for J in 1 .. n loop
 800             v (J) := Y.D (J);
 801          end loop;
 802
 803          --  [Division of nonnegative integers.] Given nonnegative integers u
 804          --  = (ul,u2..um+n) and v = (v1,v2..vn), where v1 /= 0 and n > 1, we
 805          --  form the quotient u / v = (q0,ql..qm) and the remainder u mod v =
 806          --  (r1,r2..rn).
 807
 808          pragma Assert (v1 /= 0);
 809          pragma Assert (n > 1);
 810
 811          --  Dl. [Normalize.] Set d = b/(vl + 1). Then set (u0,u1,u2..um+n)
 812          --  equal to (u1,u2..um+n) times d, and set (v1,v2..vn) equal to
 813          --  (v1,v2..vn) times d. Note the introduction of a new digit position
 814          --  u0 at the left of u1; if d = 1 all we need to do in this step is
 815          --  to set u0 = 0.
 816
 817          d := b / (DD (v1) + 1);
 818
 819          if d = 1 then
 820             u0 := 0;
 821
 822          else
 823             declare
 824                Carry : DD;
 825                Tmp   : DD;
 826
 827             begin
 828                --  Multiply Dividend (u) by d
 829
 830                Carry := 0;
 831                for J in reverse 1 .. m + n loop
 832                   Tmp   := DD (u (J)) * d + Carry;
 833                   u (J) := LSD (Tmp);
 834                   Carry := Tmp / Base;
 835                end loop;
 836
 837                u0 := SD (Carry);
 838
 839                --  Multiply Divisor (v) by d
 840
 841                Carry := 0;
 842                for J in reverse 1 .. n loop
 843                   Tmp   := DD (v (J)) * d + Carry;
 844                   v (J) := LSD (Tmp);
 845                   Carry := Tmp / Base;
 846                end loop;
 847
 848                pragma Assert (Carry = 0);
 849             end;
 850          end if;
 851
 852          --  D2. [Initialize j.] Set j = 0. The loop on j, steps D2 through D7,
 853          --  will be essentially a division of (uj, uj+1..uj+n) by (v1,v2..vn)
 854          --  to get a single quotient digit qj.
 855
 856          j := 0;
 857
 858          --  Loop through digits
 859
 860          loop
 861             --  Note: In the original printing, step D3 was as follows:
 862
 863             --  D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise
 864             --  set qhat to (uj,uj+1)/v1. Now test if v2 * qhat is greater than
 865             --  (uj*b + uj+1 - qhat*v1)*b + uj+2. If so, decrease qhat by 1 and
 866             --  repeat this test
 867
 868             --  This had a bug not discovered till 1995, see Vol 2 errata:
 869             --  http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz. Under
 870             --  rare circumstances the expression in the test could overflow.
 871             --  This version was further corrected in 2005, see Vol 2 errata:
 872             --  http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz.
 873             --  The code below is the fixed version of this step.
 874
 875             --  D3. [Calculate qhat.] Set qhat to (uj,uj+1)/v1 and rhat to
 876             --  to (uj,uj+1) mod v1.
 877
 878             temp := u (j) & u (j + 1);
 879             qhat := temp / DD (v1);
 880             rhat := temp mod DD (v1);
 881
 882             --  D3 (continued). Now test if qhat >= b or v2*qhat > (rhat,uj+2):
 883             --  if so, decrease qhat by 1, increase rhat by v1, and repeat this
 884             --  test if rhat < b. [The test on v2 determines at high speed
 885             --  most of the cases in which the trial value qhat is one too
 886             --  large, and eliminates all cases where qhat is two too large.]
 887
 888             while qhat >= b
 889               or else DD (v2) * qhat > LSD (rhat) & u (j + 2)
 890             loop
 891                qhat := qhat - 1;
 892                rhat := rhat + DD (v1);
 893                exit when rhat >= b;
 894             end loop;
 895
 896             --  D4. [Multiply and subtract.] Replace (uj,uj+1..uj+n) by
 897             --  (uj,uj+1..uj+n) minus qhat times (v1,v2..vn). This step
 898             --  consists of a simple multiplication by a one-place number,
 899             --  combined with a subtraction.
 900
 901             --  The digits (uj,uj+1..uj+n) are always kept positive; if the
 902             --  result of this step is actually negative then (uj,uj+1..uj+n)
 903             --  is left as the true value plus b**(n+1), i.e. as the b's
 904             --  complement of the true value, and a "borrow" to the left is
 905             --  remembered.
 906
 907             declare
 908                Borrow : SD;
 909                Carry  : DD;
 910                Temp   : DD;
 911
 912                Negative : Boolean;
 913                --  Records if subtraction causes a negative result, requiring
 914                --  an add back (case where qhat turned out to be 1 too large).
 915
 916             begin
 917                Borrow := 0;
 918                for K in reverse 1 .. n loop
 919                   Temp := qhat * DD (v (K)) + DD (Borrow);
 920                   Borrow := MSD (Temp);
 921
 922                   if LSD (Temp) > u (j + K) then
 923                      Borrow := Borrow + 1;
 924                   end if;
 925
 926                   u (j + K) := u (j + K) - LSD (Temp);
 927                end loop;
 928
 929                Negative := u (j) < Borrow;
 930                u (j) := u (j) - Borrow;
 931
 932                --  D5. [Test remainder.] Set qj = qhat. If the result of step
 933                --  D4 was negative, we will do the add back step (step D6).
 934
 935                q (j) := LSD (qhat);
 936
 937                if Negative then
 938
 939                   --  D6. [Add back.] Decrease qj by 1, and add (0,v1,v2..vn)
 940                   --  to (uj,uj+1,uj+2..uj+n). (A carry will occur to the left
 941                   --  of uj, and it is be ignored since it cancels with the
 942                   --  borrow that occurred in D4.)
 943
 944                   q (j) := q (j) - 1;
 945
 946                   Carry := 0;
 947                   for K in reverse 1 .. n loop
 948                      Temp := DD (v (K)) + DD (u (j + K)) + Carry;
 949                      u (j + K) := LSD (Temp);
 950                      Carry := Temp / Base;
 951                   end loop;
 952
 953                   u (j) := u (j) + SD (Carry);
 954                end if;
 955             end;
 956
 957             --  D7. [Loop on j.] Increase j by one. Now if j <= m, go back to
 958             --  D3 (the start of the loop on j).
 959
 960             j := j + 1;
 961             exit when not (j <= m);
 962          end loop;
 963
 964          --  D8. [Unnormalize.] Now (qo,ql..qm) is the desired quotient, and
 965          --  the desired remainder may be obtained by dividing (um+1..um+n)
 966          --  by d.
 967
 968          if not Discard_Quotient then
 969             Quotient := Normalize (q);
 970          end if;
 971
 972          if not Discard_Remainder then
 973             declare
 974                Remdr : DD;
 975
 976             begin
 977                Remdr := 0;
 978                for K in 1 .. n loop
 979                   Remdr := Base * Remdr + DD (u (m + K));
 980                   r (K) := SD (Remdr / d);
 981                   Remdr := Remdr rem d;
 982                end loop;
 983
 984                pragma Assert (Remdr = 0);
 985             end;
 986
 987             Remainder := Normalize (r);
 988          end if;
 989       end Algorithm_D;
 990    end Div_Rem;
 991
 992    -----------------
 993    -- From_Bignum --
 994    -----------------
 995
 996    function From_Bignum (X : Bignum) return Long_Long_Integer is
 997    begin
 998       if X.Len = 0 then
 999          return 0;
1000
1001       elsif X.Len = 1 then
1002          return (if X.Neg then -LLI (X.D (1)) else LLI (X.D (1)));
1003
1004       elsif X.Len = 2 then
1005          declare
1006             Mag : constant DD := X.D (1) & X.D (2);
1007          begin
1008             if X.Neg and then Mag <= 2 ** 63 then
1009                return -LLI (Mag);
1010             elsif Mag < 2 ** 63 then
1011                return LLI (Mag);
1012             end if;
1013          end;
1014       end if;
1015
1016       raise Constraint_Error with "expression value out of range";
1017    end From_Bignum;
1018
1019    -------------------------
1020    -- Bignum_In_LLI_Range --
1021    -------------------------
1022
1023    function Bignum_In_LLI_Range (X : Bignum) return Boolean is
1024    begin
1025       --  If length is 0 or 1, definitely fits
1026
1027       if X.Len <= 1 then
1028          return True;
1029
1030       --  If length is greater than 2, definitely does not fit
1031
1032       elsif X.Len > 2 then
1033          return False;
1034
1035       --  Length is 2, more tests needed
1036
1037       else
1038          declare
1039             Mag : constant DD := X.D (1) & X.D (2);
1040          begin
1041             return Mag < 2 ** 63 or else (X.Neg and then Mag = 2 ** 63);
1042          end;
1043       end if;
1044    end Bignum_In_LLI_Range;
1045
1046    ---------------
1047    -- Normalize --
1048    ---------------
1049
1050    function Normalize
1051      (X   : Digit_Vector;
1052       Neg : Boolean := False) return Bignum
1053    is
1054       B : Bignum;
1055       J : Length;
1056
1057    begin
1058       J := X'First;
1059       while J <= X'Last and then X (J) = 0 loop
1060          J := J + 1;
1061       end loop;
1062
1063       B := Allocate_Bignum (X'Last - J + 1);
1064       B.Neg := B.Len > 0 and then Neg;
1065       B.D := X (J .. X'Last);
1066       return B;
1067    end Normalize;
1068
1069    ---------------
1070    -- To_Bignum --
1071    ---------------
1072
1073    function To_Bignum (X : Long_Long_Integer) return Bignum is
1074       R : Bignum;
1075
1076    begin
1077       if X = 0 then
1078          R := Allocate_Bignum (0);
1079
1080       --  One word result
1081
1082       elsif X in -(2 ** 32 - 1) .. +(2 ** 32 - 1) then
1083          R := Allocate_Bignum (1);
1084          R.D (1) := SD (abs (X));
1085
1086       --  Largest negative number annoyance
1087
1088       elsif X = Long_Long_Integer'First then
1089          R := Allocate_Bignum (2);
1090          R.D (1) := 2 ** 31;
1091          R.D (2) := 0;
1092
1093       --  Normal two word case
1094
1095       else
1096          R := Allocate_Bignum (2);
1097          R.D (2) := SD (abs (X) mod Base);
1098          R.D (1) := SD (abs (X) / Base);
1099       end if;
1100
1101       R.Neg := X < 0;
1102       return R;
1103    end To_Bignum;
1104
1105 end System.Bignums;