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