projects / gcc.git / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
history | raw | HEAD
[Ada] Fix rounding of fixed-point arithmetic operation
[gcc.git] / gcc / ada / libgnat / 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-2019, 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 -- Set final signs (RM 4.5.5(27-30))
151
152 Den_Pos := (Y < 0) = (Z < 0);
153
154 -- Compute Y * Z. Note that if the result overflows 64 bits unsigned,
155 -- then the rounded result is zero, except for the very special case
156 -- where X = -2**63 and abs(Y*Z) = 2**64, when Round is True.
157
158 if Yhi /= 0 then
159 if Zhi /= 0 then
160
161 -- Handle the special case when Round is True
162
163 if Yhi = 1
164 and then Zhi = 1
165 and then Ylo = 0
166 and then Zlo = 0
167 and then X = Int64'First
168 and then Round
169 then
170 Q := (if Den_Pos then -1 else 1);
171 else
172 Q := 0;
173 end if;
174
175 R := X;
176 return;
177 else
178 T2 := Yhi * Zlo;
179 end if;
180
181 else
182 T2 := Ylo * Zhi;
183 end if;
184
185 T1 := Ylo * Zlo;
186 T2 := T2 + Hi (T1);
187
188 if Hi (T2) /= 0 then
189
190 -- Handle the special case when Round is True
191
192 if Hi (T2) = 1
193 and then Lo (T2) = 0
194 and then Lo (T1) = 0
195 and then X = Int64'First
196 and then Round
197 then
198 Q := (if Den_Pos then -1 else 1);
199 else
200 Q := 0;
201 end if;
202
203 R := X;
204 return;
205 end if;
206
207 Du := Lo (T2) & Lo (T1);
208
209 -- Check overflow case of largest negative number divided by -1
210
211 if X = Int64'First and then Du = 1 and then not Den_Pos then
212 Raise_Error;
213 end if;
214
215 -- Perform the actual division
216
217 Qu := Xu / Du;
218 Ru := Xu rem Du;
219
220 -- Deal with rounding case
221
222 if Round and then Ru > (Du - Uns64'(1)) / Uns64'(2) then
223 Qu := Qu + Uns64'(1);
224 end if;
225
226 -- Case of dividend (X) sign positive
227
228 if X >= 0 then
229 R := To_Int (Ru);
230 Q := (if Den_Pos then To_Int (Qu) else -To_Int (Qu));
231
232 -- Case of dividend (X) sign negative
233
234 -- We perform the unary minus operation on the unsigned value
235 -- before conversion to signed, to avoid a possible overflow for
236 -- value -2**63, both for computing R and Q.
237
238 else
239 R := To_Int (-Ru);
240 Q := (if Den_Pos then To_Int (-Qu) else To_Int (Qu));
241 end if;
242 end Double_Divide;
243
244 ---------
245 -- Le3 --
246 ---------
247
248 function Le3 (X1, X2, X3 : Uns32; Y1, Y2, Y3 : Uns32) return Boolean is
249 begin
250 if X1 < Y1 then
251 return True;
252 elsif X1 > Y1 then
253 return False;
254 elsif X2 < Y2 then
255 return True;
256 elsif X2 > Y2 then
257 return False;
258 else
259 return X3 <= Y3;
260 end if;
261 end Le3;
262
263 -------------------------------
264 -- Multiply_With_Ovflo_Check --
265 -------------------------------
266
267 function Multiply_With_Ovflo_Check (X, Y : Int64) return Int64 is
268 Xu : constant Uns64 := abs X;
269 Xhi : constant Uns32 := Hi (Xu);
270 Xlo : constant Uns32 := Lo (Xu);
271
272 Yu : constant Uns64 := abs Y;
273 Yhi : constant Uns32 := Hi (Yu);
274 Ylo : constant Uns32 := Lo (Yu);
275
276 T1, T2 : Uns64;
277
278 begin
279 if Xhi /= 0 then
280 if Yhi /= 0 then
281 Raise_Error;
282 else
283 T2 := Xhi * Ylo;
284 end if;
285
286 elsif Yhi /= 0 then
287 T2 := Xlo * Yhi;
288
289 else -- Yhi = Xhi = 0
290 T2 := 0;
291 end if;
292
293 -- Here we have T2 set to the contribution to the upper half of the
294 -- result from the upper halves of the input values.
295
296 T1 := Xlo * Ylo;
297 T2 := T2 + Hi (T1);
298
299 if Hi (T2) /= 0 then
300 Raise_Error;
301 end if;
302
303 T2 := Lo (T2) & Lo (T1);
304
305 if X >= 0 then
306 if Y >= 0 then
307 return To_Pos_Int (T2);
308 else
309 return To_Neg_Int (T2);
310 end if;
311 else -- X < 0
312 if Y < 0 then
313 return To_Pos_Int (T2);
314 else
315 return To_Neg_Int (T2);
316 end if;
317 end if;
318
319 end Multiply_With_Ovflo_Check;
320
321 -----------------
322 -- Raise_Error --
323 -----------------
324
325 procedure Raise_Error is
326 begin
327 raise Constraint_Error with "64-bit arithmetic overflow";
328 end Raise_Error;
329
330 -------------------
331 -- Scaled_Divide --
332 -------------------
333
334 procedure Scaled_Divide
335 (X, Y, Z : Int64;
336 Q, R : out Int64;
337 Round : Boolean)
338 is
339 Xu : constant Uns64 := abs X;
340 Xhi : constant Uns32 := Hi (Xu);
341 Xlo : constant Uns32 := Lo (Xu);
342
343 Yu : constant Uns64 := abs Y;
344 Yhi : constant Uns32 := Hi (Yu);
345 Ylo : constant Uns32 := Lo (Yu);
346
347 Zu : Uns64 := abs Z;
348 Zhi : Uns32 := Hi (Zu);
349 Zlo : Uns32 := Lo (Zu);
350
351 D : array (1 .. 4) of Uns32;
352 -- The dividend, four digits (D(1) is high order)
353
354 Qd : array (1 .. 2) of Uns32;
355 -- The quotient digits, two digits (Qd(1) is high order)
356
357 S1, S2, S3 : Uns32;
358 -- Value to subtract, three digits (S1 is high order)
359
360 Qu : Uns64;
361 Ru : Uns64;
362 -- Unsigned quotient and remainder
363
364 Scale : Natural;
365 -- Scaling factor used for multiple-precision divide. Dividend and
366 -- Divisor are multiplied by 2 ** Scale, and the final remainder is
367 -- divided by the scaling factor. The reason for this scaling is to
368 -- allow more accurate estimation of quotient digits.
369
370 T1, T2, T3 : Uns64;
371 -- Temporary values
372
373 begin
374 -- First do the multiplication, giving the four digit dividend
375
376 T1 := Xlo * Ylo;
377 D (4) := Lo (T1);
378 D (3) := Hi (T1);
379
380 if Yhi /= 0 then
381 T1 := Xlo * Yhi;
382 T2 := D (3) + Lo (T1);
383 D (3) := Lo (T2);
384 D (2) := Hi (T1) + Hi (T2);
385
386 if Xhi /= 0 then
387 T1 := Xhi * Ylo;
388 T2 := D (3) + Lo (T1);
389 D (3) := Lo (T2);
390 T3 := D (2) + Hi (T1);
391 T3 := T3 + Hi (T2);
392 D (2) := Lo (T3);
393 D (1) := Hi (T3);
394
395 T1 := (D (1) & D (2)) + Uns64'(Xhi * Yhi);
396 D (1) := Hi (T1);
397 D (2) := Lo (T1);
398
399 else
400 D (1) := 0;
401 end if;
402
403 else
404 if Xhi /= 0 then
405 T1 := Xhi * Ylo;
406 T2 := D (3) + Lo (T1);
407 D (3) := Lo (T2);
408 D (2) := Hi (T1) + Hi (T2);
409
410 else
411 D (2) := 0;
412 end if;
413
414 D (1) := 0;
415 end if;
416
417 -- Now it is time for the dreaded multiple precision division. First an
418 -- easy case, check for the simple case of a one digit divisor.
419
420 if Zhi = 0 then
421 if D (1) /= 0 or else D (2) >= Zlo then
422 Raise_Error;
423
424 -- Here we are dividing at most three digits by one digit
425
426 else
427 T1 := D (2) & D (3);
428 T2 := Lo (T1 rem Zlo) & D (4);
429
430 Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo);
431 Ru := T2 rem Zlo;
432 end if;
433
434 -- If divisor is double digit and dividend is too large, raise error
435
436 elsif (D (1) & D (2)) >= Zu then
437 Raise_Error;
438
439 -- This is the complex case where we definitely have a double digit
440 -- divisor and a dividend of at least three digits. We use the classical
441 -- multiple-precision division algorithm (see section (4.3.1) of Knuth's
442 -- "The Art of Computer Programming", Vol. 2 for a description
443 -- (algorithm D).
444
445 else
446 -- First normalize the divisor so that it has the leading bit on.
447 -- We do this by finding the appropriate left shift amount.
448
449 Scale := 0;
450
451 if (Zhi and 16#FFFF0000#) = 0 then
452 Scale := 16;
453 Zu := Shift_Left (Zu, 16);
454 end if;
455
456 if (Hi (Zu) and 16#FF00_0000#) = 0 then
457 Scale := Scale + 8;
458 Zu := Shift_Left (Zu, 8);
459 end if;
460
461 if (Hi (Zu) and 16#F000_0000#) = 0 then
462 Scale := Scale + 4;
463 Zu := Shift_Left (Zu, 4);
464 end if;
465
466 if (Hi (Zu) and 16#C000_0000#) = 0 then
467 Scale := Scale + 2;
468 Zu := Shift_Left (Zu, 2);
469 end if;
470
471 if (Hi (Zu) and 16#8000_0000#) = 0 then
472 Scale := Scale + 1;
473 Zu := Shift_Left (Zu, 1);
474 end if;
475
476 Zhi := Hi (Zu);
477 Zlo := Lo (Zu);
478
479 -- Note that when we scale up the dividend, it still fits in four
480 -- digits, since we already tested for overflow, and scaling does
481 -- not change the invariant that (D (1) & D (2)) < Zu.
482
483 T1 := Shift_Left (D (1) & D (2), Scale);
484 D (1) := Hi (T1);
485 T2 := Shift_Left (0 & D (3), Scale);
486 D (2) := Lo (T1) or Hi (T2);
487 T3 := Shift_Left (0 & D (4), Scale);
488 D (3) := Lo (T2) or Hi (T3);
489 D (4) := Lo (T3);
490
491 -- Loop to compute quotient digits, runs twice for Qd(1) and Qd(2)
492
493 for J in 0 .. 1 loop
494
495 -- Compute next quotient digit. We have to divide three digits by
496 -- two digits. We estimate the quotient by dividing the leading
497 -- two digits by the leading digit. Given the scaling we did above
498 -- which ensured the first bit of the divisor is set, this gives
499 -- an estimate of the quotient that is at most two too high.
500
501 Qd (J + 1) := (if D (J + 1) = Zhi
502 then 2 ** 32 - 1
503 else Lo ((D (J + 1) & D (J + 2)) / Zhi));
504
505 -- Compute amount to subtract
506
507 T1 := Qd (J + 1) * Zlo;
508 T2 := Qd (J + 1) * Zhi;
509 S3 := Lo (T1);
510 T1 := Hi (T1) + Lo (T2);
511 S2 := Lo (T1);
512 S1 := Hi (T1) + Hi (T2);
513
514 -- Adjust quotient digit if it was too high
515
516 -- We use the version of the algorithm in the 2nd Edition of
517 -- "The Art of Computer Programming". This had a bug not
518 -- discovered till 1995, see Vol 2 errata:
519 -- http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz.
520 -- Under rare circumstances the expression in the test could
521 -- overflow. This version was further corrected in 2005, see
522 -- Vol 2 errata:
523 -- http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz.
524 -- This implementation is not impacted by these bugs, due to the
525 -- use of a word-size comparison done in function Le3 instead of
526 -- a comparison on two-word integer quantities in the original
527 -- algorithm.
528
529 loop
530 exit when Le3 (S1, S2, S3, D (J + 1), D (J + 2), D (J + 3));
531 Qd (J + 1) := Qd (J + 1) - 1;
532 Sub3 (S1, S2, S3, 0, Zhi, Zlo);
533 end loop;
534
535 -- Now subtract S1&S2&S3 from D1&D2&D3 ready for next step
536
537 Sub3 (D (J + 1), D (J + 2), D (J + 3), S1, S2, S3);
538 end loop;
539
540 -- The two quotient digits are now set, and the remainder of the
541 -- scaled division is in D3&D4. To get the remainder for the
542 -- original unscaled division, we rescale this dividend.
543
544 -- We rescale the divisor as well, to make the proper comparison
545 -- for rounding below.
546
547 Qu := Qd (1) & Qd (2);
548 Ru := Shift_Right (D (3) & D (4), Scale);
549 Zu := Shift_Right (Zu, Scale);
550 end if;
551
552 -- Deal with rounding case
553
554 if Round and then Ru > (Zu - Uns64'(1)) / Uns64'(2) then
555
556 -- Protect against wrapping around when rounding, by signaling
557 -- an overflow when the quotient is too large.
558
559 if Qu = Uns64'Last then
560 Raise_Error;
561 end if;
562
563 Qu := Qu + Uns64 (1);
564 end if;
565
566 -- Set final signs (RM 4.5.5(27-30))
567
568 -- Case of dividend (X * Y) sign positive
569
570 if (X >= 0 and then Y >= 0) or else (X < 0 and then Y < 0) then
571 R := To_Pos_Int (Ru);
572 Q := (if Z > 0 then To_Pos_Int (Qu) else To_Neg_Int (Qu));
573
574 -- Case of dividend (X * Y) sign negative
575
576 else
577 R := To_Neg_Int (Ru);
578 Q := (if Z > 0 then To_Neg_Int (Qu) else To_Pos_Int (Qu));
579 end if;
580 end Scaled_Divide;
581
582 ----------
583 -- Sub3 --
584 ----------
585
586 procedure Sub3 (X1, X2, X3 : in out Uns32; Y1, Y2, Y3 : Uns32) is
587 begin
588 if Y3 > X3 then
589 if X2 = 0 then
590 X1 := X1 - 1;
591 end if;
592
593 X2 := X2 - 1;
594 end if;
595
596 X3 := X3 - Y3;
597
598 if Y2 > X2 then
599 X1 := X1 - 1;
600 end if;
601
602 X2 := X2 - Y2;
603 X1 := X1 - Y1;
604 end Sub3;
605
606 -------------------------------
607 -- Subtract_With_Ovflo_Check --
608 -------------------------------
609
610 function Subtract_With_Ovflo_Check (X, Y : Int64) return Int64 is
611 R : constant Int64 := To_Int (To_Uns (X) - To_Uns (Y));
612
613 begin
614 if X >= 0 then
615 if Y > 0 or else R >= 0 then
616 return R;
617 end if;
618
619 else -- X < 0
620 if Y <= 0 or else R < 0 then
621 return R;
622 end if;
623 end if;
624
625 Raise_Error;
626 end Subtract_With_Ovflo_Check;
627
628 ----------------
629 -- To_Neg_Int --
630 ----------------
631
632 function To_Neg_Int (A : Uns64) return Int64 is
633 R : constant Int64 := (if A = 2**63 then Int64'First else -To_Int (A));
634 -- Note that we can't just use the expression of the Else, because it
635 -- overflows for A = 2**63.
636 begin
637 if R <= 0 then
638 return R;
639 else
640 Raise_Error;
641 end if;
642 end To_Neg_Int;
643
644 ----------------
645 -- To_Pos_Int --
646 ----------------
647
648 function To_Pos_Int (A : Uns64) return Int64 is
649 R : constant Int64 := To_Int (A);
650 begin
651 if R >= 0 then
652 return R;
653 else
654 Raise_Error;
655 end if;
656 end To_Pos_Int;
657
658 end System.Arith_64;