File : eval_fat.adb
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT COMPILER COMPONENTS --
4 -- --
5 -- E V A L _ F A T --
6 -- --
7 -- B o d y --
8 -- --
9 -- Copyright (C) 1992-2016, 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. 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 COPYING3. If not, go to --
19 -- http://www.gnu.org/licenses for a complete copy of the license. --
20 -- --
21 -- GNAT was originally developed by the GNAT team at New York University. --
22 -- Extensive contributions were provided by Ada Core Technologies Inc. --
23 -- --
24 ------------------------------------------------------------------------------
25
26 with Einfo; use Einfo;
27 with Errout; use Errout;
28 with Sem_Util; use Sem_Util;
29
30 package body Eval_Fat is
31
32 Radix : constant Int := 2;
33 -- This code is currently only correct for the radix 2 case. We use the
34 -- symbolic value Radix where possible to help in the unlikely case of
35 -- anyone ever having to adjust this code for another value, and for
36 -- documentation purposes.
37
38 -- Another assumption is that the range of the floating-point type is
39 -- symmetric around zero.
40
41 type Radix_Power_Table is array (Int range 1 .. 4) of Int;
42
43 Radix_Powers : constant Radix_Power_Table :=
44 (Radix ** 1, Radix ** 2, Radix ** 3, Radix ** 4);
45
46 -----------------------
47 -- Local Subprograms --
48 -----------------------
49
50 procedure Decompose
51 (RT : R;
52 X : T;
53 Fraction : out T;
54 Exponent : out UI;
55 Mode : Rounding_Mode := Round);
56 -- Decomposes a non-zero floating-point number into fraction and exponent
57 -- parts. The fraction is in the interval 1.0 / Radix .. T'Pred (1.0) and
58 -- uses Rbase = Radix. The result is rounded to a nearest machine number.
59
60 --------------
61 -- Adjacent --
62 --------------
63
64 function Adjacent (RT : R; X, Towards : T) return T is
65 begin
66 if Towards = X then
67 return X;
68 elsif Towards > X then
69 return Succ (RT, X);
70 else
71 return Pred (RT, X);
72 end if;
73 end Adjacent;
74
75 -------------
76 -- Ceiling --
77 -------------
78
79 function Ceiling (RT : R; X : T) return T is
80 XT : constant T := Truncation (RT, X);
81 begin
82 if UR_Is_Negative (X) then
83 return XT;
84 elsif X = XT then
85 return X;
86 else
87 return XT + Ureal_1;
88 end if;
89 end Ceiling;
90
91 -------------
92 -- Compose --
93 -------------
94
95 function Compose (RT : R; Fraction : T; Exponent : UI) return T is
96 Arg_Frac : T;
97 Arg_Exp : UI;
98 pragma Warnings (Off, Arg_Exp);
99 begin
100 Decompose (RT, Fraction, Arg_Frac, Arg_Exp);
101 return Scaling (RT, Arg_Frac, Exponent);
102 end Compose;
103
104 ---------------
105 -- Copy_Sign --
106 ---------------
107
108 function Copy_Sign (RT : R; Value, Sign : T) return T is
109 pragma Warnings (Off, RT);
110 Result : T;
111
112 begin
113 Result := abs Value;
114
115 if UR_Is_Negative (Sign) then
116 return -Result;
117 else
118 return Result;
119 end if;
120 end Copy_Sign;
121
122 ---------------
123 -- Decompose --
124 ---------------
125
126 procedure Decompose
127 (RT : R;
128 X : T;
129 Fraction : out T;
130 Exponent : out UI;
131 Mode : Rounding_Mode := Round)
132 is
133 Int_F : UI;
134
135 begin
136 Decompose_Int (RT, abs X, Int_F, Exponent, Mode);
137
138 Fraction := UR_From_Components
139 (Num => Int_F,
140 Den => Machine_Mantissa_Value (RT),
141 Rbase => Radix,
142 Negative => False);
143
144 if UR_Is_Negative (X) then
145 Fraction := -Fraction;
146 end if;
147
148 return;
149 end Decompose;
150
151 -------------------
152 -- Decompose_Int --
153 -------------------
154
155 -- This procedure should be modified with care, as there are many non-
156 -- obvious details that may cause problems that are hard to detect. For
157 -- zero arguments, Fraction and Exponent are set to zero. Note that sign
158 -- of zero cannot be preserved.
159
160 procedure Decompose_Int
161 (RT : R;
162 X : T;
163 Fraction : out UI;
164 Exponent : out UI;
165 Mode : Rounding_Mode)
166 is
167 Base : Int := Rbase (X);
168 N : UI := abs Numerator (X);
169 D : UI := Denominator (X);
170
171 N_Times_Radix : UI;
172
173 Even : Boolean;
174 -- True iff Fraction is even
175
176 Most_Significant_Digit : constant UI :=
177 Radix ** (Machine_Mantissa_Value (RT) - 1);
178
179 Uintp_Mark : Uintp.Save_Mark;
180 -- The code is divided into blocks that systematically release
181 -- intermediate values (this routine generates lots of junk).
182
183 begin
184 if N = Uint_0 then
185 Fraction := Uint_0;
186 Exponent := Uint_0;
187 return;
188 end if;
189
190 Calculate_D_And_Exponent_1 : begin
191 Uintp_Mark := Mark;
192 Exponent := Uint_0;
193
194 -- In cases where Base > 1, the actual denominator is Base**D. For
195 -- cases where Base is a power of Radix, use the value 1 for the
196 -- Denominator and adjust the exponent.
197
198 -- Note: Exponent has different sign from D, because D is a divisor
199
200 for Power in 1 .. Radix_Powers'Last loop
201 if Base = Radix_Powers (Power) then
202 Exponent := -D * Power;
203 Base := 0;
204 D := Uint_1;
205 exit;
206 end if;
207 end loop;
208
209 Release_And_Save (Uintp_Mark, D, Exponent);
210 end Calculate_D_And_Exponent_1;
211
212 if Base > 0 then
213 Calculate_Exponent : begin
214 Uintp_Mark := Mark;
215
216 -- For bases that are a multiple of the Radix, divide the base by
217 -- Radix and adjust the Exponent. This will help because D will be
218 -- much smaller and faster to process.
219
220 -- This occurs for decimal bases on machines with binary floating-
221 -- point for example. When calculating 1E40, with Radix = 2, N
222 -- will be 93 bits instead of 133.
223
224 -- N E
225 -- ------ * Radix
226 -- D
227 -- Base
228
229 -- N E
230 -- = -------------------------- * Radix
231 -- D D
232 -- (Base/Radix) * Radix
233
234 -- N E-D
235 -- = --------------- * Radix
236 -- D
237 -- (Base/Radix)
238
239 -- This code is commented out, because it causes numerous
240 -- failures in the regression suite. To be studied ???
241
242 while False and then Base > 0 and then Base mod Radix = 0 loop
243 Base := Base / Radix;
244 Exponent := Exponent + D;
245 end loop;
246
247 Release_And_Save (Uintp_Mark, Exponent);
248 end Calculate_Exponent;
249
250 -- For remaining bases we must actually compute the exponentiation
251
252 -- Because the exponentiation can be negative, and D must be integer,
253 -- the numerator is corrected instead.
254
255 Calculate_N_And_D : begin
256 Uintp_Mark := Mark;
257
258 if D < 0 then
259 N := N * Base ** (-D);
260 D := Uint_1;
261 else
262 D := Base ** D;
263 end if;
264
265 Release_And_Save (Uintp_Mark, N, D);
266 end Calculate_N_And_D;
267
268 Base := 0;
269 end if;
270
271 -- Now scale N and D so that N / D is a value in the interval [1.0 /
272 -- Radix, 1.0) and adjust Exponent accordingly, so the value N / D *
273 -- Radix ** Exponent remains unchanged.
274
275 -- Step 1 - Adjust N so N / D >= 1 / Radix, or N = 0
276
277 -- N and D are positive, so N / D >= 1 / Radix implies N * Radix >= D.
278 -- As this scaling is not possible for N is Uint_0, zero is handled
279 -- explicitly at the start of this subprogram.
280
281 Calculate_N_And_Exponent : begin
282 Uintp_Mark := Mark;
283
284 N_Times_Radix := N * Radix;
285 while not (N_Times_Radix >= D) loop
286 N := N_Times_Radix;
287 Exponent := Exponent - 1;
288 N_Times_Radix := N * Radix;
289 end loop;
290
291 Release_And_Save (Uintp_Mark, N, Exponent);
292 end Calculate_N_And_Exponent;
293
294 -- Step 2 - Adjust D so N / D < 1
295
296 -- Scale up D so N / D < 1, so N < D
297
298 Calculate_D_And_Exponent_2 : begin
299 Uintp_Mark := Mark;
300
301 while not (N < D) loop
302
303 -- As N / D >= 1, N / (D * Radix) will be at least 1 / Radix, so
304 -- the result of Step 1 stays valid
305
306 D := D * Radix;
307 Exponent := Exponent + 1;
308 end loop;
309
310 Release_And_Save (Uintp_Mark, D, Exponent);
311 end Calculate_D_And_Exponent_2;
312
313 -- Here the value N / D is in the range [1.0 / Radix .. 1.0)
314
315 -- Now find the fraction by doing a very simple-minded division until
316 -- enough digits have been computed.
317
318 -- This division works for all radices, but is only efficient for a
319 -- binary radix. It is just like a manual division algorithm, but
320 -- instead of moving the denominator one digit right, we move the
321 -- numerator one digit left so the numerator and denominator remain
322 -- integral.
323
324 Fraction := Uint_0;
325 Even := True;
326
327 Calculate_Fraction_And_N : begin
328 Uintp_Mark := Mark;
329
330 loop
331 while N >= D loop
332 N := N - D;
333 Fraction := Fraction + 1;
334 Even := not Even;
335 end loop;
336
337 -- Stop when the result is in [1.0 / Radix, 1.0)
338
339 exit when Fraction >= Most_Significant_Digit;
340
341 N := N * Radix;
342 Fraction := Fraction * Radix;
343 Even := True;
344 end loop;
345
346 Release_And_Save (Uintp_Mark, Fraction, N);
347 end Calculate_Fraction_And_N;
348
349 Calculate_Fraction_And_Exponent : begin
350 Uintp_Mark := Mark;
351
352 -- Determine correct rounding based on the remainder which is in
353 -- N and the divisor D. The rounding is performed on the absolute
354 -- value of X, so Ceiling and Floor need to check for the sign of
355 -- X explicitly.
356
357 case Mode is
358 when Round_Even =>
359
360 -- This rounding mode corresponds to the unbiased rounding
361 -- method that is used at run time. When the real value is
362 -- exactly between two machine numbers, choose the machine
363 -- number with its least significant bit equal to zero.
364
365 -- The recommendation advice in RM 4.9(38) is that static
366 -- expressions are rounded to machine numbers in the same
367 -- way as the target machine does.
368
369 if (Even and then N * 2 > D)
370 or else
371 (not Even and then N * 2 >= D)
372 then
373 Fraction := Fraction + 1;
374 end if;
375
376 when Round =>
377
378 -- Do not round to even as is done with IEEE arithmetic, but
379 -- instead round away from zero when the result is exactly
380 -- between two machine numbers. This biased rounding method
381 -- should not be used to convert static expressions to
382 -- machine numbers, see AI95-268.
383
384 if N * 2 >= D then
385 Fraction := Fraction + 1;
386 end if;
387
388 when Ceiling =>
389 if N > Uint_0 and then not UR_Is_Negative (X) then
390 Fraction := Fraction + 1;
391 end if;
392
393 when Floor =>
394 if N > Uint_0 and then UR_Is_Negative (X) then
395 Fraction := Fraction + 1;
396 end if;
397 end case;
398
399 -- The result must be normalized to [1.0/Radix, 1.0), so adjust if
400 -- the result is 1.0 because of rounding.
401
402 if Fraction = Most_Significant_Digit * Radix then
403 Fraction := Most_Significant_Digit;
404 Exponent := Exponent + 1;
405 end if;
406
407 -- Put back sign after applying the rounding
408
409 if UR_Is_Negative (X) then
410 Fraction := -Fraction;
411 end if;
412
413 Release_And_Save (Uintp_Mark, Fraction, Exponent);
414 end Calculate_Fraction_And_Exponent;
415 end Decompose_Int;
416
417 --------------
418 -- Exponent --
419 --------------
420
421 function Exponent (RT : R; X : T) return UI is
422 X_Frac : UI;
423 X_Exp : UI;
424 pragma Warnings (Off, X_Frac);
425 begin
426 Decompose_Int (RT, X, X_Frac, X_Exp, Round_Even);
427 return X_Exp;
428 end Exponent;
429
430 -----------
431 -- Floor --
432 -----------
433
434 function Floor (RT : R; X : T) return T is
435 XT : constant T := Truncation (RT, X);
436
437 begin
438 if UR_Is_Positive (X) then
439 return XT;
440
441 elsif XT = X then
442 return X;
443
444 else
445 return XT - Ureal_1;
446 end if;
447 end Floor;
448
449 --------------
450 -- Fraction --
451 --------------
452
453 function Fraction (RT : R; X : T) return T is
454 X_Frac : T;
455 X_Exp : UI;
456 pragma Warnings (Off, X_Exp);
457 begin
458 Decompose (RT, X, X_Frac, X_Exp);
459 return X_Frac;
460 end Fraction;
461
462 ------------------
463 -- Leading_Part --
464 ------------------
465
466 function Leading_Part (RT : R; X : T; Radix_Digits : UI) return T is
467 RD : constant UI := UI_Min (Radix_Digits, Machine_Mantissa_Value (RT));
468 L : UI;
469 Y : T;
470 begin
471 L := Exponent (RT, X) - RD;
472 Y := UR_From_Uint (UR_Trunc (Scaling (RT, X, -L)));
473 return Scaling (RT, Y, L);
474 end Leading_Part;
475
476 -------------
477 -- Machine --
478 -------------
479
480 function Machine
481 (RT : R;
482 X : T;
483 Mode : Rounding_Mode;
484 Enode : Node_Id) return T
485 is
486 X_Frac : T;
487 X_Exp : UI;
488 Emin : constant UI := Machine_Emin_Value (RT);
489
490 begin
491 Decompose (RT, X, X_Frac, X_Exp, Mode);
492
493 -- Case of denormalized number or (gradual) underflow
494
495 -- A denormalized number is one with the minimum exponent Emin, but that
496 -- breaks the assumption that the first digit of the mantissa is a one.
497 -- This allows the first non-zero digit to be in any of the remaining
498 -- Mant - 1 spots. The gap between subsequent denormalized numbers is
499 -- the same as for the smallest normalized numbers. However, the number
500 -- of significant digits left decreases as a result of the mantissa now
501 -- having leading seros.
502
503 if X_Exp < Emin then
504 declare
505 Emin_Den : constant UI := Machine_Emin_Value (RT)
506 - Machine_Mantissa_Value (RT) + Uint_1;
507 begin
508 if X_Exp < Emin_Den or not Has_Denormals (RT) then
509 if Has_Signed_Zeros (RT) and then UR_Is_Negative (X) then
510 Error_Msg_N
511 ("floating-point value underflows to -0.0??", Enode);
512 return Ureal_M_0;
513
514 else
515 Error_Msg_N
516 ("floating-point value underflows to 0.0??", Enode);
517 return Ureal_0;
518 end if;
519
520 elsif Has_Denormals (RT) then
521
522 -- Emin - Mant <= X_Exp < Emin, so result is denormal. Handle
523 -- gradual underflow by first computing the number of
524 -- significant bits still available for the mantissa and
525 -- then truncating the fraction to this number of bits.
526
527 -- If this value is different from the original fraction,
528 -- precision is lost due to gradual underflow.
529
530 -- We probably should round here and prevent double rounding as
531 -- a result of first rounding to a model number and then to a
532 -- machine number. However, this is an extremely rare case that
533 -- is not worth the extra complexity. In any case, a warning is
534 -- issued in cases where gradual underflow occurs.
535
536 declare
537 Denorm_Sig_Bits : constant UI := X_Exp - Emin_Den + 1;
538
539 X_Frac_Denorm : constant T := UR_From_Components
540 (UR_Trunc (Scaling (RT, abs X_Frac, Denorm_Sig_Bits)),
541 Denorm_Sig_Bits,
542 Radix,
543 UR_Is_Negative (X));
544
545 begin
546 if X_Frac_Denorm /= X_Frac then
547 Error_Msg_N
548 ("gradual underflow causes loss of precision??",
549 Enode);
550 X_Frac := X_Frac_Denorm;
551 end if;
552 end;
553 end if;
554 end;
555 end if;
556
557 return Scaling (RT, X_Frac, X_Exp);
558 end Machine;
559
560 -----------
561 -- Model --
562 -----------
563
564 function Model (RT : R; X : T) return T is
565 X_Frac : T;
566 X_Exp : UI;
567 begin
568 Decompose (RT, X, X_Frac, X_Exp);
569 return Compose (RT, X_Frac, X_Exp);
570 end Model;
571
572 ----------
573 -- Pred --
574 ----------
575
576 function Pred (RT : R; X : T) return T is
577 begin
578 return -Succ (RT, -X);
579 end Pred;
580
581 ---------------
582 -- Remainder --
583 ---------------
584
585 function Remainder (RT : R; X, Y : T) return T is
586 A : T;
587 B : T;
588 Arg : T;
589 P : T;
590 Arg_Frac : T;
591 P_Frac : T;
592 Sign_X : T;
593 IEEE_Rem : T;
594 Arg_Exp : UI;
595 P_Exp : UI;
596 K : UI;
597 P_Even : Boolean;
598
599 pragma Warnings (Off, Arg_Frac);
600
601 begin
602 if UR_Is_Positive (X) then
603 Sign_X := Ureal_1;
604 else
605 Sign_X := -Ureal_1;
606 end if;
607
608 Arg := abs X;
609 P := abs Y;
610
611 if Arg < P then
612 P_Even := True;
613 IEEE_Rem := Arg;
614 P_Exp := Exponent (RT, P);
615
616 else
617 -- ??? what about zero cases?
618 Decompose (RT, Arg, Arg_Frac, Arg_Exp);
619 Decompose (RT, P, P_Frac, P_Exp);
620
621 P := Compose (RT, P_Frac, Arg_Exp);
622 K := Arg_Exp - P_Exp;
623 P_Even := True;
624 IEEE_Rem := Arg;
625
626 for Cnt in reverse 0 .. UI_To_Int (K) loop
627 if IEEE_Rem >= P then
628 P_Even := False;
629 IEEE_Rem := IEEE_Rem - P;
630 else
631 P_Even := True;
632 end if;
633
634 P := P * Ureal_Half;
635 end loop;
636 end if;
637
638 -- That completes the calculation of modulus remainder. The final step
639 -- is get the IEEE remainder. Here we compare Rem with (abs Y) / 2.
640
641 if P_Exp >= 0 then
642 A := IEEE_Rem;
643 B := abs Y * Ureal_Half;
644
645 else
646 A := IEEE_Rem * Ureal_2;
647 B := abs Y;
648 end if;
649
650 if A > B or else (A = B and then not P_Even) then
651 IEEE_Rem := IEEE_Rem - abs Y;
652 end if;
653
654 return Sign_X * IEEE_Rem;
655 end Remainder;
656
657 --------------
658 -- Rounding --
659 --------------
660
661 function Rounding (RT : R; X : T) return T is
662 Result : T;
663 Tail : T;
664
665 begin
666 Result := Truncation (RT, abs X);
667 Tail := abs X - Result;
668
669 if Tail >= Ureal_Half then
670 Result := Result + Ureal_1;
671 end if;
672
673 if UR_Is_Negative (X) then
674 return -Result;
675 else
676 return Result;
677 end if;
678 end Rounding;
679
680 -------------
681 -- Scaling --
682 -------------
683
684 function Scaling (RT : R; X : T; Adjustment : UI) return T is
685 pragma Warnings (Off, RT);
686
687 begin
688 if Rbase (X) = Radix then
689 return UR_From_Components
690 (Num => Numerator (X),
691 Den => Denominator (X) - Adjustment,
692 Rbase => Radix,
693 Negative => UR_Is_Negative (X));
694
695 elsif Adjustment >= 0 then
696 return X * Radix ** Adjustment;
697 else
698 return X / Radix ** (-Adjustment);
699 end if;
700 end Scaling;
701
702 ----------
703 -- Succ --
704 ----------
705
706 function Succ (RT : R; X : T) return T is
707 Emin : constant UI := Machine_Emin_Value (RT);
708 Mantissa : constant UI := Machine_Mantissa_Value (RT);
709 Exp : UI := UI_Max (Emin, Exponent (RT, X));
710 Frac : T;
711 New_Frac : T;
712
713 begin
714 if UR_Is_Zero (X) then
715 Exp := Emin;
716 end if;
717
718 -- Set exponent such that the radix point will be directly following the
719 -- mantissa after scaling.
720
721 if Has_Denormals (RT) or Exp /= Emin then
722 Exp := Exp - Mantissa;
723 else
724 Exp := Exp - 1;
725 end if;
726
727 Frac := Scaling (RT, X, -Exp);
728 New_Frac := Ceiling (RT, Frac);
729
730 if New_Frac = Frac then
731 if New_Frac = Scaling (RT, -Ureal_1, Mantissa - 1) then
732 New_Frac := New_Frac + Scaling (RT, Ureal_1, Uint_Minus_1);
733 else
734 New_Frac := New_Frac + Ureal_1;
735 end if;
736 end if;
737
738 return Scaling (RT, New_Frac, Exp);
739 end Succ;
740
741 ----------------
742 -- Truncation --
743 ----------------
744
745 function Truncation (RT : R; X : T) return T is
746 pragma Warnings (Off, RT);
747 begin
748 return UR_From_Uint (UR_Trunc (X));
749 end Truncation;
750
751 -----------------------
752 -- Unbiased_Rounding --
753 -----------------------
754
755 function Unbiased_Rounding (RT : R; X : T) return T is
756 Abs_X : constant T := abs X;
757 Result : T;
758 Tail : T;
759
760 begin
761 Result := Truncation (RT, Abs_X);
762 Tail := Abs_X - Result;
763
764 if Tail > Ureal_Half then
765 Result := Result + Ureal_1;
766
767 elsif Tail = Ureal_Half then
768 Result := Ureal_2 *
769 Truncation (RT, (Result / Ureal_2) + Ureal_Half);
770 end if;
771
772 if UR_Is_Negative (X) then
773 return -Result;
774 elsif UR_Is_Positive (X) then
775 return Result;
776
777 -- For zero case, make sure sign of zero is preserved
778
779 else
780 return X;
781 end if;
782 end Unbiased_Rounding;
783
784 end Eval_Fat;