|
111
|
1 ------------------------------------------------------------------------------
|
|
|
2 -- --
|
|
|
3 -- GNAT COMPILER COMPONENTS --
|
|
|
4 -- --
|
|
|
5 -- S Y S T E M . F A T _ G E N --
|
|
|
6 -- --
|
|
|
7 -- B o d y --
|
|
|
8 -- --
|
|
|
9 -- Copyright (C) 1992-2017, 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 -- The implementation here is portable to any IEEE implementation. It does
|
|
|
33 -- not handle nonbinary radix, and also assumes that model numbers and
|
|
|
34 -- machine numbers are basically identical, which is not true of all possible
|
|
|
35 -- floating-point implementations. On a non-IEEE machine, this body must be
|
|
|
36 -- specialized appropriately, or better still, its generic instantiations
|
|
|
37 -- should be replaced by efficient machine-specific code.
|
|
|
38
|
|
|
39 with Ada.Unchecked_Conversion;
|
|
|
40 with System;
|
|
|
41 package body System.Fat_Gen is
|
|
|
42
|
|
|
43 Float_Radix : constant T := T (T'Machine_Radix);
|
|
|
44 Radix_To_M_Minus_1 : constant T := Float_Radix ** (T'Machine_Mantissa - 1);
|
|
|
45
|
|
|
46 pragma Assert (T'Machine_Radix = 2);
|
|
|
47 -- This version does not handle radix 16
|
|
|
48
|
|
|
49 -- Constants for Decompose and Scaling
|
|
|
50
|
|
|
51 Rad : constant T := T (T'Machine_Radix);
|
|
|
52 Invrad : constant T := 1.0 / Rad;
|
|
|
53
|
|
|
54 subtype Expbits is Integer range 0 .. 6;
|
|
|
55 -- 2 ** (2 ** 7) might overflow. How big can radix-16 exponents get?
|
|
|
56
|
|
|
57 Log_Power : constant array (Expbits) of Integer := (1, 2, 4, 8, 16, 32, 64);
|
|
|
58
|
|
|
59 R_Power : constant array (Expbits) of T :=
|
|
|
60 (Rad ** 1,
|
|
|
61 Rad ** 2,
|
|
|
62 Rad ** 4,
|
|
|
63 Rad ** 8,
|
|
|
64 Rad ** 16,
|
|
|
65 Rad ** 32,
|
|
|
66 Rad ** 64);
|
|
|
67
|
|
|
68 R_Neg_Power : constant array (Expbits) of T :=
|
|
|
69 (Invrad ** 1,
|
|
|
70 Invrad ** 2,
|
|
|
71 Invrad ** 4,
|
|
|
72 Invrad ** 8,
|
|
|
73 Invrad ** 16,
|
|
|
74 Invrad ** 32,
|
|
|
75 Invrad ** 64);
|
|
|
76
|
|
|
77 -----------------------
|
|
|
78 -- Local Subprograms --
|
|
|
79 -----------------------
|
|
|
80
|
|
|
81 procedure Decompose (XX : T; Frac : out T; Expo : out UI);
|
|
|
82 -- Decomposes a floating-point number into fraction and exponent parts.
|
|
|
83 -- Both results are signed, with Frac having the sign of XX, and UI has
|
|
|
84 -- the sign of the exponent. The absolute value of Frac is in the range
|
|
|
85 -- 0.0 <= Frac < 1.0. If Frac = 0.0 or -0.0, then Expo is always zero.
|
|
|
86
|
|
|
87 function Gradual_Scaling (Adjustment : UI) return T;
|
|
|
88 -- Like Scaling with a first argument of 1.0, but returns the smallest
|
|
|
89 -- denormal rather than zero when the adjustment is smaller than
|
|
|
90 -- Machine_Emin. Used for Succ and Pred.
|
|
|
91
|
|
|
92 --------------
|
|
|
93 -- Adjacent --
|
|
|
94 --------------
|
|
|
95
|
|
|
96 function Adjacent (X, Towards : T) return T is
|
|
|
97 begin
|
|
|
98 if Towards = X then
|
|
|
99 return X;
|
|
|
100 elsif Towards > X then
|
|
|
101 return Succ (X);
|
|
|
102 else
|
|
|
103 return Pred (X);
|
|
|
104 end if;
|
|
|
105 end Adjacent;
|
|
|
106
|
|
|
107 -------------
|
|
|
108 -- Ceiling --
|
|
|
109 -------------
|
|
|
110
|
|
|
111 function Ceiling (X : T) return T is
|
|
|
112 XT : constant T := Truncation (X);
|
|
|
113 begin
|
|
|
114 if X <= 0.0 then
|
|
|
115 return XT;
|
|
|
116 elsif X = XT then
|
|
|
117 return X;
|
|
|
118 else
|
|
|
119 return XT + 1.0;
|
|
|
120 end if;
|
|
|
121 end Ceiling;
|
|
|
122
|
|
|
123 -------------
|
|
|
124 -- Compose --
|
|
|
125 -------------
|
|
|
126
|
|
|
127 function Compose (Fraction : T; Exponent : UI) return T is
|
|
|
128 Arg_Frac : T;
|
|
|
129 Arg_Exp : UI;
|
|
|
130 pragma Unreferenced (Arg_Exp);
|
|
|
131 begin
|
|
|
132 Decompose (Fraction, Arg_Frac, Arg_Exp);
|
|
|
133 return Scaling (Arg_Frac, Exponent);
|
|
|
134 end Compose;
|
|
|
135
|
|
|
136 ---------------
|
|
|
137 -- Copy_Sign --
|
|
|
138 ---------------
|
|
|
139
|
|
|
140 function Copy_Sign (Value, Sign : T) return T is
|
|
|
141 Result : T;
|
|
|
142
|
|
|
143 function Is_Negative (V : T) return Boolean;
|
|
|
144 pragma Import (Intrinsic, Is_Negative);
|
|
|
145
|
|
|
146 begin
|
|
|
147 Result := abs Value;
|
|
|
148
|
|
|
149 if Is_Negative (Sign) then
|
|
|
150 return -Result;
|
|
|
151 else
|
|
|
152 return Result;
|
|
|
153 end if;
|
|
|
154 end Copy_Sign;
|
|
|
155
|
|
|
156 ---------------
|
|
|
157 -- Decompose --
|
|
|
158 ---------------
|
|
|
159
|
|
|
160 procedure Decompose (XX : T; Frac : out T; Expo : out UI) is
|
|
|
161 X : constant T := T'Machine (XX);
|
|
|
162
|
|
|
163 begin
|
|
|
164 if X = 0.0 then
|
|
|
165
|
|
|
166 -- The normalized exponent of zero is zero, see RM A.5.2(15)
|
|
|
167
|
|
|
168 Frac := X;
|
|
|
169 Expo := 0;
|
|
|
170
|
|
|
171 -- Check for infinities, transfinites, whatnot
|
|
|
172
|
|
|
173 elsif X > T'Safe_Last then
|
|
|
174 Frac := Invrad;
|
|
|
175 Expo := T'Machine_Emax + 1;
|
|
|
176
|
|
|
177 elsif X < T'Safe_First then
|
|
|
178 Frac := -Invrad;
|
|
|
179 Expo := T'Machine_Emax + 2; -- how many extra negative values?
|
|
|
180
|
|
|
181 else
|
|
|
182 -- Case of nonzero finite x. Essentially, we just multiply
|
|
|
183 -- by Rad ** (+-2**N) to reduce the range.
|
|
|
184
|
|
|
185 declare
|
|
|
186 Ax : T := abs X;
|
|
|
187 Ex : UI := 0;
|
|
|
188
|
|
|
189 -- Ax * Rad ** Ex is invariant
|
|
|
190
|
|
|
191 begin
|
|
|
192 if Ax >= 1.0 then
|
|
|
193 while Ax >= R_Power (Expbits'Last) loop
|
|
|
194 Ax := Ax * R_Neg_Power (Expbits'Last);
|
|
|
195 Ex := Ex + Log_Power (Expbits'Last);
|
|
|
196 end loop;
|
|
|
197
|
|
|
198 -- Ax < Rad ** 64
|
|
|
199
|
|
|
200 for N in reverse Expbits'First .. Expbits'Last - 1 loop
|
|
|
201 if Ax >= R_Power (N) then
|
|
|
202 Ax := Ax * R_Neg_Power (N);
|
|
|
203 Ex := Ex + Log_Power (N);
|
|
|
204 end if;
|
|
|
205
|
|
|
206 -- Ax < R_Power (N)
|
|
|
207
|
|
|
208 end loop;
|
|
|
209
|
|
|
210 -- 1 <= Ax < Rad
|
|
|
211
|
|
|
212 Ax := Ax * Invrad;
|
|
|
213 Ex := Ex + 1;
|
|
|
214
|
|
|
215 else
|
|
|
216 -- 0 < ax < 1
|
|
|
217
|
|
|
218 while Ax < R_Neg_Power (Expbits'Last) loop
|
|
|
219 Ax := Ax * R_Power (Expbits'Last);
|
|
|
220 Ex := Ex - Log_Power (Expbits'Last);
|
|
|
221 end loop;
|
|
|
222
|
|
|
223 -- Rad ** -64 <= Ax < 1
|
|
|
224
|
|
|
225 for N in reverse Expbits'First .. Expbits'Last - 1 loop
|
|
|
226 if Ax < R_Neg_Power (N) then
|
|
|
227 Ax := Ax * R_Power (N);
|
|
|
228 Ex := Ex - Log_Power (N);
|
|
|
229 end if;
|
|
|
230
|
|
|
231 -- R_Neg_Power (N) <= Ax < 1
|
|
|
232
|
|
|
233 end loop;
|
|
|
234 end if;
|
|
|
235
|
|
|
236 Frac := (if X > 0.0 then Ax else -Ax);
|
|
|
237 Expo := Ex;
|
|
|
238 end;
|
|
|
239 end if;
|
|
|
240 end Decompose;
|
|
|
241
|
|
|
242 --------------
|
|
|
243 -- Exponent --
|
|
|
244 --------------
|
|
|
245
|
|
|
246 function Exponent (X : T) return UI is
|
|
|
247 X_Frac : T;
|
|
|
248 X_Exp : UI;
|
|
|
249 pragma Unreferenced (X_Frac);
|
|
|
250 begin
|
|
|
251 Decompose (X, X_Frac, X_Exp);
|
|
|
252 return X_Exp;
|
|
|
253 end Exponent;
|
|
|
254
|
|
|
255 -----------
|
|
|
256 -- Floor --
|
|
|
257 -----------
|
|
|
258
|
|
|
259 function Floor (X : T) return T is
|
|
|
260 XT : constant T := Truncation (X);
|
|
|
261 begin
|
|
|
262 if X >= 0.0 then
|
|
|
263 return XT;
|
|
|
264 elsif XT = X then
|
|
|
265 return X;
|
|
|
266 else
|
|
|
267 return XT - 1.0;
|
|
|
268 end if;
|
|
|
269 end Floor;
|
|
|
270
|
|
|
271 --------------
|
|
|
272 -- Fraction --
|
|
|
273 --------------
|
|
|
274
|
|
|
275 function Fraction (X : T) return T is
|
|
|
276 X_Frac : T;
|
|
|
277 X_Exp : UI;
|
|
|
278 pragma Unreferenced (X_Exp);
|
|
|
279 begin
|
|
|
280 Decompose (X, X_Frac, X_Exp);
|
|
|
281 return X_Frac;
|
|
|
282 end Fraction;
|
|
|
283
|
|
|
284 ---------------------
|
|
|
285 -- Gradual_Scaling --
|
|
|
286 ---------------------
|
|
|
287
|
|
|
288 function Gradual_Scaling (Adjustment : UI) return T is
|
|
|
289 Y : T;
|
|
|
290 Y1 : T;
|
|
|
291 Ex : UI := Adjustment;
|
|
|
292
|
|
|
293 begin
|
|
|
294 if Adjustment < T'Machine_Emin - 1 then
|
|
|
295 Y := 2.0 ** T'Machine_Emin;
|
|
|
296 Y1 := Y;
|
|
|
297 Ex := Ex - T'Machine_Emin;
|
|
|
298 while Ex < 0 loop
|
|
|
299 Y := T'Machine (Y / 2.0);
|
|
|
300
|
|
|
301 if Y = 0.0 then
|
|
|
302 return Y1;
|
|
|
303 end if;
|
|
|
304
|
|
|
305 Ex := Ex + 1;
|
|
|
306 Y1 := Y;
|
|
|
307 end loop;
|
|
|
308
|
|
|
309 return Y1;
|
|
|
310
|
|
|
311 else
|
|
|
312 return Scaling (1.0, Adjustment);
|
|
|
313 end if;
|
|
|
314 end Gradual_Scaling;
|
|
|
315
|
|
|
316 ------------------
|
|
|
317 -- Leading_Part --
|
|
|
318 ------------------
|
|
|
319
|
|
|
320 function Leading_Part (X : T; Radix_Digits : UI) return T is
|
|
|
321 L : UI;
|
|
|
322 Y, Z : T;
|
|
|
323
|
|
|
324 begin
|
|
|
325 if Radix_Digits >= T'Machine_Mantissa then
|
|
|
326 return X;
|
|
|
327
|
|
|
328 elsif Radix_Digits <= 0 then
|
|
|
329 raise Constraint_Error;
|
|
|
330
|
|
|
331 else
|
|
|
332 L := Exponent (X) - Radix_Digits;
|
|
|
333 Y := Truncation (Scaling (X, -L));
|
|
|
334 Z := Scaling (Y, L);
|
|
|
335 return Z;
|
|
|
336 end if;
|
|
|
337 end Leading_Part;
|
|
|
338
|
|
|
339 -------------
|
|
|
340 -- Machine --
|
|
|
341 -------------
|
|
|
342
|
|
|
343 -- The trick with Machine is to force the compiler to store the result
|
|
|
344 -- in memory so that we do not have extra precision used. The compiler
|
|
|
345 -- is clever, so we have to outwit its possible optimizations. We do
|
|
|
346 -- this by using an intermediate pragma Volatile location.
|
|
|
347
|
|
|
348 function Machine (X : T) return T is
|
|
|
349 Temp : T;
|
|
|
350 pragma Volatile (Temp);
|
|
|
351 begin
|
|
|
352 Temp := X;
|
|
|
353 return Temp;
|
|
|
354 end Machine;
|
|
|
355
|
|
|
356 ----------------------
|
|
|
357 -- Machine_Rounding --
|
|
|
358 ----------------------
|
|
|
359
|
|
|
360 -- For now, the implementation is identical to that of Rounding, which is
|
|
|
361 -- a permissible behavior, but is not the most efficient possible approach.
|
|
|
362
|
|
|
363 function Machine_Rounding (X : T) return T is
|
|
|
364 Result : T;
|
|
|
365 Tail : T;
|
|
|
366
|
|
|
367 begin
|
|
|
368 Result := Truncation (abs X);
|
|
|
369 Tail := abs X - Result;
|
|
|
370
|
|
|
371 if Tail >= 0.5 then
|
|
|
372 Result := Result + 1.0;
|
|
|
373 end if;
|
|
|
374
|
|
|
375 if X > 0.0 then
|
|
|
376 return Result;
|
|
|
377
|
|
|
378 elsif X < 0.0 then
|
|
|
379 return -Result;
|
|
|
380
|
|
|
381 -- For zero case, make sure sign of zero is preserved
|
|
|
382
|
|
|
383 else
|
|
|
384 return X;
|
|
|
385 end if;
|
|
|
386 end Machine_Rounding;
|
|
|
387
|
|
|
388 -----------
|
|
|
389 -- Model --
|
|
|
390 -----------
|
|
|
391
|
|
|
392 -- We treat Model as identical to Machine. This is true of IEEE and other
|
|
|
393 -- nice floating-point systems, but not necessarily true of all systems.
|
|
|
394
|
|
|
395 function Model (X : T) return T is
|
|
|
396 begin
|
|
|
397 return Machine (X);
|
|
|
398 end Model;
|
|
|
399
|
|
|
400 ----------
|
|
|
401 -- Pred --
|
|
|
402 ----------
|
|
|
403
|
|
|
404 function Pred (X : T) return T is
|
|
|
405 X_Frac : T;
|
|
|
406 X_Exp : UI;
|
|
|
407
|
|
|
408 begin
|
|
|
409 -- Zero has to be treated specially, since its exponent is zero
|
|
|
410
|
|
|
411 if X = 0.0 then
|
|
|
412 return -Succ (X);
|
|
|
413
|
|
|
414 -- Special treatment for most negative number
|
|
|
415
|
|
|
416 elsif X = T'First then
|
|
|
417
|
|
|
418 -- If not generating infinities, we raise a constraint error
|
|
|
419
|
|
|
420 if T'Machine_Overflows then
|
|
|
421 raise Constraint_Error with "Pred of largest negative number";
|
|
|
422
|
|
|
423 -- Otherwise generate a negative infinity
|
|
|
424
|
|
|
425 else
|
|
|
426 return X / (X - X);
|
|
|
427 end if;
|
|
|
428
|
|
|
429 -- For infinities, return unchanged
|
|
|
430
|
|
|
431 elsif X < T'First or else X > T'Last then
|
|
|
432 return X;
|
|
|
433
|
|
|
434 -- Subtract from the given number a number equivalent to the value
|
|
|
435 -- of its least significant bit. Given that the most significant bit
|
|
|
436 -- represents a value of 1.0 * radix ** (exp - 1), the value we want
|
|
|
437 -- is obtained by shifting this by (mantissa-1) bits to the right,
|
|
|
438 -- i.e. decreasing the exponent by that amount.
|
|
|
439
|
|
|
440 else
|
|
|
441 Decompose (X, X_Frac, X_Exp);
|
|
|
442
|
|
|
443 -- A special case, if the number we had was a positive power of
|
|
|
444 -- two, then we want to subtract half of what we would otherwise
|
|
|
445 -- subtract, since the exponent is going to be reduced.
|
|
|
446
|
|
|
447 -- Note that X_Frac has the same sign as X, so if X_Frac is 0.5,
|
|
|
448 -- then we know that we have a positive number (and hence a
|
|
|
449 -- positive power of 2).
|
|
|
450
|
|
|
451 if X_Frac = 0.5 then
|
|
|
452 return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
|
|
|
453
|
|
|
454 -- Otherwise the exponent is unchanged
|
|
|
455
|
|
|
456 else
|
|
|
457 return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa);
|
|
|
458 end if;
|
|
|
459 end if;
|
|
|
460 end Pred;
|
|
|
461
|
|
|
462 ---------------
|
|
|
463 -- Remainder --
|
|
|
464 ---------------
|
|
|
465
|
|
|
466 function Remainder (X, Y : T) return T is
|
|
|
467 A : T;
|
|
|
468 B : T;
|
|
|
469 Arg : T;
|
|
|
470 P : T;
|
|
|
471 P_Frac : T;
|
|
|
472 Sign_X : T;
|
|
|
473 IEEE_Rem : T;
|
|
|
474 Arg_Exp : UI;
|
|
|
475 P_Exp : UI;
|
|
|
476 K : UI;
|
|
|
477 P_Even : Boolean;
|
|
|
478
|
|
|
479 Arg_Frac : T;
|
|
|
480 pragma Unreferenced (Arg_Frac);
|
|
|
481
|
|
|
482 begin
|
|
|
483 if Y = 0.0 then
|
|
|
484 raise Constraint_Error;
|
|
|
485 end if;
|
|
|
486
|
|
|
487 if X > 0.0 then
|
|
|
488 Sign_X := 1.0;
|
|
|
489 Arg := X;
|
|
|
490 else
|
|
|
491 Sign_X := -1.0;
|
|
|
492 Arg := -X;
|
|
|
493 end if;
|
|
|
494
|
|
|
495 P := abs Y;
|
|
|
496
|
|
|
497 if Arg < P then
|
|
|
498 P_Even := True;
|
|
|
499 IEEE_Rem := Arg;
|
|
|
500 P_Exp := Exponent (P);
|
|
|
501
|
|
|
502 else
|
|
|
503 Decompose (Arg, Arg_Frac, Arg_Exp);
|
|
|
504 Decompose (P, P_Frac, P_Exp);
|
|
|
505
|
|
|
506 P := Compose (P_Frac, Arg_Exp);
|
|
|
507 K := Arg_Exp - P_Exp;
|
|
|
508 P_Even := True;
|
|
|
509 IEEE_Rem := Arg;
|
|
|
510
|
|
|
511 for Cnt in reverse 0 .. K loop
|
|
|
512 if IEEE_Rem >= P then
|
|
|
513 P_Even := False;
|
|
|
514 IEEE_Rem := IEEE_Rem - P;
|
|
|
515 else
|
|
|
516 P_Even := True;
|
|
|
517 end if;
|
|
|
518
|
|
|
519 P := P * 0.5;
|
|
|
520 end loop;
|
|
|
521 end if;
|
|
|
522
|
|
|
523 -- That completes the calculation of modulus remainder. The final
|
|
|
524 -- step is get the IEEE remainder. Here we need to compare Rem with
|
|
|
525 -- (abs Y) / 2. We must be careful of unrepresentable Y/2 value
|
|
|
526 -- caused by subnormal numbers
|
|
|
527
|
|
|
528 if P_Exp >= 0 then
|
|
|
529 A := IEEE_Rem;
|
|
|
530 B := abs Y * 0.5;
|
|
|
531
|
|
|
532 else
|
|
|
533 A := IEEE_Rem * 2.0;
|
|
|
534 B := abs Y;
|
|
|
535 end if;
|
|
|
536
|
|
|
537 if A > B or else (A = B and then not P_Even) then
|
|
|
538 IEEE_Rem := IEEE_Rem - abs Y;
|
|
|
539 end if;
|
|
|
540
|
|
|
541 return Sign_X * IEEE_Rem;
|
|
|
542 end Remainder;
|
|
|
543
|
|
|
544 --------------
|
|
|
545 -- Rounding --
|
|
|
546 --------------
|
|
|
547
|
|
|
548 function Rounding (X : T) return T is
|
|
|
549 Result : T;
|
|
|
550 Tail : T;
|
|
|
551
|
|
|
552 begin
|
|
|
553 Result := Truncation (abs X);
|
|
|
554 Tail := abs X - Result;
|
|
|
555
|
|
|
556 if Tail >= 0.5 then
|
|
|
557 Result := Result + 1.0;
|
|
|
558 end if;
|
|
|
559
|
|
|
560 if X > 0.0 then
|
|
|
561 return Result;
|
|
|
562
|
|
|
563 elsif X < 0.0 then
|
|
|
564 return -Result;
|
|
|
565
|
|
|
566 -- For zero case, make sure sign of zero is preserved
|
|
|
567
|
|
|
568 else
|
|
|
569 return X;
|
|
|
570 end if;
|
|
|
571 end Rounding;
|
|
|
572
|
|
|
573 -------------
|
|
|
574 -- Scaling --
|
|
|
575 -------------
|
|
|
576
|
|
|
577 -- Return x * rad ** adjustment quickly, or quietly underflow to zero,
|
|
|
578 -- or overflow naturally.
|
|
|
579
|
|
|
580 function Scaling (X : T; Adjustment : UI) return T is
|
|
|
581 begin
|
|
|
582 if X = 0.0 or else Adjustment = 0 then
|
|
|
583 return X;
|
|
|
584 end if;
|
|
|
585
|
|
|
586 -- Nonzero x essentially, just multiply repeatedly by Rad ** (+-2**n)
|
|
|
587
|
|
|
588 declare
|
|
|
589 Y : T := X;
|
|
|
590 Ex : UI := Adjustment;
|
|
|
591
|
|
|
592 -- Y * Rad ** Ex is invariant
|
|
|
593
|
|
|
594 begin
|
|
|
595 if Ex < 0 then
|
|
|
596 while Ex <= -Log_Power (Expbits'Last) loop
|
|
|
597 Y := Y * R_Neg_Power (Expbits'Last);
|
|
|
598 Ex := Ex + Log_Power (Expbits'Last);
|
|
|
599 end loop;
|
|
|
600
|
|
|
601 -- -64 < Ex <= 0
|
|
|
602
|
|
|
603 for N in reverse Expbits'First .. Expbits'Last - 1 loop
|
|
|
604 if Ex <= -Log_Power (N) then
|
|
|
605 Y := Y * R_Neg_Power (N);
|
|
|
606 Ex := Ex + Log_Power (N);
|
|
|
607 end if;
|
|
|
608
|
|
|
609 -- -Log_Power (N) < Ex <= 0
|
|
|
610
|
|
|
611 end loop;
|
|
|
612
|
|
|
613 -- Ex = 0
|
|
|
614
|
|
|
615 else
|
|
|
616 -- Ex >= 0
|
|
|
617
|
|
|
618 while Ex >= Log_Power (Expbits'Last) loop
|
|
|
619 Y := Y * R_Power (Expbits'Last);
|
|
|
620 Ex := Ex - Log_Power (Expbits'Last);
|
|
|
621 end loop;
|
|
|
622
|
|
|
623 -- 0 <= Ex < 64
|
|
|
624
|
|
|
625 for N in reverse Expbits'First .. Expbits'Last - 1 loop
|
|
|
626 if Ex >= Log_Power (N) then
|
|
|
627 Y := Y * R_Power (N);
|
|
|
628 Ex := Ex - Log_Power (N);
|
|
|
629 end if;
|
|
|
630
|
|
|
631 -- 0 <= Ex < Log_Power (N)
|
|
|
632
|
|
|
633 end loop;
|
|
|
634
|
|
|
635 -- Ex = 0
|
|
|
636
|
|
|
637 end if;
|
|
|
638
|
|
|
639 return Y;
|
|
|
640 end;
|
|
|
641 end Scaling;
|
|
|
642
|
|
|
643 ----------
|
|
|
644 -- Succ --
|
|
|
645 ----------
|
|
|
646
|
|
|
647 function Succ (X : T) return T is
|
|
|
648 X_Frac : T;
|
|
|
649 X_Exp : UI;
|
|
|
650 X1, X2 : T;
|
|
|
651
|
|
|
652 begin
|
|
|
653 -- Treat zero specially since it has a zero exponent
|
|
|
654
|
|
|
655 if X = 0.0 then
|
|
|
656 X1 := 2.0 ** T'Machine_Emin;
|
|
|
657
|
|
|
658 -- Following loop generates smallest denormal
|
|
|
659
|
|
|
660 loop
|
|
|
661 X2 := T'Machine (X1 / 2.0);
|
|
|
662 exit when X2 = 0.0;
|
|
|
663 X1 := X2;
|
|
|
664 end loop;
|
|
|
665
|
|
|
666 return X1;
|
|
|
667
|
|
|
668 -- Special treatment for largest positive number
|
|
|
669
|
|
|
670 elsif X = T'Last then
|
|
|
671
|
|
|
672 -- If not generating infinities, we raise a constraint error
|
|
|
673
|
|
|
674 if T'Machine_Overflows then
|
|
|
675 raise Constraint_Error with "Succ of largest negative number";
|
|
|
676
|
|
|
677 -- Otherwise generate a positive infinity
|
|
|
678
|
|
|
679 else
|
|
|
680 return X / (X - X);
|
|
|
681 end if;
|
|
|
682
|
|
|
683 -- For infinities, return unchanged
|
|
|
684
|
|
|
685 elsif X < T'First or else X > T'Last then
|
|
|
686 return X;
|
|
|
687
|
|
|
688 -- Add to the given number a number equivalent to the value
|
|
|
689 -- of its least significant bit. Given that the most significant bit
|
|
|
690 -- represents a value of 1.0 * radix ** (exp - 1), the value we want
|
|
|
691 -- is obtained by shifting this by (mantissa-1) bits to the right,
|
|
|
692 -- i.e. decreasing the exponent by that amount.
|
|
|
693
|
|
|
694 else
|
|
|
695 Decompose (X, X_Frac, X_Exp);
|
|
|
696
|
|
|
697 -- A special case, if the number we had was a negative power of two,
|
|
|
698 -- then we want to add half of what we would otherwise add, since the
|
|
|
699 -- exponent is going to be reduced.
|
|
|
700
|
|
|
701 -- Note that X_Frac has the same sign as X, so if X_Frac is -0.5,
|
|
|
702 -- then we know that we have a negative number (and hence a negative
|
|
|
703 -- power of 2).
|
|
|
704
|
|
|
705 if X_Frac = -0.5 then
|
|
|
706 return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1);
|
|
|
707
|
|
|
708 -- Otherwise the exponent is unchanged
|
|
|
709
|
|
|
710 else
|
|
|
711 return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa);
|
|
|
712 end if;
|
|
|
713 end if;
|
|
|
714 end Succ;
|
|
|
715
|
|
|
716 ----------------
|
|
|
717 -- Truncation --
|
|
|
718 ----------------
|
|
|
719
|
|
|
720 -- The basic approach is to compute
|
|
|
721
|
|
|
722 -- T'Machine (RM1 + N) - RM1
|
|
|
723
|
|
|
724 -- where N >= 0.0 and RM1 = radix ** (mantissa - 1)
|
|
|
725
|
|
|
726 -- This works provided that the intermediate result (RM1 + N) does not
|
|
|
727 -- have extra precision (which is why we call Machine). When we compute
|
|
|
728 -- RM1 + N, the exponent of N will be normalized and the mantissa shifted
|
|
|
729 -- appropriately so the lower order bits, which cannot contribute to the
|
|
|
730 -- integer part of N, fall off on the right. When we subtract RM1 again,
|
|
|
731 -- the significant bits of N are shifted to the left, and what we have is
|
|
|
732 -- an integer, because only the first e bits are different from zero
|
|
|
733 -- (assuming binary radix here).
|
|
|
734
|
|
|
735 function Truncation (X : T) return T is
|
|
|
736 Result : T;
|
|
|
737
|
|
|
738 begin
|
|
|
739 Result := abs X;
|
|
|
740
|
|
|
741 if Result >= Radix_To_M_Minus_1 then
|
|
|
742 return Machine (X);
|
|
|
743
|
|
|
744 else
|
|
|
745 Result := Machine (Radix_To_M_Minus_1 + Result) - Radix_To_M_Minus_1;
|
|
|
746
|
|
|
747 if Result > abs X then
|
|
|
748 Result := Result - 1.0;
|
|
|
749 end if;
|
|
|
750
|
|
|
751 if X > 0.0 then
|
|
|
752 return Result;
|
|
|
753
|
|
|
754 elsif X < 0.0 then
|
|
|
755 return -Result;
|
|
|
756
|
|
|
757 -- For zero case, make sure sign of zero is preserved
|
|
|
758
|
|
|
759 else
|
|
|
760 return X;
|
|
|
761 end if;
|
|
|
762 end if;
|
|
|
763 end Truncation;
|
|
|
764
|
|
|
765 -----------------------
|
|
|
766 -- Unbiased_Rounding --
|
|
|
767 -----------------------
|
|
|
768
|
|
|
769 function Unbiased_Rounding (X : T) return T is
|
|
|
770 Abs_X : constant T := abs X;
|
|
|
771 Result : T;
|
|
|
772 Tail : T;
|
|
|
773
|
|
|
774 begin
|
|
|
775 Result := Truncation (Abs_X);
|
|
|
776 Tail := Abs_X - Result;
|
|
|
777
|
|
|
778 if Tail > 0.5 then
|
|
|
779 Result := Result + 1.0;
|
|
|
780
|
|
|
781 elsif Tail = 0.5 then
|
|
|
782 Result := 2.0 * Truncation ((Result / 2.0) + 0.5);
|
|
|
783 end if;
|
|
|
784
|
|
|
785 if X > 0.0 then
|
|
|
786 return Result;
|
|
|
787
|
|
|
788 elsif X < 0.0 then
|
|
|
789 return -Result;
|
|
|
790
|
|
|
791 -- For zero case, make sure sign of zero is preserved
|
|
|
792
|
|
|
793 else
|
|
|
794 return X;
|
|
|
795 end if;
|
|
|
796 end Unbiased_Rounding;
|
|
|
797
|
|
|
798 -----------
|
|
|
799 -- Valid --
|
|
|
800 -----------
|
|
|
801
|
|
|
802 function Valid (X : not null access T) return Boolean is
|
|
|
803 IEEE_Emin : constant Integer := T'Machine_Emin - 1;
|
|
|
804 IEEE_Emax : constant Integer := T'Machine_Emax - 1;
|
|
|
805
|
|
|
806 IEEE_Bias : constant Integer := -(IEEE_Emin - 1);
|
|
|
807
|
|
|
808 subtype IEEE_Exponent_Range is
|
|
|
809 Integer range IEEE_Emin - 1 .. IEEE_Emax + 1;
|
|
|
810
|
|
|
811 -- The implementation of this floating point attribute uses a
|
|
|
812 -- representation type Float_Rep that allows direct access to the
|
|
|
813 -- exponent and mantissa parts of a floating point number.
|
|
|
814
|
|
|
815 -- The Float_Rep type is an array of Float_Word elements. This
|
|
|
816 -- representation is chosen to make it possible to size the type based
|
|
|
817 -- on a generic parameter. Since the array size is known at compile
|
|
|
818 -- time, efficient code can still be generated. The size of Float_Word
|
|
|
819 -- elements should be large enough to allow accessing the exponent in
|
|
|
820 -- one read, but small enough so that all floating point object sizes
|
|
|
821 -- are a multiple of the Float_Word'Size.
|
|
|
822
|
|
|
823 -- The following conditions must be met for all possible instantiations
|
|
|
824 -- of the attributes package:
|
|
|
825
|
|
|
826 -- - T'Size is an integral multiple of Float_Word'Size
|
|
|
827
|
|
|
828 -- - The exponent and sign are completely contained in a single
|
|
|
829 -- component of Float_Rep, named Most_Significant_Word (MSW).
|
|
|
830
|
|
|
831 -- - The sign occupies the most significant bit of the MSW and the
|
|
|
832 -- exponent is in the following bits. Unused bits (if any) are in
|
|
|
833 -- the least significant part.
|
|
|
834
|
|
|
835 type Float_Word is mod 2**Positive'Min (System.Word_Size, 32);
|
|
|
836 type Rep_Index is range 0 .. 7;
|
|
|
837
|
|
|
838 Rep_Words : constant Positive :=
|
|
|
839 (T'Size + Float_Word'Size - 1) / Float_Word'Size;
|
|
|
840 Rep_Last : constant Rep_Index :=
|
|
|
841 Rep_Index'Min
|
|
|
842 (Rep_Index (Rep_Words - 1),
|
|
|
843 (T'Mantissa + 16) / Float_Word'Size);
|
|
|
844 -- Determine the number of Float_Words needed for representing the
|
|
|
845 -- entire floating-point value. Do not take into account excessive
|
|
|
846 -- padding, as occurs on IA-64 where 80 bits floats get padded to 128
|
|
|
847 -- bits. In general, the exponent field cannot be larger than 15 bits,
|
|
|
848 -- even for 128-bit floating-point types, so the final format size
|
|
|
849 -- won't be larger than T'Mantissa + 16.
|
|
|
850
|
|
|
851 type Float_Rep is
|
|
|
852 array (Rep_Index range 0 .. Rep_Index (Rep_Words - 1)) of Float_Word;
|
|
|
853
|
|
|
854 pragma Suppress_Initialization (Float_Rep);
|
|
|
855 -- This pragma suppresses the generation of an initialization procedure
|
|
|
856 -- for type Float_Rep when operating in Initialize/Normalize_Scalars
|
|
|
857 -- mode. This is not just a matter of efficiency, but of functionality,
|
|
|
858 -- since Valid has a pragma Inline_Always, which is not permitted if
|
|
|
859 -- there are nested subprograms present.
|
|
|
860
|
|
|
861 Most_Significant_Word : constant Rep_Index :=
|
|
|
862 Rep_Last * Standard'Default_Bit_Order;
|
|
|
863 -- Finding the location of the Exponent_Word is a bit tricky. In general
|
|
|
864 -- we assume Word_Order = Bit_Order.
|
|
|
865
|
|
|
866 Exponent_Factor : constant Float_Word :=
|
|
|
867 2**(Float_Word'Size - 1) /
|
|
|
868 Float_Word (IEEE_Emax - IEEE_Emin + 3) *
|
|
|
869 Boolean'Pos (Most_Significant_Word /= 2) +
|
|
|
870 Boolean'Pos (Most_Significant_Word = 2);
|
|
|
871 -- Factor that the extracted exponent needs to be divided by to be in
|
|
|
872 -- range 0 .. IEEE_Emax - IEEE_Emin + 2. Special case: Exponent_Factor
|
|
|
873 -- is 1 for x86/IA64 double extended (GCC adds unused bits to the type).
|
|
|
874
|
|
|
875 Exponent_Mask : constant Float_Word :=
|
|
|
876 Float_Word (IEEE_Emax - IEEE_Emin + 2) *
|
|
|
877 Exponent_Factor;
|
|
|
878 -- Value needed to mask out the exponent field. This assumes that the
|
|
|
879 -- range IEEE_Emin - 1 .. IEEE_Emax + contains 2**N values, for some N
|
|
|
880 -- in Natural.
|
|
|
881
|
|
|
882 function To_Float is new Ada.Unchecked_Conversion (Float_Rep, T);
|
|
|
883
|
|
|
884 type Float_Access is access all T;
|
|
|
885 function To_Address is
|
|
|
886 new Ada.Unchecked_Conversion (Float_Access, System.Address);
|
|
|
887
|
|
|
888 XA : constant System.Address := To_Address (Float_Access (X));
|
|
|
889
|
|
|
890 R : Float_Rep;
|
|
|
891 pragma Import (Ada, R);
|
|
|
892 for R'Address use XA;
|
|
|
893 -- R is a view of the input floating-point parameter. Note that we
|
|
|
894 -- must avoid copying the actual bits of this parameter in float
|
|
|
895 -- form (since it may be a signalling NaN).
|
|
|
896
|
|
|
897 E : constant IEEE_Exponent_Range :=
|
|
|
898 Integer ((R (Most_Significant_Word) and Exponent_Mask) /
|
|
|
899 Exponent_Factor)
|
|
|
900 - IEEE_Bias;
|
|
|
901 -- Mask/Shift T to only get bits from the exponent. Then convert biased
|
|
|
902 -- value to integer value.
|
|
|
903
|
|
|
904 SR : Float_Rep;
|
|
|
905 -- Float_Rep representation of significant of X.all
|
|
|
906
|
|
|
907 begin
|
|
|
908 if T'Denorm then
|
|
|
909
|
|
|
910 -- All denormalized numbers are valid, so the only invalid numbers
|
|
|
911 -- are overflows and NaNs, both with exponent = Emax + 1.
|
|
|
912
|
|
|
913 return E /= IEEE_Emax + 1;
|
|
|
914
|
|
|
915 end if;
|
|
|
916
|
|
|
917 -- All denormalized numbers except 0.0 are invalid
|
|
|
918
|
|
|
919 -- Set exponent of X to zero, so we end up with the significand, which
|
|
|
920 -- definitely is a valid number and can be converted back to a float.
|
|
|
921
|
|
|
922 SR := R;
|
|
|
923 SR (Most_Significant_Word) :=
|
|
|
924 (SR (Most_Significant_Word)
|
|
|
925 and not Exponent_Mask) + Float_Word (IEEE_Bias) * Exponent_Factor;
|
|
|
926
|
|
|
927 return (E in IEEE_Emin .. IEEE_Emax) or else
|
|
|
928 ((E = IEEE_Emin - 1) and then abs To_Float (SR) = 1.0);
|
|
|
929 end Valid;
|
|
|
930
|
|
|
931 end System.Fat_Gen;
|
