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