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-2015, 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 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 -- shifted appropriately so the lower order bits, which cannot contribute 730 -- to the integer part of N, fall off on the right. When we subtract RM1 731 -- again, the significant bits of N are shifted to the left, and what we 732 -- have is an integer, because only the first e bits are different from 733 -- zero (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 931end System.Fat_Gen; 932