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