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-2018, 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 39with Ada.Unchecked_Conversion; 40with System; 41package 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 T'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 T'Machine (X); 743 744 else 745 Result := 746 T'Machine (Radix_To_M_Minus_1 + Result) - Radix_To_M_Minus_1; 747 748 if Result > abs X then 749 Result := Result - 1.0; 750 end if; 751 752 if X > 0.0 then 753 return Result; 754 755 elsif X < 0.0 then 756 return -Result; 757 758 -- For zero case, make sure sign of zero is preserved 759 760 else 761 return X; 762 end if; 763 end if; 764 end Truncation; 765 766 ----------------------- 767 -- Unbiased_Rounding -- 768 ----------------------- 769 770 function Unbiased_Rounding (X : T) return T is 771 Abs_X : constant T := abs X; 772 Result : T; 773 Tail : T; 774 775 begin 776 Result := Truncation (Abs_X); 777 Tail := Abs_X - Result; 778 779 if Tail > 0.5 then 780 Result := Result + 1.0; 781 782 elsif Tail = 0.5 then 783 Result := 2.0 * Truncation ((Result / 2.0) + 0.5); 784 end if; 785 786 if X > 0.0 then 787 return Result; 788 789 elsif X < 0.0 then 790 return -Result; 791 792 -- For zero case, make sure sign of zero is preserved 793 794 else 795 return X; 796 end if; 797 end Unbiased_Rounding; 798 799 ----------- 800 -- Valid -- 801 ----------- 802 803 function Valid (X : not null access T) return Boolean is 804 IEEE_Emin : constant Integer := T'Machine_Emin - 1; 805 IEEE_Emax : constant Integer := T'Machine_Emax - 1; 806 807 IEEE_Bias : constant Integer := -(IEEE_Emin - 1); 808 809 subtype IEEE_Exponent_Range is 810 Integer range IEEE_Emin - 1 .. IEEE_Emax + 1; 811 812 -- The implementation of this floating point attribute uses a 813 -- representation type Float_Rep that allows direct access to the 814 -- exponent and mantissa parts of a floating point number. 815 816 -- The Float_Rep type is an array of Float_Word elements. This 817 -- representation is chosen to make it possible to size the type based 818 -- on a generic parameter. Since the array size is known at compile 819 -- time, efficient code can still be generated. The size of Float_Word 820 -- elements should be large enough to allow accessing the exponent in 821 -- one read, but small enough so that all floating point object sizes 822 -- are a multiple of the Float_Word'Size. 823 824 -- The following conditions must be met for all possible instantiations 825 -- of the attributes package: 826 827 -- - T'Size is an integral multiple of Float_Word'Size 828 829 -- - The exponent and sign are completely contained in a single 830 -- component of Float_Rep, named Most_Significant_Word (MSW). 831 832 -- - The sign occupies the most significant bit of the MSW and the 833 -- exponent is in the following bits. Unused bits (if any) are in 834 -- the least significant part. 835 836 type Float_Word is mod 2**Positive'Min (System.Word_Size, 32); 837 type Rep_Index is range 0 .. 7; 838 839 Rep_Words : constant Positive := 840 (T'Size + Float_Word'Size - 1) / Float_Word'Size; 841 Rep_Last : constant Rep_Index := 842 Rep_Index'Min 843 (Rep_Index (Rep_Words - 1), 844 (T'Mantissa + 16) / Float_Word'Size); 845 -- Determine the number of Float_Words needed for representing the 846 -- entire floating-point value. Do not take into account excessive 847 -- padding, as occurs on IA-64 where 80 bits floats get padded to 128 848 -- bits. In general, the exponent field cannot be larger than 15 bits, 849 -- even for 128-bit floating-point types, so the final format size 850 -- won't be larger than T'Mantissa + 16. 851 852 type Float_Rep is 853 array (Rep_Index range 0 .. Rep_Index (Rep_Words - 1)) of Float_Word; 854 855 pragma Suppress_Initialization (Float_Rep); 856 -- This pragma suppresses the generation of an initialization procedure 857 -- for type Float_Rep when operating in Initialize/Normalize_Scalars 858 -- mode. This is not just a matter of efficiency, but of functionality, 859 -- since Valid has a pragma Inline_Always, which is not permitted if 860 -- there are nested subprograms present. 861 862 Most_Significant_Word : constant Rep_Index := 863 Rep_Last * Standard'Default_Bit_Order; 864 -- Finding the location of the Exponent_Word is a bit tricky. In general 865 -- we assume Word_Order = Bit_Order. 866 867 Exponent_Factor : constant Float_Word := 868 2**(Float_Word'Size - 1) / 869 Float_Word (IEEE_Emax - IEEE_Emin + 3) * 870 Boolean'Pos (Most_Significant_Word /= 2) + 871 Boolean'Pos (Most_Significant_Word = 2); 872 -- Factor that the extracted exponent needs to be divided by to be in 873 -- range 0 .. IEEE_Emax - IEEE_Emin + 2. Special case: Exponent_Factor 874 -- is 1 for x86/IA64 double extended (GCC adds unused bits to the type). 875 876 Exponent_Mask : constant Float_Word := 877 Float_Word (IEEE_Emax - IEEE_Emin + 2) * 878 Exponent_Factor; 879 -- Value needed to mask out the exponent field. This assumes that the 880 -- range IEEE_Emin - 1 .. IEEE_Emax + contains 2**N values, for some N 881 -- in Natural. 882 883 function To_Float is new Ada.Unchecked_Conversion (Float_Rep, T); 884 885 type Float_Access is access all T; 886 function To_Address is 887 new Ada.Unchecked_Conversion (Float_Access, System.Address); 888 889 XA : constant System.Address := To_Address (Float_Access (X)); 890 891 R : Float_Rep; 892 pragma Import (Ada, R); 893 for R'Address use XA; 894 -- R is a view of the input floating-point parameter. Note that we 895 -- must avoid copying the actual bits of this parameter in float 896 -- form (since it may be a signalling NaN). 897 898 E : constant IEEE_Exponent_Range := 899 Integer ((R (Most_Significant_Word) and Exponent_Mask) / 900 Exponent_Factor) 901 - IEEE_Bias; 902 -- Mask/Shift T to only get bits from the exponent. Then convert biased 903 -- value to integer value. 904 905 SR : Float_Rep; 906 -- Float_Rep representation of significant of X.all 907 908 begin 909 if T'Denorm then 910 911 -- All denormalized numbers are valid, so the only invalid numbers 912 -- are overflows and NaNs, both with exponent = Emax + 1. 913 914 return E /= IEEE_Emax + 1; 915 916 end if; 917 918 -- All denormalized numbers except 0.0 are invalid 919 920 -- Set exponent of X to zero, so we end up with the significand, which 921 -- definitely is a valid number and can be converted back to a float. 922 923 SR := R; 924 SR (Most_Significant_Word) := 925 (SR (Most_Significant_Word) 926 and not Exponent_Mask) + Float_Word (IEEE_Bias) * Exponent_Factor; 927 928 return (E in IEEE_Emin .. IEEE_Emax) or else 929 ((E = IEEE_Emin - 1) and then abs To_Float (SR) = 1.0); 930 end Valid; 931 932end System.Fat_Gen; 933