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