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