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