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