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-2015, 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 32with Interfaces; use Interfaces; 33 34with Ada.Unchecked_Conversion; 35 36package 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 := (if A = 2**63 then Int64'First else -To_Int (A)); 581 -- Note that we can't just use the expression of the Else, because it 582 -- overflows for A = 2**63. 583 begin 584 if R <= 0 then 585 return R; 586 else 587 Raise_Error; 588 end if; 589 end To_Neg_Int; 590 591 ---------------- 592 -- To_Pos_Int -- 593 ---------------- 594 595 function To_Pos_Int (A : Uns64) return Int64 is 596 R : constant Int64 := To_Int (A); 597 begin 598 if R >= 0 then 599 return R; 600 else 601 Raise_Error; 602 end if; 603 end To_Pos_Int; 604 605end System.Arith_64; 606