gcc/ada/s-arit64.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUN-TIME COMPONENTS                         --
   4 --                                                                          --
   5 --                      S Y S T E M . A R I T H _ 6 4                       --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --          Copyright (C) 1992-2004 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 2,  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.  See the GNU General Public License --
  17 -- for  more details.  You should have  received  a copy of the GNU General --
  18 -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  19 -- to  the Free Software Foundation,  59 Temple Place - Suite 330,  Boston, --
  20 -- MA 02111-1307, USA.                                                      --
  21 --                                                                          --
  22 -- As a special exception,  if other files  instantiate  generics from this --
  23 -- unit, or you link  this unit with other files  to produce an executable, --
  24 -- this  unit  does not  by itself cause  the resulting  executable  to  be --
  25 -- covered  by the  GNU  General  Public  License.  This exception does not --
  26 -- however invalidate  any other reasons why  the executable file  might be --
  27 -- covered by the  GNU Public License.                                      --
  28 --                                                                          --
  29 -- GNAT was originally developed  by the GNAT team at  New York University. --
  30 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  31 --                                                                          --
  32 ------------------------------------------------------------------------------
  33
  34 with System.Pure_Exceptions; use System.Pure_Exceptions;
  35
  36 with Interfaces; use Interfaces;
  37 with Unchecked_Conversion;
  38
  39 package body System.Arith_64 is
  40
  41    pragma Suppress (Overflow_Check);
  42    pragma Suppress (Range_Check);
  43
  44    subtype Uns64 is Unsigned_64;
  45    function To_Uns is new Unchecked_Conversion (Int64, Uns64);
  46    function To_Int is new Unchecked_Conversion (Uns64, Int64);
  47
  48    subtype Uns32 is Unsigned_32;
  49
  50    -----------------------
  51    -- Local Subprograms --
  52    -----------------------
  53
  54    function "+" (A, B : Uns32) return Uns64;
  55    function "+" (A : Uns64; B : Uns32) return Uns64;
  56    pragma Inline ("+");
  57    --  Length doubling additions
  58
  59    function "*" (A, B : Uns32) return Uns64;
  60    pragma Inline ("*");
  61    --  Length doubling multiplication
  62
  63    function "/" (A : Uns64; B : Uns32) return Uns64;
  64    pragma Inline ("/");
  65    --  Length doubling division
  66
  67    function "rem" (A : Uns64; B : Uns32) return Uns64;
  68    pragma Inline ("rem");
  69    --  Length doubling remainder
  70
  71    function "&" (Hi, Lo : Uns32) return Uns64;
  72    pragma Inline ("&");
  73    --  Concatenate hi, lo values to form 64-bit result
  74
  75    function Le3 (X1, X2, X3 : Uns32; Y1, Y2, Y3 : Uns32) return Boolean;
  76    --  Determines if 96 bit value X1&X2&X3 <= Y1&Y2&Y3
  77
  78    function Lo (A : Uns64) return Uns32;
  79    pragma Inline (Lo);
  80    --  Low order half of 64-bit value
  81
  82    function Hi (A : Uns64) return Uns32;
  83    pragma Inline (Hi);
  84    --  High order half of 64-bit value
  85
  86    procedure Sub3 (X1, X2, X3 : in out Uns32; Y1, Y2, Y3 : in Uns32);
  87    --  Computes X1&X2&X3 := X1&X2&X3 - Y1&Y1&Y3 with mod 2**96 wrap
  88
  89    function To_Neg_Int (A : Uns64) return Int64;
  90    --  Convert to negative integer equivalent. If the input is in the range
  91    --  0 .. 2 ** 63, then the corresponding negative signed integer (obtained
  92    --  by negating the given value) is returned, otherwise constraint error
  93    --  is raised.
  94
  95    function To_Pos_Int (A : Uns64) return Int64;
  96    --  Convert to positive integer equivalent. If the input is in the range
  97    --  0 .. 2 ** 63-1, then the corresponding non-negative signed integer is
  98    --  returned, otherwise constraint error is raised.
  99
 100    procedure Raise_Error;
 101    pragma No_Return (Raise_Error);
 102    --  Raise constraint error with appropriate message
 103
 104    ---------
 105    -- "&" --
 106    ---------
 107
 108    function "&" (Hi, Lo : Uns32) return Uns64 is
 109    begin
 110       return Shift_Left (Uns64 (Hi), 32) or Uns64 (Lo);
 111    end "&";
 112
 113    ---------
 114    -- "*" --
 115    ---------
 116
 117    function "*" (A, B : Uns32) return Uns64 is
 118    begin
 119       return Uns64 (A) * Uns64 (B);
 120    end "*";
 121
 122    ---------
 123    -- "+" --
 124    ---------
 125
 126    function "+" (A, B : Uns32) return Uns64 is
 127    begin
 128       return Uns64 (A) + Uns64 (B);
 129    end "+";
 130
 131    function "+" (A : Uns64; B : Uns32) return Uns64 is
 132    begin
 133       return A + Uns64 (B);
 134    end "+";
 135
 136    ---------
 137    -- "/" --
 138    ---------
 139
 140    function "/" (A : Uns64; B : Uns32) return Uns64 is
 141    begin
 142       return A / Uns64 (B);
 143    end "/";
 144
 145    -----------
 146    -- "rem" --
 147    -----------
 148
 149    function "rem" (A : Uns64; B : Uns32) return Uns64 is
 150    begin
 151       return A rem Uns64 (B);
 152    end "rem";
 153
 154    --------------------------
 155    -- Add_With_Ovflo_Check --
 156    --------------------------
 157
 158    function Add_With_Ovflo_Check (X, Y : Int64) return Int64 is
 159       R : constant Int64 := To_Int (To_Uns (X) + To_Uns (Y));
 160
 161    begin
 162       if X >= 0 then
 163          if Y < 0 or else R >= 0 then
 164             return R;
 165          end if;
 166
 167       else -- X < 0
 168          if Y > 0 or else R < 0 then
 169             return R;
 170          end if;
 171       end if;
 172
 173       Raise_Error;
 174    end Add_With_Ovflo_Check;
 175
 176    -------------------
 177    -- Double_Divide --
 178    -------------------
 179
 180    procedure Double_Divide
 181      (X, Y, Z : Int64;
 182       Q, R    : out Int64;
 183       Round   : Boolean)
 184    is
 185       Xu  : constant Uns64 := To_Uns (abs X);
 186       Yu  : constant Uns64 := To_Uns (abs Y);
 187
 188       Yhi : constant Uns32 := Hi (Yu);
 189       Ylo : constant Uns32 := Lo (Yu);
 190
 191       Zu  : constant Uns64 := To_Uns (abs Z);
 192       Zhi : constant Uns32 := Hi (Zu);
 193       Zlo : constant Uns32 := Lo (Zu);
 194
 195       T1, T2     : Uns64;
 196       Du, Qu, Ru : Uns64;
 197       Den_Pos    : Boolean;
 198
 199    begin
 200       if Yu = 0 or else Zu = 0 then
 201          Raise_Error;
 202       end if;
 203
 204       --  Compute Y * Z. Note that if the result overflows 64 bits unsigned,
 205       --  then the rounded result is clearly zero (since the dividend is at
 206       --  most 2**63 - 1, the extra bit of precision is nice here!)
 207
 208       if Yhi /= 0 then
 209          if Zhi /= 0 then
 210             Q := 0;
 211             R := X;
 212             return;
 213          else
 214             T2 := Yhi * Zlo;
 215          end if;
 216
 217       else
 218          if Zhi /= 0 then
 219             T2 := Ylo * Zhi;
 220          else
 221             T2 := 0;
 222          end if;
 223       end if;
 224
 225       T1 := Ylo * Zlo;
 226       T2 := T2 + Hi (T1);
 227
 228       if Hi (T2) /= 0 then
 229          Q := 0;
 230          R := X;
 231          return;
 232       end if;
 233
 234       Du := Lo (T2) & Lo (T1);
 235       Qu := Xu / Du;
 236       Ru := Xu rem Du;
 237
 238       --  Deal with rounding case
 239
 240       if Round and then Ru > (Du - Uns64'(1)) / Uns64'(2) then
 241          Qu := Qu + Uns64'(1);
 242       end if;
 243
 244       --  Set final signs (RM 4.5.5(27-30))
 245
 246       Den_Pos := (Y < 0) = (Z < 0);
 247
 248       --  Case of dividend (X) sign positive
 249
 250       if X >= 0 then
 251          R := To_Int (Ru);
 252
 253          if Den_Pos then
 254             Q := To_Int (Qu);
 255          else
 256             Q := -To_Int (Qu);
 257          end if;
 258
 259       --  Case of dividend (X) sign negative
 260
 261       else
 262          R := -To_Int (Ru);
 263
 264          if Den_Pos then
 265             Q := -To_Int (Qu);
 266          else
 267             Q := To_Int (Qu);
 268          end if;
 269       end if;
 270    end Double_Divide;
 271
 272    --------
 273    -- Hi --
 274    --------
 275
 276    function Hi (A : Uns64) return Uns32 is
 277    begin
 278       return Uns32 (Shift_Right (A, 32));
 279    end Hi;
 280
 281    ---------
 282    -- Le3 --
 283    ---------
 284
 285    function Le3 (X1, X2, X3 : Uns32; Y1, Y2, Y3 : Uns32) return Boolean is
 286    begin
 287       if X1 < Y1 then
 288          return True;
 289       elsif X1 > Y1 then
 290          return False;
 291       elsif X2 < Y2 then
 292          return True;
 293       elsif X2 > Y2 then
 294          return False;
 295       else
 296          return X3 <= Y3;
 297       end if;
 298    end Le3;
 299
 300    --------
 301    -- Lo --
 302    --------
 303
 304    function Lo (A : Uns64) return Uns32 is
 305    begin
 306       return Uns32 (A and 16#FFFF_FFFF#);
 307    end Lo;
 308
 309    -------------------------------
 310    -- Multiply_With_Ovflo_Check --
 311    -------------------------------
 312
 313    function Multiply_With_Ovflo_Check (X, Y : Int64) return Int64 is
 314       Xu  : constant Uns64 := To_Uns (abs X);
 315       Xhi : constant Uns32 := Hi (Xu);
 316       Xlo : constant Uns32 := Lo (Xu);
 317
 318       Yu  : constant Uns64 := To_Uns (abs Y);
 319       Yhi : constant Uns32 := Hi (Yu);
 320       Ylo : constant Uns32 := Lo (Yu);
 321
 322       T1, T2 : Uns64;
 323
 324    begin
 325       if Xhi /= 0 then
 326          if Yhi /= 0 then
 327             Raise_Error;
 328          else
 329             T2 := Xhi * Ylo;
 330          end if;
 331
 332       elsif Yhi /= 0 then
 333          T2 := Xlo * Yhi;
 334
 335       else -- Yhi = Xhi = 0
 336          T2 := 0;
 337       end if;
 338
 339       --  Here we have T2 set to the contribution to the upper half
 340       --  of the result from the upper halves of the input values.
 341
 342       T1 := Xlo * Ylo;
 343       T2 := T2 + Hi (T1);
 344
 345       if Hi (T2) /= 0 then
 346          Raise_Error;
 347       end if;
 348
 349       T2 := Lo (T2) & Lo (T1);
 350
 351       if X >= 0 then
 352          if Y >= 0 then
 353             return To_Pos_Int (T2);
 354          else
 355             return To_Neg_Int (T2);
 356          end if;
 357       else -- X < 0
 358          if Y < 0 then
 359             return To_Pos_Int (T2);
 360          else
 361             return To_Neg_Int (T2);
 362          end if;
 363       end if;
 364
 365    end Multiply_With_Ovflo_Check;
 366
 367    -----------------
 368    -- Raise_Error --
 369    -----------------
 370
 371    procedure Raise_Error is
 372    begin
 373       Raise_Exception (CE, "64-bit arithmetic overflow");
 374    end Raise_Error;
 375
 376    -------------------
 377    -- Scaled_Divide --
 378    -------------------
 379
 380    procedure Scaled_Divide
 381      (X, Y, Z : Int64;
 382       Q, R    : out Int64;
 383       Round   : Boolean)
 384    is
 385       Xu  : constant Uns64 := To_Uns (abs X);
 386       Xhi : constant Uns32 := Hi (Xu);
 387       Xlo : constant Uns32 := Lo (Xu);
 388
 389       Yu  : constant Uns64 := To_Uns (abs Y);
 390       Yhi : constant Uns32 := Hi (Yu);
 391       Ylo : constant Uns32 := Lo (Yu);
 392
 393       Zu  : Uns64 := To_Uns (abs Z);
 394       Zhi : Uns32 := Hi (Zu);
 395       Zlo : Uns32 := Lo (Zu);
 396
 397       D : array (1 .. 4) of Uns32;
 398       --  The dividend, four digits (D(1) is high order)
 399
 400       Qd : array (1 .. 2) of Uns32;
 401       --  The quotient digits, two digits (Qd(1) is high order)
 402
 403       S1, S2, S3 : Uns32;
 404       --  Value to subtract, three digits (S1 is high order)
 405
 406       Qu : Uns64;
 407       Ru : Uns64;
 408       --  Unsigned quotient and remainder
 409
 410       Scale : Natural;
 411       --  Scaling factor used for multiple-precision divide. Dividend and
 412       --  Divisor are multiplied by 2 ** Scale, and the final remainder
 413       --  is divided by the scaling factor. The reason for this scaling
 414       --  is to allow more accurate estimation of quotient digits.
 415
 416       T1, T2, T3 : Uns64;
 417       --  Temporary values
 418
 419    begin
 420       --  First do the multiplication, giving the four digit dividend
 421
 422       T1 := Xlo * Ylo;
 423       D (4) := Lo (T1);
 424       D (3) := Hi (T1);
 425
 426       if Yhi /= 0 then
 427          T1 := Xlo * Yhi;
 428          T2 := D (3) + Lo (T1);
 429          D (3) := Lo (T2);
 430          D (2) := Hi (T1) + Hi (T2);
 431
 432          if Xhi /= 0 then
 433             T1 := Xhi * Ylo;
 434             T2 := D (3) + Lo (T1);
 435             D (3) := Lo (T2);
 436             T3 := D (2) + Hi (T1);
 437             T3 := T3 + Hi (T2);
 438             D (2) := Lo (T3);
 439             D (1) := Hi (T3);
 440
 441             T1 := (D (1) & D (2)) + Uns64'(Xhi * Yhi);
 442             D (1) := Hi (T1);
 443             D (2) := Lo (T1);
 444
 445          else
 446             D (1) := 0;
 447          end if;
 448
 449       else
 450          if Xhi /= 0 then
 451             T1 := Xhi * Ylo;
 452             T2 := D (3) + Lo (T1);
 453             D (3) := Lo (T2);
 454             D (2) := Hi (T1) + Hi (T2);
 455
 456          else
 457             D (2) := 0;
 458          end if;
 459
 460          D (1) := 0;
 461       end if;
 462
 463       --  Now it is time for the dreaded multiple precision division. First
 464       --  an easy case, check for the simple case of a one digit divisor.
 465
 466       if Zhi = 0 then
 467          if D (1) /= 0 or else D (2) >= Zlo then
 468             Raise_Error;
 469
 470          --  Here we are dividing at most three digits by one digit
 471
 472          else
 473             T1 := D (2) & D (3);
 474             T2 := Lo (T1 rem Zlo) & D (4);
 475
 476             Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo);
 477             Ru := T2 rem Zlo;
 478          end if;
 479
 480       --  If divisor is double digit and too large, raise error
 481
 482       elsif (D (1) & D (2)) >= Zu then
 483          Raise_Error;
 484
 485       --  This is the complex case where we definitely have a double digit
 486       --  divisor and a dividend of at least three digits. We use the classical
 487       --  multiple division algorithm (see section (4.3.1) of Knuth's "The Art
 488       --  of Computer Programming", Vol. 2 for a description (algorithm D).
 489
 490       else
 491          --  First normalize the divisor so that it has the leading bit on.
 492          --  We do this by finding the appropriate left shift amount.
 493
 494          Scale := 0;
 495
 496          if (Zhi and 16#FFFF0000#) = 0 then
 497             Scale := 16;
 498             Zu := Shift_Left (Zu, 16);
 499          end if;
 500
 501          if (Hi (Zu) and 16#FF00_0000#) = 0 then
 502             Scale := Scale + 8;
 503             Zu := Shift_Left (Zu, 8);
 504          end if;
 505
 506          if (Hi (Zu) and 16#F000_0000#) = 0 then
 507             Scale := Scale + 4;
 508             Zu := Shift_Left (Zu, 4);
 509          end if;
 510
 511          if (Hi (Zu) and 16#C000_0000#) = 0 then
 512             Scale := Scale + 2;
 513             Zu := Shift_Left (Zu, 2);
 514          end if;
 515
 516          if (Hi (Zu) and 16#8000_0000#) = 0 then
 517             Scale := Scale + 1;
 518             Zu := Shift_Left (Zu, 1);
 519          end if;
 520
 521          Zhi := Hi (Zu);
 522          Zlo := Lo (Zu);
 523
 524          --  Note that when we scale up the dividend, it still fits in four
 525          --  digits, since we already tested for overflow, and scaling does
 526          --  not change the invariant that (D (1) & D (2)) >= Zu.
 527
 528          T1 := Shift_Left (D (1) & D (2), Scale);
 529          D (1) := Hi (T1);
 530          T2 := Shift_Left (0 & D (3), Scale);
 531          D (2) := Lo (T1) or Hi (T2);
 532          T3 := Shift_Left (0 & D (4), Scale);
 533          D (3) := Lo (T2) or Hi (T3);
 534          D (4) := Lo (T3);
 535
 536          --  Loop to compute quotient digits, runs twice for Qd(1) and Qd(2).
 537
 538          for J in 0 .. 1 loop
 539
 540             --  Compute next quotient digit. We have to divide three digits by
 541             --  two digits. We estimate the quotient by dividing the leading
 542             --  two digits by the leading digit. Given the scaling we did above
 543             --  which ensured the first bit of the divisor is set, this gives
 544             --  an estimate of the quotient that is at most two too high.
 545
 546             if D (J + 1) = Zhi then
 547                Qd (J + 1) := 2 ** 32 - 1;
 548             else
 549                Qd (J + 1) := Lo ((D (J + 1) & D (J + 2)) / Zhi);
 550             end if;
 551
 552             --  Compute amount to subtract
 553
 554             T1 := Qd (J + 1) * Zlo;
 555             T2 := Qd (J + 1) * Zhi;
 556             S3 := Lo (T1);
 557             T1 := Hi (T1) + Lo (T2);
 558             S2 := Lo (T1);
 559             S1 := Hi (T1) + Hi (T2);
 560
 561             --  Adjust quotient digit if it was too high
 562
 563             loop
 564                exit when Le3 (S1, S2, S3, D (J + 1), D (J + 2), D (J + 3));
 565                Qd (J + 1) := Qd (J + 1) - 1;
 566                Sub3 (S1, S2, S3, 0, Zhi, Zlo);
 567             end loop;
 568
 569             --  Now subtract S1&S2&S3 from D1&D2&D3 ready for next step
 570
 571             Sub3 (D (J + 1), D (J + 2), D (J + 3), S1, S2, S3);
 572          end loop;
 573
 574          --  The two quotient digits are now set, and the remainder of the
 575          --  scaled division is in D3&D4. To get the remainder for the
 576          --  original unscaled division, we rescale this dividend.
 577
 578          --  We rescale the divisor as well, to make the proper comparison
 579          --  for rounding below.
 580
 581          Qu := Qd (1) & Qd (2);
 582          Ru := Shift_Right (D (3) & D (4), Scale);
 583          Zu := Shift_Right (Zu, Scale);
 584       end if;
 585
 586       --  Deal with rounding case
 587
 588       if Round and then Ru > (Zu - Uns64'(1)) / Uns64'(2) then
 589          Qu := Qu + Uns64 (1);
 590       end if;
 591
 592       --  Set final signs (RM 4.5.5(27-30))
 593
 594       --  Case of dividend (X * Y) sign positive
 595
 596       if (X >= 0 and then Y >= 0)
 597         or else (X < 0 and then Y < 0)
 598       then
 599          R := To_Pos_Int (Ru);
 600
 601          if Z > 0 then
 602             Q := To_Pos_Int (Qu);
 603          else
 604             Q := To_Neg_Int (Qu);
 605          end if;
 606
 607       --  Case of dividend (X * Y) sign negative
 608
 609       else
 610          R := To_Neg_Int (Ru);
 611
 612          if Z > 0 then
 613             Q := To_Neg_Int (Qu);
 614          else
 615             Q := To_Pos_Int (Qu);
 616          end if;
 617       end if;
 618    end Scaled_Divide;
 619
 620    ----------
 621    -- Sub3 --
 622    ----------
 623
 624    procedure Sub3 (X1, X2, X3 : in out Uns32; Y1, Y2, Y3 : in Uns32) is
 625    begin
 626       if Y3 > X3 then
 627          if X2 = 0 then
 628             X1 := X1 - 1;
 629          end if;
 630
 631          X2 := X2 - 1;
 632       end if;
 633
 634       X3 := X3 - Y3;
 635
 636       if Y2 > X2 then
 637          X1 := X1 - 1;
 638       end if;
 639
 640       X2 := X2 - Y2;
 641       X1 := X1 - Y1;
 642    end Sub3;
 643
 644    -------------------------------
 645    -- Subtract_With_Ovflo_Check --
 646    -------------------------------
 647
 648    function Subtract_With_Ovflo_Check (X, Y : Int64) return Int64 is
 649       R : constant Int64 := To_Int (To_Uns (X) - To_Uns (Y));
 650
 651    begin
 652       if X >= 0 then
 653          if Y > 0 or else R >= 0 then
 654             return R;
 655          end if;
 656
 657       else -- X < 0
 658          if Y <= 0 or else R < 0 then
 659             return R;
 660          end if;
 661       end if;
 662
 663       Raise_Error;
 664    end Subtract_With_Ovflo_Check;
 665
 666    ----------------
 667    -- To_Neg_Int --
 668    ----------------
 669
 670    function To_Neg_Int (A : Uns64) return Int64 is
 671       R : constant Int64 := -To_Int (A);
 672
 673    begin
 674       if R <= 0 then
 675          return R;
 676       else
 677          Raise_Error;
 678       end if;
 679    end To_Neg_Int;
 680
 681    ----------------
 682    -- To_Pos_Int --
 683    ----------------
 684
 685    function To_Pos_Int (A : Uns64) return Int64 is
 686       R : constant Int64 := To_Int (A);
 687
 688    begin
 689       if R >= 0 then
 690          return R;
 691       else
 692          Raise_Error;
 693       end if;
 694    end To_Pos_Int;
 695
 696 end System.Arith_64;