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