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