gcc/ada/libgnat/s-genbig.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT COMPILER COMPONENTS                         --
   4 --                                                                          --
   5 --               S Y S T E M . G E N E R I C _ B I G N U M S                --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --          Copyright (C) 2012-2024, 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.
  33
  34 package body System.Generic_Bignums is
  35
  36    use Interfaces;
  37    --  So that operations on Unsigned_32/Unsigned_64 are available
  38
  39    use Shared_Bignums;
  40
  41    type DD is mod Base ** 2;
  42    --  Double length digit used for intermediate computations
  43
  44    function MSD (X : DD) return SD is (SD (X / Base));
  45    function LSD (X : DD) return SD is (SD (X mod Base));
  46    --  Most significant and least significant digit of double digit value
  47
  48    function "&" (X, Y : SD) return DD is (DD (X) * Base + DD (Y));
  49    --  Compose double digit value from two single digit values
  50
  51    subtype LLI is Long_Long_Integer;
  52    subtype LLLI is Long_Long_Long_Integer;
  53
  54    LLLI_Is_128 : constant Boolean := Long_Long_Long_Integer'Size = 128;
  55    --  True if Long_Long_Long_Integer is 128-bit large
  56
  57    One_Data : constant Digit_Vector (1 .. 1) := [1];
  58    --  Constant one
  59
  60    Zero_Data : constant Digit_Vector (1 .. 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 Big_Integer
  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    type Compare_Result is (LT, EQ, GT);
  78    --  Indicates result of comparison in following call
  79
  80    function Compare
  81      (X, Y         : Digit_Vector;
  82       X_Neg, Y_Neg : Boolean) return Compare_Result
  83    with
  84      Pre => X'First = 1 and then Y'First = 1;
  85    --  Compare (X with sign X_Neg) with (Y with sign Y_Neg), and return the
  86    --  result of the signed comparison.
  87
  88    procedure Div_Rem
  89      (X, Y              : Bignum;
  90       Quotient          : out Big_Integer;
  91       Remainder         : out Big_Integer;
  92       Discard_Quotient  : Boolean := False;
  93       Discard_Remainder : Boolean := False);
  94    --  Returns the Quotient and Remainder from dividing abs (X) by abs (Y). The
  95    --  values of X and Y are not modified. If Discard_Quotient is True, then
  96    --  Quotient is undefined on return, and if Discard_Remainder is True, then
  97    --  Remainder is undefined on return. Service routine for Big_Div/Rem/Mod.
  98
  99    function Normalize
 100      (X   : Digit_Vector;
 101       Neg : Boolean := False) return Big_Integer;
 102    --  Given a digit vector and sign, allocate and construct a big integer
 103    --  value. Note that X may have leading zeroes which must be removed, and if
 104    --  the result is zero, the sign is forced positive.
 105    --  If X is too big, Storage_Error is raised.
 106
 107    function "**" (X : Bignum; Y : SD) return Big_Integer;
 108    --  Exponentiation routine where we know right operand is one word
 109
 110    ---------
 111    -- Add --
 112    ---------
 113
 114    function Add
 115      (X, Y  : Digit_Vector;
 116       X_Neg : Boolean;
 117       Y_Neg : Boolean) return Big_Integer
 118    is
 119    begin
 120       --  If signs are the same, we are doing an addition, it is convenient to
 121       --  ensure that the first operand is the longer of the two.
 122
 123       if X_Neg = Y_Neg then
 124          if X'Last < Y'Last then
 125             return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg);
 126
 127          --  Here signs are the same, and the first operand is the longer
 128
 129          else
 130             pragma Assert (X_Neg = Y_Neg and then X'Last >= Y'Last);
 131
 132             --  Do addition, putting result in Sum (allowing for carry)
 133
 134             declare
 135                Sum : Digit_Vector (0 .. X'Last);
 136                RD  : DD;
 137
 138             begin
 139                RD := 0;
 140                for J in reverse 1 .. X'Last loop
 141                   RD := RD + DD (X (J));
 142
 143                   if J >= 1 + (X'Last - Y'Last) then
 144                      RD := RD + DD (Y (J - (X'Last - Y'Last)));
 145                   end if;
 146
 147                   Sum (J) := LSD (RD);
 148                   RD := RD / Base;
 149                end loop;
 150
 151                Sum (0) := SD (RD);
 152                return Normalize (Sum, X_Neg);
 153             end;
 154          end if;
 155
 156       --  Signs are different so really this is a subtraction, we want to make
 157       --  sure that the largest magnitude operand is the first one, and then
 158       --  the result will have the sign of the first operand.
 159
 160       else
 161          declare
 162             CR : constant Compare_Result := Compare (X, Y, False, False);
 163
 164          begin
 165             if CR = EQ then
 166                return Normalize (Zero_Data);
 167
 168             elsif CR = LT then
 169                return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg);
 170
 171             else
 172                pragma Assert (X_Neg /= Y_Neg and then CR = GT);
 173
 174                --  Do subtraction, putting result in Diff
 175
 176                declare
 177                   Diff : Digit_Vector (1 .. X'Length);
 178                   RD   : DD;
 179
 180                begin
 181                   RD := 0;
 182                   for J in reverse 1 .. X'Last loop
 183                      RD := RD + DD (X (J));
 184
 185                      if J >= 1 + (X'Last - Y'Last) then
 186                         RD := RD - DD (Y (J - (X'Last - Y'Last)));
 187                      end if;
 188
 189                      Diff (J) := LSD (RD);
 190                      RD := (if RD < Base then 0 else -1);
 191                   end loop;
 192
 193                   return Normalize (Diff, X_Neg);
 194                end;
 195             end if;
 196          end;
 197       end if;
 198    end Add;
 199
 200    -------------
 201    -- Big_Abs --
 202    -------------
 203
 204    function Big_Abs (X : Bignum) return Big_Integer is
 205    begin
 206       return Normalize (X.D);
 207    end Big_Abs;
 208
 209    -------------
 210    -- Big_Add --
 211    -------------
 212
 213    function Big_Add  (X, Y : Bignum) return Big_Integer is
 214    begin
 215       return Add (X.D, Y.D, X.Neg, Y.Neg);
 216    end Big_Add;
 217
 218    -------------
 219    -- Big_Div --
 220    -------------
 221
 222    --  This table is excerpted from RM 4.5.5(28-30) and shows how the result
 223    --  varies with the signs of the operands.
 224
 225    --   A      B   A/B      A     B    A/B
 226    --
 227    --   10     5    2      -10    5    -2
 228    --   11     5    2      -11    5    -2
 229    --   12     5    2      -12    5    -2
 230    --   13     5    2      -13    5    -2
 231    --   14     5    2      -14    5    -2
 232    --
 233    --   A      B   A/B      A     B    A/B
 234    --
 235    --   10    -5   -2      -10   -5     2
 236    --   11    -5   -2      -11   -5     2
 237    --   12    -5   -2      -12   -5     2
 238    --   13    -5   -2      -13   -5     2
 239    --   14    -5   -2      -14   -5     2
 240
 241    function Big_Div  (X, Y : Bignum) return Big_Integer is
 242       Q, R : aliased Big_Integer;
 243    begin
 244       Div_Rem (X, Y, Q, R, Discard_Remainder => True);
 245       To_Bignum (Q).Neg := To_Bignum (Q).Len > 0 and then (X.Neg xor Y.Neg);
 246       return Q;
 247    end Big_Div;
 248
 249    ----------
 250    -- "**" --
 251    ----------
 252
 253    function "**" (X : Bignum; Y : SD) return Big_Integer is
 254    begin
 255       case Y is
 256
 257          --  X ** 0 is 1
 258
 259          when 0 =>
 260             return Normalize (One_Data);
 261
 262          --  X ** 1 is X
 263
 264          when 1 =>
 265             return Normalize (X.D);
 266
 267          --  X ** 2 is X * X
 268
 269          when 2 =>
 270             return Big_Mul (X, X);
 271
 272          --  For X greater than 2, use the recursion
 273
 274          --  X even, X ** Y = (X ** (Y/2)) ** 2;
 275          --  X odd,  X ** Y = (X ** (Y/2)) ** 2 * X;
 276
 277          when others =>
 278             declare
 279                XY2  : aliased Big_Integer := X ** (Y / 2);
 280                XY2S : aliased Big_Integer :=
 281                  Big_Mul (To_Bignum (XY2), To_Bignum (XY2));
 282
 283             begin
 284                Free_Big_Integer (XY2);
 285
 286                if (Y and 1) = 0 then
 287                   return XY2S;
 288                else
 289                   return Res : constant Big_Integer :=
 290                     Big_Mul (To_Bignum (XY2S), X)
 291                   do
 292                      Free_Big_Integer (XY2S);
 293                   end return;
 294                end if;
 295             end;
 296       end case;
 297    end "**";
 298
 299    -------------
 300    -- Big_Exp --
 301    -------------
 302
 303    function Big_Exp  (X, Y : Bignum) return Big_Integer is
 304    begin
 305       --  Error if right operand negative
 306
 307       if Y.Neg then
 308          raise Constraint_Error with "exponentiation to negative power";
 309
 310       --  X ** 0 is always 1 (including 0 ** 0, so do this test first)
 311
 312       elsif Y.Len = 0 then
 313          return Normalize (One_Data);
 314
 315       --  0 ** X is always 0 (for X non-zero)
 316
 317       elsif X.Len = 0 then
 318          return Normalize (Zero_Data);
 319
 320       --  (+1) ** Y = 1
 321       --  (-1) ** Y = +/-1 depending on whether Y is even or odd
 322
 323       elsif X.Len = 1 and then X.D (1) = 1 then
 324          return Normalize
 325            (X.D, Neg => X.Neg and then (Y.D (Y.Len) and 1) = 1);
 326
 327       --  If the absolute value of the base is greater than 1, then the
 328       --  exponent must not be bigger than one word, otherwise the result
 329       --  is ludicrously large, and we just signal Storage_Error right away.
 330
 331       elsif Y.Len > 1 then
 332          raise Storage_Error with "exponentiation result is too large";
 333
 334       --  Special case (+/-)2 ** K, where K is 1 .. 31 using a shift
 335
 336       elsif X.Len = 1 and then X.D (1) = 2 and then Y.D (1) < 32 then
 337          declare
 338             D : constant Digit_Vector (1 .. 1) :=
 339                   [Shift_Left (SD'(1), Natural (Y.D (1)))];
 340          begin
 341             return Normalize (D, X.Neg);
 342          end;
 343
 344       --  Remaining cases have right operand of one word
 345
 346       else
 347          return X ** Y.D (1);
 348       end if;
 349    end Big_Exp;
 350
 351    -------------
 352    -- Big_And --
 353    -------------
 354
 355    function Big_And (X, Y : Bignum) return Big_Integer is
 356    begin
 357       if X.Len > Y.Len then
 358          return Big_And (X => Y, Y => X);
 359       end if;
 360
 361       --  X is the smallest integer
 362
 363       declare
 364          Result : Digit_Vector (1 .. X.Len);
 365          Diff   : constant Length := Y.Len - X.Len;
 366       begin
 367          for J in 1 .. X.Len loop
 368             Result (J) := X.D (J) and Y.D (J + Diff);
 369          end loop;
 370
 371          return Normalize (Result, X.Neg and Y.Neg);
 372       end;
 373    end Big_And;
 374
 375    ------------
 376    -- Big_Or --
 377    ------------
 378
 379    function Big_Or  (X, Y : Bignum) return Big_Integer is
 380    begin
 381       if X.Len < Y.Len then
 382          return Big_Or (X => Y, Y => X);
 383       end if;
 384
 385       --  X is the largest integer
 386
 387       declare
 388          Result : Digit_Vector (1 .. X.Len);
 389          Index  : Length;
 390          Diff   : constant Length := X.Len - Y.Len;
 391
 392       begin
 393          Index := 1;
 394
 395          while Index <= Diff loop
 396             Result (Index) := X.D (Index);
 397             Index := Index + 1;
 398          end loop;
 399
 400          for J in 1 .. Y.Len loop
 401             Result (Index) := X.D (Index) or Y.D (J);
 402             Index := Index + 1;
 403          end loop;
 404
 405          return Normalize (Result, X.Neg or Y.Neg);
 406       end;
 407    end Big_Or;
 408
 409    --------------------
 410    -- Big_Shift_Left --
 411    --------------------
 412
 413    function Big_Shift_Left (X : Bignum; Amount : Natural) return Big_Integer is
 414    begin
 415       if X.Neg then
 416          raise Constraint_Error;
 417       elsif Amount = 0 then
 418          return Allocate_Big_Integer (X.D, False);
 419       end if;
 420
 421       declare
 422          Shift  : constant Natural := Amount rem SD'Size;
 423          Result : Digit_Vector (0 .. X.Len + Amount / SD'Size);
 424          Carry  : SD := 0;
 425
 426       begin
 427          for J in X.Len + 1 .. Result'Last loop
 428             Result (J) := 0;
 429          end loop;
 430
 431          for J in reverse 1 .. X.Len loop
 432             Result (J) := Shift_Left (X.D (J), Shift) or Carry;
 433             Carry := Shift_Right (X.D (J), SD'Size - Shift);
 434          end loop;
 435
 436          Result (0) := Carry;
 437          return Normalize (Result, False);
 438       end;
 439    end Big_Shift_Left;
 440
 441    ---------------------
 442    -- Big_Shift_Right --
 443    ---------------------
 444
 445    function Big_Shift_Right
 446      (X : Bignum; Amount : Natural) return Big_Integer is
 447    begin
 448       if X.Neg then
 449          raise Constraint_Error;
 450       elsif Amount = 0 then
 451          return Allocate_Big_Integer (X.D, False);
 452       end if;
 453
 454       declare
 455          Shift  : constant Natural := Amount rem SD'Size;
 456          Result : Digit_Vector (1 .. X.Len - Amount / SD'Size);
 457          Carry  : SD := 0;
 458
 459       begin
 460          for J in 1 .. Result'Last - 1 loop
 461             Result (J) := Shift_Right (X.D (J), Shift) or Carry;
 462             Carry := Shift_Left (X.D (J), SD'Size - Shift);
 463          end loop;
 464
 465          Result (Result'Last) :=
 466            Shift_Right (X.D (Result'Last), Shift) or Carry;
 467
 468          return Normalize (Result, False);
 469       end;
 470    end Big_Shift_Right;
 471
 472    ------------
 473    -- Big_EQ --
 474    ------------
 475
 476    function Big_EQ (X, Y : Bignum) return Boolean is
 477    begin
 478       return Compare (X.D, Y.D, X.Neg, Y.Neg) = EQ;
 479    end Big_EQ;
 480
 481    ------------
 482    -- Big_GE --
 483    ------------
 484
 485    function Big_GE (X, Y : Bignum) return Boolean is
 486    begin
 487       return Compare (X.D, Y.D, X.Neg, Y.Neg) /= LT;
 488    end Big_GE;
 489
 490    ------------
 491    -- Big_GT --
 492    ------------
 493
 494    function Big_GT (X, Y : Bignum) return Boolean is
 495    begin
 496       return Compare (X.D, Y.D, X.Neg, Y.Neg) = GT;
 497    end Big_GT;
 498
 499    ------------
 500    -- Big_LE --
 501    ------------
 502
 503    function Big_LE (X, Y : Bignum) return Boolean is
 504    begin
 505       return Compare (X.D, Y.D, X.Neg, Y.Neg) /= GT;
 506    end Big_LE;
 507
 508    ------------
 509    -- Big_LT --
 510    ------------
 511
 512    function Big_LT (X, Y : Bignum) return Boolean is
 513    begin
 514       return Compare (X.D, Y.D, X.Neg, Y.Neg) = LT;
 515    end Big_LT;
 516
 517    -------------
 518    -- Big_Mod --
 519    -------------
 520
 521    --  This table is excerpted from RM 4.5.5(28-30) and shows how the result
 522    --  of Rem and Mod vary with the signs of the operands.
 523
 524    --   A      B    A mod B  A rem B     A     B    A mod B  A rem B
 525
 526    --   10     5       0        0       -10    5       0        0
 527    --   11     5       1        1       -11    5       4       -1
 528    --   12     5       2        2       -12    5       3       -2
 529    --   13     5       3        3       -13    5       2       -3
 530    --   14     5       4        4       -14    5       1       -4
 531
 532    --   A      B    A mod B  A rem B     A     B    A mod B  A rem B
 533
 534    --   10    -5       0        0       -10   -5       0        0
 535    --   11    -5      -4        1       -11   -5      -1       -1
 536    --   12    -5      -3        2       -12   -5      -2       -2
 537    --   13    -5      -2        3       -13   -5      -3       -3
 538    --   14    -5      -1        4       -14   -5      -4       -4
 539
 540    function Big_Mod (X, Y : Bignum) return Big_Integer is
 541       Q, R : aliased Big_Integer;
 542
 543    begin
 544       --  If signs are same, result is same as Rem
 545
 546       if X.Neg = Y.Neg then
 547          return Big_Rem (X, Y);
 548
 549       --  Case where Mod is different
 550
 551       else
 552          --  Do division
 553
 554          Div_Rem (X, Y, Q, R, Discard_Quotient => True);
 555
 556          --  Zero result is unchanged
 557
 558          if To_Bignum (R).Len = 0 then
 559             return R;
 560
 561          --  Otherwise adjust result
 562
 563          else
 564             declare
 565                T1 : aliased Big_Integer := Big_Sub (Y, To_Bignum (R));
 566             begin
 567                To_Bignum (T1).Neg := Y.Neg;
 568                Free_Big_Integer (R);
 569                return T1;
 570             end;
 571          end if;
 572       end if;
 573    end Big_Mod;
 574
 575    -------------
 576    -- Big_Mul --
 577    -------------
 578
 579    function Big_Mul (X, Y : Bignum) return Big_Integer is
 580       Result : Digit_Vector (1 .. X.Len + Y.Len) := [others => 0];
 581       --  Accumulate result (max length of result is sum of operand lengths)
 582
 583       L : Length;
 584       --  Current result digit
 585
 586       D : DD;
 587       --  Result digit
 588
 589    begin
 590       for J in 1 .. X.Len loop
 591          for K in 1 .. Y.Len loop
 592             L := Result'Last - (X.Len - J) - (Y.Len - K);
 593             D := DD (X.D (J)) * DD (Y.D (K)) + DD (Result (L));
 594             Result (L) := LSD (D);
 595             D := D / Base;
 596
 597             --  D is carry which must be propagated
 598
 599             while D /= 0 and then L >= 1 loop
 600                L := L - 1;
 601                D := D + DD (Result (L));
 602                Result (L) := LSD (D);
 603                D := D / Base;
 604             end loop;
 605
 606             --  Must not have a carry trying to extend max length
 607
 608             pragma Assert (D = 0);
 609          end loop;
 610       end loop;
 611
 612       --  Return result
 613
 614       return Normalize (Result, X.Neg xor Y.Neg);
 615    end Big_Mul;
 616
 617    ------------
 618    -- Big_NE --
 619    ------------
 620
 621    function Big_NE (X, Y : Bignum) return Boolean is
 622    begin
 623       return Compare (X.D, Y.D, X.Neg, Y.Neg) /= EQ;
 624    end Big_NE;
 625
 626    -------------
 627    -- Big_Neg --
 628    -------------
 629
 630    function Big_Neg (X : Bignum) return Big_Integer is
 631    begin
 632       return Normalize (X.D, not X.Neg);
 633    end Big_Neg;
 634
 635    -------------
 636    -- Big_Rem --
 637    -------------
 638
 639    --  This table is excerpted from RM 4.5.5(28-30) and shows how the result
 640    --  varies with the signs of the operands.
 641
 642    --   A      B   A rem B   A     B   A rem B
 643
 644    --   10     5      0     -10    5      0
 645    --   11     5      1     -11    5     -1
 646    --   12     5      2     -12    5     -2
 647    --   13     5      3     -13    5     -3
 648    --   14     5      4     -14    5     -4
 649
 650    --   A      B  A rem B    A     B   A rem B
 651
 652    --   10    -5     0      -10   -5      0
 653    --   11    -5     1      -11   -5     -1
 654    --   12    -5     2      -12   -5     -2
 655    --   13    -5     3      -13   -5     -3
 656    --   14    -5     4      -14   -5     -4
 657
 658    function Big_Rem (X, Y : Bignum) return Big_Integer is
 659       Q, R : aliased Big_Integer;
 660    begin
 661       Div_Rem (X, Y, Q, R, Discard_Quotient => True);
 662       To_Bignum (R).Neg := To_Bignum (R).Len > 0 and then X.Neg;
 663       return R;
 664    end Big_Rem;
 665
 666    -------------
 667    -- Big_Sub --
 668    -------------
 669
 670    function Big_Sub (X, Y : Bignum) return Big_Integer is
 671    begin
 672       --  If right operand zero, return left operand (avoiding sharing)
 673
 674       if Y.Len = 0 then
 675          return Normalize (X.D, X.Neg);
 676
 677       --  Otherwise add negative of right operand
 678
 679       else
 680          return Add (X.D, Y.D, X.Neg, not Y.Neg);
 681       end if;
 682    end Big_Sub;
 683
 684    -------------
 685    -- Compare --
 686    -------------
 687
 688    function Compare
 689      (X, Y         : Digit_Vector;
 690       X_Neg, Y_Neg : Boolean) return Compare_Result
 691    is
 692    begin
 693       --  Signs are different, that's decisive, since 0 is always plus
 694
 695       if X_Neg /= Y_Neg then
 696          return (if X_Neg then LT else GT);
 697
 698       --  Lengths are different, that's decisive since no leading zeroes
 699
 700       elsif X'Last /= Y'Last then
 701          return (if X'Last > Y'Last xor X_Neg then GT else LT);
 702
 703       --  Need to compare data
 704
 705       else
 706          for J in X'Range loop
 707             if X (J) /= Y (J) then
 708                return (if X (J) > Y (J) xor X_Neg then GT else LT);
 709             end if;
 710          end loop;
 711
 712          return EQ;
 713       end if;
 714    end Compare;
 715
 716    -------------
 717    -- Div_Rem --
 718    -------------
 719
 720    procedure Div_Rem
 721      (X, Y              : Bignum;
 722       Quotient          : out Big_Integer;
 723       Remainder         : out Big_Integer;
 724       Discard_Quotient  : Boolean := False;
 725       Discard_Remainder : Boolean := False) is
 726    begin
 727       --  Error if division by zero
 728
 729       if Y.Len = 0 then
 730          raise Constraint_Error with "division by zero";
 731       end if;
 732
 733       --  Handle simple cases with special tests
 734
 735       --  If X < Y then quotient is zero and remainder is X
 736
 737       if Compare (X.D, Y.D, False, False) = LT then
 738          if not Discard_Quotient then
 739             Quotient := Normalize (Zero_Data);
 740          end if;
 741
 742          if not Discard_Remainder then
 743             Remainder := Normalize (X.D);
 744          end if;
 745
 746          return;
 747
 748       --  If both X and Y are less than 2**63-1, we can use Long_Long_Integer
 749       --  arithmetic. Note it is good not to do an accurate range check against
 750       --  Long_Long_Integer since -2**63 / -1 overflows.
 751
 752       elsif (X.Len <= 1 or else (X.Len = 2 and then X.D (1) < 2**31))
 753               and then
 754             (Y.Len <= 1 or else (Y.Len = 2 and then Y.D (1) < 2**31))
 755       then
 756          declare
 757             A : constant LLI := abs (From_Bignum (X));
 758             B : constant LLI := abs (From_Bignum (Y));
 759          begin
 760             if not Discard_Quotient then
 761                Quotient := To_Bignum (A / B);
 762             end if;
 763
 764             if not Discard_Remainder then
 765                Remainder := To_Bignum (A rem B);
 766             end if;
 767
 768             return;
 769          end;
 770
 771       --  Easy case if divisor is one digit
 772
 773       elsif Y.Len = 1 then
 774          declare
 775             ND  : DD;
 776             Div : constant DD := DD (Y.D (1));
 777
 778             Result : Digit_Vector (1 .. X.Len);
 779             Remdr  : Digit_Vector (1 .. 1);
 780
 781          begin
 782             ND := 0;
 783             for J in 1 .. X.Len loop
 784                ND := Base * ND + DD (X.D (J));
 785                pragma Assert (Div /= 0);
 786                Result (J) := SD (ND / Div);
 787                ND := ND rem Div;
 788             end loop;
 789
 790             if not Discard_Quotient then
 791                Quotient  := Normalize (Result);
 792             end if;
 793
 794             if not Discard_Remainder then
 795                Remdr (1) := SD (ND);
 796                Remainder := Normalize (Remdr);
 797             end if;
 798
 799             return;
 800          end;
 801       end if;
 802
 803       --  The complex full multi-precision case. We will employ algorithm
 804       --  D defined in the section "The Classical Algorithms" (sec. 4.3.1)
 805       --  of Donald Knuth's "The Art of Computer Programming", Vol. 2, 2nd
 806       --  edition. The terminology is adjusted for this section to match that
 807       --  reference.
 808
 809       --  We are dividing X.Len digits of X (called u here) by Y.Len digits
 810       --  of Y (called v here), developing the quotient and remainder. The
 811       --  numbers are represented using Base, which was chosen so that we have
 812       --  the operations of multiplying to single digits (SD) to form a double
 813       --  digit (DD), and dividing a double digit (DD) by a single digit (SD)
 814       --  to give a single digit quotient and a single digit remainder.
 815
 816       --  Algorithm D from Knuth
 817
 818       --  Comments here with square brackets are directly from Knuth
 819
 820       Algorithm_D : declare
 821
 822          --  The following lower case variables correspond exactly to the
 823          --  terminology used in algorithm D.
 824
 825          m : constant Length := X.Len - Y.Len;
 826          n : constant Length := Y.Len;
 827          b : constant DD     := Base;
 828
 829          u : Digit_Vector (0 .. m + n);
 830          v : Digit_Vector (1 .. n);
 831          q : Digit_Vector (0 .. m);
 832          r : Digit_Vector (1 .. n);
 833
 834          u0 : SD renames u (0);
 835          v1 : SD renames v (1);
 836          v2 : SD renames v (2);
 837
 838          d    : DD;
 839          j    : Length;
 840          qhat : DD;
 841          rhat : DD;
 842          temp : DD;
 843
 844       begin
 845          --  Initialize data of left and right operands
 846
 847          for J in 1 .. m + n loop
 848             u (J) := X.D (J);
 849          end loop;
 850
 851          for J in 1 .. n loop
 852             v (J) := Y.D (J);
 853          end loop;
 854
 855          --  [Division of nonnegative integers.] Given nonnegative integers u
 856          --  = (ul,u2..um+n) and v = (v1,v2..vn), where v1 /= 0 and n > 1, we
 857          --  form the quotient u / v = (q0,ql..qm) and the remainder u mod v =
 858          --  (r1,r2..rn).
 859
 860          pragma Assert (v1 /= 0);
 861          pragma Assert (n > 1);
 862
 863          --  Dl. [Normalize.] Set d = b/(vl + 1). Then set (u0,u1,u2..um+n)
 864          --  equal to (u1,u2..um+n) times d, and set (v1,v2..vn) equal to
 865          --  (v1,v2..vn) times d. Note the introduction of a new digit position
 866          --  u0 at the left of u1; if d = 1 all we need to do in this step is
 867          --  to set u0 = 0.
 868
 869          d := b / (DD (v1) + 1);
 870
 871          if d = 1 then
 872             u0 := 0;
 873
 874          else
 875             declare
 876                Carry : DD;
 877                Tmp   : DD;
 878
 879             begin
 880                --  Multiply Dividend (u) by d
 881
 882                Carry := 0;
 883                for J in reverse 1 .. m + n loop
 884                   Tmp   := DD (u (J)) * d + Carry;
 885                   u (J) := LSD (Tmp);
 886                   Carry := Tmp / Base;
 887                end loop;
 888
 889                u0 := SD (Carry);
 890
 891                --  Multiply Divisor (v) by d
 892
 893                Carry := 0;
 894                for J in reverse 1 .. n loop
 895                   Tmp   := DD (v (J)) * d + Carry;
 896                   v (J) := LSD (Tmp);
 897                   Carry := Tmp / Base;
 898                end loop;
 899
 900                pragma Assert (Carry = 0);
 901             end;
 902          end if;
 903
 904          --  D2. [Initialize j.] Set j = 0. The loop on j, steps D2 through D7,
 905          --  will be essentially a division of (uj, uj+1..uj+n) by (v1,v2..vn)
 906          --  to get a single quotient digit qj.
 907
 908          j := 0;
 909
 910          --  Loop through digits
 911
 912          loop
 913             --  Note: In the original printing, step D3 was as follows:
 914
 915             --  D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise
 916             --  set qhat to (uj,uj+1)/v1. Now test if v2 * qhat is greater than
 917             --  (uj*b + uj+1 - qhat*v1)*b + uj+2. If so, decrease qhat by 1 and
 918             --  repeat this test
 919
 920             --  This had a bug not discovered till 1995, see Vol 2 errata:
 921             --  http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz. Under
 922             --  rare circumstances the expression in the test could overflow.
 923             --  This version was further corrected in 2005, see Vol 2 errata:
 924             --  http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz.
 925             --  The code below is the fixed version of this step.
 926
 927             --  D3. [Calculate qhat.] Set qhat to (uj,uj+1)/v1 and rhat to
 928             --  to (uj,uj+1) mod v1.
 929
 930             temp := u (j) & u (j + 1);
 931             qhat := temp / DD (v1);
 932             rhat := temp mod DD (v1);
 933
 934             --  D3 (continued). Now test if qhat >= b or v2*qhat > (rhat,uj+2):
 935             --  if so, decrease qhat by 1, increase rhat by v1, and repeat this
 936             --  test if rhat < b. [The test on v2 determines at high speed
 937             --  most of the cases in which the trial value qhat is one too
 938             --  large, and eliminates all cases where qhat is two too large.]
 939
 940             while qhat >= b
 941               or else DD (v2) * qhat > LSD (rhat) & u (j + 2)
 942             loop
 943                qhat := qhat - 1;
 944                rhat := rhat + DD (v1);
 945                exit when rhat >= b;
 946             end loop;
 947
 948             --  D4. [Multiply and subtract.] Replace (uj,uj+1..uj+n) by
 949             --  (uj,uj+1..uj+n) minus qhat times (v1,v2..vn). This step
 950             --  consists of a simple multiplication by a one-place number,
 951             --  combined with a subtraction.
 952
 953             --  The digits (uj,uj+1..uj+n) are always kept positive; if the
 954             --  result of this step is actually negative then (uj,uj+1..uj+n)
 955             --  is left as the true value plus b**(n+1), i.e. as the b's
 956             --  complement of the true value, and a "borrow" to the left is
 957             --  remembered.
 958
 959             declare
 960                Borrow : SD;
 961                Carry  : DD;
 962                Temp   : DD;
 963
 964                Negative : Boolean;
 965                --  Records if subtraction causes a negative result, requiring
 966                --  an add back (case where qhat turned out to be 1 too large).
 967
 968             begin
 969                Borrow := 0;
 970                for K in reverse 1 .. n loop
 971                   Temp := qhat * DD (v (K)) + DD (Borrow);
 972                   Borrow := MSD (Temp);
 973
 974                   if LSD (Temp) > u (j + K) then
 975                      Borrow := Borrow + 1;
 976                   end if;
 977
 978                   u (j + K) := u (j + K) - LSD (Temp);
 979                end loop;
 980
 981                Negative := u (j) < Borrow;
 982                u (j) := u (j) - Borrow;
 983
 984                --  D5. [Test remainder.] Set qj = qhat. If the result of step
 985                --  D4 was negative, we will do the add back step (step D6).
 986
 987                q (j) := LSD (qhat);
 988
 989                if Negative then
 990
 991                   --  D6. [Add back.] Decrease qj by 1, and add (0,v1,v2..vn)
 992                   --  to (uj,uj+1,uj+2..uj+n). (A carry will occur to the left
 993                   --  of uj, and it is be ignored since it cancels with the
 994                   --  borrow that occurred in D4.)
 995
 996                   q (j) := q (j) - 1;
 997
 998                   Carry := 0;
 999                   for K in reverse 1 .. n loop
1000                      Temp := DD (v (K)) + DD (u (j + K)) + Carry;
1001                      u (j + K) := LSD (Temp);
1002                      Carry := Temp / Base;
1003                   end loop;
1004
1005                   u (j) := u (j) + SD (Carry);
1006                end if;
1007             end;
1008
1009             --  D7. [Loop on j.] Increase j by one. Now if j <= m, go back to
1010             --  D3 (the start of the loop on j).
1011
1012             j := j + 1;
1013             exit when not (j <= m);
1014          end loop;
1015
1016          --  D8. [Unnormalize.] Now (qo,ql..qm) is the desired quotient, and
1017          --  the desired remainder may be obtained by dividing (um+1..um+n)
1018          --  by d.
1019
1020          if not Discard_Quotient then
1021             Quotient := Normalize (q);
1022          end if;
1023
1024          if not Discard_Remainder then
1025             declare
1026                Remdr : DD;
1027             begin
1028                Remdr := 0;
1029
1030                for K in 1 .. n loop
1031                   Remdr := Base * Remdr + DD (u (m + K));
1032                   r (K) := SD (Remdr / d);
1033                   Remdr := Remdr rem d;
1034                end loop;
1035
1036                pragma Assert (Remdr = 0);
1037             end;
1038
1039             Remainder := Normalize (r);
1040          end if;
1041       end Algorithm_D;
1042    end Div_Rem;
1043
1044    -----------------
1045    -- From_Bignum --
1046    -----------------
1047
1048    function From_Bignum (X : Bignum) return Long_Long_Long_Integer is
1049    begin
1050       if X.Len = 0 then
1051          return 0;
1052
1053       elsif X.Len = 1 then
1054          return (if X.Neg then -LLLI (X.D (1)) else LLLI (X.D (1)));
1055
1056       elsif X.Len = 2 then
1057          declare
1058             Mag : constant DD := X.D (1) & X.D (2);
1059          begin
1060             if X.Neg and then (Mag <= 2 ** 63 or else LLLI_Is_128) then
1061                return -LLLI (Mag);
1062             elsif Mag < 2 ** 63 or else LLLI_Is_128 then
1063                return LLLI (Mag);
1064             end if;
1065          end;
1066
1067       elsif X.Len = 3 and then LLLI_Is_128 then
1068          declare
1069             Hi  : constant SD := X.D (1);
1070             Lo  : constant DD := X.D (2) & X.D (3);
1071             Mag : constant Unsigned_128 :=
1072                     Shift_Left (Unsigned_128 (Hi), 64) + Unsigned_128 (Lo);
1073          begin
1074             return (if X.Neg then -LLLI (Mag) else LLLI (Mag));
1075          end;
1076
1077       elsif X.Len = 4 and then LLLI_Is_128 then
1078          declare
1079             Hi  : constant DD := X.D (1) & X.D (2);
1080             Lo  : constant DD := X.D (3) & X.D (4);
1081             Mag : constant Unsigned_128 :=
1082                     Shift_Left (Unsigned_128 (Hi), 64) + Unsigned_128 (Lo);
1083          begin
1084             if X.Neg
1085               and then (Hi < 2 ** 63 or else (Hi = 2 ** 63 and then Lo = 0))
1086             then
1087                return -LLLI (Mag);
1088             elsif Hi < 2 ** 63 then
1089                return LLLI (Mag);
1090             end if;
1091          end;
1092       end if;
1093
1094       raise Constraint_Error with "expression value out of range";
1095    end From_Bignum;
1096
1097    function From_Bignum (X : Bignum) return Long_Long_Integer is
1098    begin
1099       return Long_Long_Integer (Long_Long_Long_Integer'(From_Bignum (X)));
1100    end From_Bignum;
1101
1102    function From_Bignum (X : Bignum) return Unsigned_128 is
1103    begin
1104       if X.Neg then
1105          null;
1106
1107       elsif X.Len = 0 then
1108          return 0;
1109
1110       elsif X.Len = 1 then
1111          return Unsigned_128 (X.D (1));
1112
1113       elsif X.Len = 2 then
1114          return Unsigned_128 (DD'(X.D (1) & X.D (2)));
1115
1116       elsif X.Len = 3 and then LLLI_Is_128 then
1117          return
1118            Shift_Left (Unsigned_128 (X.D (1)), 64) +
1119              Unsigned_128 (DD'(X.D (2) & X.D (3)));
1120
1121       elsif X.Len = 4 and then LLLI_Is_128 then
1122          return
1123            Shift_Left (Unsigned_128 (DD'(X.D (1) & X.D (2))), 64) +
1124              Unsigned_128 (DD'(X.D (3) & X.D (4)));
1125       end if;
1126
1127       raise Constraint_Error with "expression value out of range";
1128    end From_Bignum;
1129
1130    function From_Bignum (X : Bignum) return Unsigned_64 is
1131    begin
1132       return Unsigned_64 (Unsigned_128'(From_Bignum (X)));
1133    end From_Bignum;
1134
1135    -------------------------
1136    -- Bignum_In_LLI_Range --
1137    -------------------------
1138
1139    function Bignum_In_LLI_Range (X : Bignum) return Boolean is
1140    begin
1141       --  If length is 0 or 1, definitely fits
1142
1143       if X.Len <= 1 then
1144          return True;
1145
1146       --  If length is greater than 2, definitely does not fit
1147
1148       elsif X.Len > 2 then
1149          return False;
1150
1151       --  Length is 2, more tests needed
1152
1153       else
1154          declare
1155             Mag : constant DD := X.D (1) & X.D (2);
1156          begin
1157             return Mag < 2 ** 63 or else (X.Neg and then Mag = 2 ** 63);
1158          end;
1159       end if;
1160    end Bignum_In_LLI_Range;
1161
1162    ---------------
1163    -- Normalize --
1164    ---------------
1165
1166    Bignum_Limit : constant := 200;
1167
1168    function Normalize
1169      (X   : Digit_Vector;
1170       Neg : Boolean := False) return Big_Integer
1171    is
1172       J : Length;
1173
1174    begin
1175       J := X'First;
1176       while J <= X'Last and then X (J) = 0 loop
1177          J := J + 1;
1178       end loop;
1179
1180       if X'Last - J > Bignum_Limit then
1181          raise Storage_Error with "big integer limit exceeded";
1182       end if;
1183
1184       return Allocate_Big_Integer (X (J .. X'Last), J <= X'Last and then Neg);
1185    end Normalize;
1186
1187    ---------------
1188    -- To_Bignum --
1189    ---------------
1190
1191    function To_Bignum (X : Long_Long_Long_Integer) return Big_Integer is
1192
1193       function Convert_128
1194         (X : Long_Long_Long_Integer; Neg : Boolean) return Big_Integer;
1195       --  Convert a 128 bits natural integer to a Big_Integer
1196
1197       -----------------
1198       -- Convert_128 --
1199       -----------------
1200
1201       function Convert_128
1202         (X : Long_Long_Long_Integer; Neg : Boolean) return Big_Integer
1203       is
1204          Vector : Digit_Vector (1 .. 4);
1205          High   : constant Unsigned_64 :=
1206            Unsigned_64 (Shift_Right (Unsigned_128 (X), 64));
1207          Low    : constant Unsigned_64 :=
1208            Unsigned_64 (Unsigned_128 (X) and 16#FFFF_FFFF_FFFF_FFFF#);
1209
1210       begin
1211          Vector (1) := SD (High / Base);
1212          Vector (2) := SD (High mod Base);
1213          Vector (3) := SD (Low / Base);
1214          Vector (4) := SD (Low mod Base);
1215          return Normalize (Vector, Neg);
1216       end Convert_128;
1217
1218    begin
1219       if X = 0 then
1220          return Allocate_Big_Integer ([], False);
1221
1222       --  One word result
1223
1224       elsif X in -(2 ** 32 - 1) .. +(2 ** 32 - 1) then
1225          return Allocate_Big_Integer ([SD (abs X)], X < 0);
1226
1227       --  Large negative number annoyance
1228
1229       elsif X = -2 ** 63 then
1230          return Allocate_Big_Integer ([2 ** 31, 0], True);
1231
1232       elsif LLLI_Is_128 and then X = Long_Long_Long_Integer'First then
1233          return Allocate_Big_Integer ([2 ** 31, 0, 0, 0], True);
1234
1235       --  Other negative numbers
1236
1237       elsif X < 0 then
1238          if LLLI_Is_128 then
1239             return Convert_128 (-X, True);
1240          else
1241             return Allocate_Big_Integer
1242                      ((SD ((-X) / Base), SD ((-X) mod Base)), True);
1243          end if;
1244
1245       --  Positive numbers
1246
1247       else
1248          if LLLI_Is_128 then
1249             return Convert_128 (X, False);
1250          else
1251             return Allocate_Big_Integer
1252                      ((SD (X / Base), SD (X mod Base)), False);
1253          end if;
1254       end if;
1255    end To_Bignum;
1256
1257    function To_Bignum (X : Long_Long_Integer) return Big_Integer is
1258    begin
1259       return To_Bignum (Long_Long_Long_Integer (X));
1260    end To_Bignum;
1261
1262    function To_Bignum (X : Unsigned_128) return Big_Integer is
1263    begin
1264       if X = 0 then
1265          return Allocate_Big_Integer ([], False);
1266
1267       --  One word result
1268
1269       elsif X < 2 ** 32 then
1270          return Allocate_Big_Integer ([SD (X)], False);
1271
1272       --  Two word result
1273
1274       elsif Shift_Right (X, 32) < 2 ** 32 then
1275          return Allocate_Big_Integer ([SD (X / Base), SD (X mod Base)], False);
1276
1277       --  Three or four word result
1278
1279       else
1280          declare
1281             Vector : Digit_Vector (1 .. 4);
1282             High   : constant Unsigned_64 := Unsigned_64 (Shift_Right (X, 64));
1283             Low    : constant Unsigned_64 :=
1284               Unsigned_64 (X and 16#FFFF_FFFF_FFFF_FFFF#);
1285
1286          begin
1287             Vector (1) := SD (High / Base);
1288             Vector (2) := SD (High mod Base);
1289             Vector (3) := SD (Low / Base);
1290             Vector (4) := SD (Low mod Base);
1291             return Normalize (Vector, False);
1292          end;
1293       end if;
1294    end To_Bignum;
1295
1296    function To_Bignum (X : Unsigned_64) return Big_Integer is
1297    begin
1298       return To_Bignum (Unsigned_128 (X));
1299    end To_Bignum;
1300
1301    ---------------
1302    -- To_String --
1303    ---------------
1304
1305    Hex_Chars : constant array (0 .. 15) of Character := "0123456789ABCDEF";
1306
1307    function To_String
1308      (X : Bignum; Width : Natural := 0; Base : Positive := 10) return String
1309    is
1310       Big_Base : aliased Bignum_Data := (1, False, [SD (Base)]);
1311
1312       function Add_Base (S : String) return String;
1313       --  Add base information if Base /= 10
1314
1315       function Leading_Padding
1316         (Str        : String;
1317          Min_Length : Natural;
1318          Char       : Character := ' ') return String;
1319       --  Return padding of Char concatenated with Str so that the resulting
1320       --  string is at least Min_Length long.
1321
1322       function Image (Arg : Bignum) return String;
1323       --  Return image of Arg, assuming Arg is positive.
1324
1325       function Image (N : Natural) return String;
1326       --  Return image of N, with no leading space.
1327
1328       --------------
1329       -- Add_Base --
1330       --------------
1331
1332       function Add_Base (S : String) return String is
1333       begin
1334          if Base = 10 then
1335             return S;
1336          else
1337             return Image (Base) & "#" & S & "#";
1338          end if;
1339       end Add_Base;
1340
1341       -----------
1342       -- Image --
1343       -----------
1344
1345       function Image (N : Natural) return String is
1346          S : constant String := Natural'Image (N);
1347       begin
1348          return S (2 .. S'Last);
1349       end Image;
1350
1351       function Image (Arg : Bignum) return String is
1352       begin
1353          if Big_LT (Arg, Big_Base'Unchecked_Access) then
1354             return [Hex_Chars (Natural (LLI'(From_Bignum (Arg))))];
1355          else
1356             declare
1357                Div    : aliased Big_Integer;
1358                Remain : aliased Big_Integer;
1359                R      : Natural;
1360
1361             begin
1362                Div_Rem (Arg, Big_Base'Unchecked_Access, Div, Remain);
1363                R := Natural (LLI'(From_Bignum (To_Bignum (Remain))));
1364                Free_Big_Integer (Remain);
1365
1366                return S : constant String :=
1367                  Image (To_Bignum (Div)) & Hex_Chars (R)
1368                do
1369                   Free_Big_Integer (Div);
1370                end return;
1371             end;
1372          end if;
1373       end Image;
1374
1375       ---------------------
1376       -- Leading_Padding --
1377       ---------------------
1378
1379       function Leading_Padding
1380         (Str        : String;
1381          Min_Length : Natural;
1382          Char       : Character := ' ') return String is
1383       begin
1384          return [1 .. Integer'Max (Integer (Min_Length) - Str'Length, 0)
1385                         => Char] & Str;
1386       end Leading_Padding;
1387
1388       Zero : aliased Bignum_Data := (0, False, D => Zero_Data);
1389
1390    begin
1391       if Big_LT (X, Zero'Unchecked_Access) then
1392          declare
1393             X_Pos : aliased Bignum_Data := (X.Len, not X.Neg, X.D);
1394          begin
1395             return Leading_Padding
1396                      ("-" & Add_Base (Image (X_Pos'Unchecked_Access)), Width);
1397          end;
1398       else
1399          return Leading_Padding (" " & Add_Base (Image (X)), Width);
1400       end if;
1401    end To_String;
1402
1403    -------------
1404    -- Is_Zero --
1405    -------------
1406
1407    function Is_Zero (X : Bignum) return Boolean is
1408      (X /= null and then X.D = Zero_Data);
1409
1410 end System.Generic_Bignums;