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-2019, 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 raise Constraint_Error with "Pred of largest negative number"; 419 420 -- For infinities, return unchanged 421 422 elsif X < T'First or else X > T'Last then 423 return X; 424 425 -- Subtract from the given number a number equivalent to the value 426 -- of its least significant bit. Given that the most significant bit 427 -- represents a value of 1.0 * radix ** (exp - 1), the value we want 428 -- is obtained by shifting this by (mantissa-1) bits to the right, 429 -- i.e. decreasing the exponent by that amount. 430 431 else 432 Decompose (X, X_Frac, X_Exp); 433 434 -- A special case, if the number we had was a positive power of 435 -- two, then we want to subtract half of what we would otherwise 436 -- subtract, since the exponent is going to be reduced. 437 438 -- Note that X_Frac has the same sign as X, so if X_Frac is 0.5, 439 -- then we know that we have a positive number (and hence a 440 -- positive power of 2). 441 442 if X_Frac = 0.5 then 443 return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1); 444 445 -- Otherwise the exponent is unchanged 446 447 else 448 return X - Gradual_Scaling (X_Exp - T'Machine_Mantissa); 449 end if; 450 end if; 451 end Pred; 452 453 --------------- 454 -- Remainder -- 455 --------------- 456 457 function Remainder (X, Y : T) return T is 458 A : T; 459 B : T; 460 Arg : T; 461 P : T; 462 P_Frac : T; 463 Sign_X : T; 464 IEEE_Rem : T; 465 Arg_Exp : UI; 466 P_Exp : UI; 467 K : UI; 468 P_Even : Boolean; 469 470 Arg_Frac : T; 471 pragma Unreferenced (Arg_Frac); 472 473 begin 474 if Y = 0.0 then 475 raise Constraint_Error; 476 end if; 477 478 if X > 0.0 then 479 Sign_X := 1.0; 480 Arg := X; 481 else 482 Sign_X := -1.0; 483 Arg := -X; 484 end if; 485 486 P := abs Y; 487 488 if Arg < P then 489 P_Even := True; 490 IEEE_Rem := Arg; 491 P_Exp := Exponent (P); 492 493 else 494 Decompose (Arg, Arg_Frac, Arg_Exp); 495 Decompose (P, P_Frac, P_Exp); 496 497 P := Compose (P_Frac, Arg_Exp); 498 K := Arg_Exp - P_Exp; 499 P_Even := True; 500 IEEE_Rem := Arg; 501 502 for Cnt in reverse 0 .. K loop 503 if IEEE_Rem >= P then 504 P_Even := False; 505 IEEE_Rem := IEEE_Rem - P; 506 else 507 P_Even := True; 508 end if; 509 510 P := P * 0.5; 511 end loop; 512 end if; 513 514 -- That completes the calculation of modulus remainder. The final 515 -- step is get the IEEE remainder. Here we need to compare Rem with 516 -- (abs Y) / 2. We must be careful of unrepresentable Y/2 value 517 -- caused by subnormal numbers 518 519 if P_Exp >= 0 then 520 A := IEEE_Rem; 521 B := abs Y * 0.5; 522 523 else 524 A := IEEE_Rem * 2.0; 525 B := abs Y; 526 end if; 527 528 if A > B or else (A = B and then not P_Even) then 529 IEEE_Rem := IEEE_Rem - abs Y; 530 end if; 531 532 return Sign_X * IEEE_Rem; 533 end Remainder; 534 535 -------------- 536 -- Rounding -- 537 -------------- 538 539 function Rounding (X : T) return T is 540 Result : T; 541 Tail : T; 542 543 begin 544 Result := Truncation (abs X); 545 Tail := abs X - Result; 546 547 if Tail >= 0.5 then 548 Result := Result + 1.0; 549 end if; 550 551 if X > 0.0 then 552 return Result; 553 554 elsif X < 0.0 then 555 return -Result; 556 557 -- For zero case, make sure sign of zero is preserved 558 559 else 560 return X; 561 end if; 562 end Rounding; 563 564 ------------- 565 -- Scaling -- 566 ------------- 567 568 -- Return x * rad ** adjustment quickly, or quietly underflow to zero, 569 -- or overflow naturally. 570 571 function Scaling (X : T; Adjustment : UI) return T is 572 begin 573 if X = 0.0 or else Adjustment = 0 then 574 return X; 575 end if; 576 577 -- Nonzero x essentially, just multiply repeatedly by Rad ** (+-2**n) 578 579 declare 580 Y : T := X; 581 Ex : UI := Adjustment; 582 583 -- Y * Rad ** Ex is invariant 584 585 begin 586 if Ex < 0 then 587 while Ex <= -Log_Power (Expbits'Last) loop 588 Y := Y * R_Neg_Power (Expbits'Last); 589 Ex := Ex + Log_Power (Expbits'Last); 590 end loop; 591 592 -- -64 < Ex <= 0 593 594 for N in reverse Expbits'First .. Expbits'Last - 1 loop 595 if Ex <= -Log_Power (N) then 596 Y := Y * R_Neg_Power (N); 597 Ex := Ex + Log_Power (N); 598 end if; 599 600 -- -Log_Power (N) < Ex <= 0 601 602 end loop; 603 604 -- Ex = 0 605 606 else 607 -- Ex >= 0 608 609 while Ex >= Log_Power (Expbits'Last) loop 610 Y := Y * R_Power (Expbits'Last); 611 Ex := Ex - Log_Power (Expbits'Last); 612 end loop; 613 614 -- 0 <= Ex < 64 615 616 for N in reverse Expbits'First .. Expbits'Last - 1 loop 617 if Ex >= Log_Power (N) then 618 Y := Y * R_Power (N); 619 Ex := Ex - Log_Power (N); 620 end if; 621 622 -- 0 <= Ex < Log_Power (N) 623 624 end loop; 625 626 -- Ex = 0 627 628 end if; 629 630 return Y; 631 end; 632 end Scaling; 633 634 ---------- 635 -- Succ -- 636 ---------- 637 638 function Succ (X : T) return T is 639 X_Frac : T; 640 X_Exp : UI; 641 X1, X2 : T; 642 643 begin 644 -- Treat zero specially since it has a zero exponent 645 646 if X = 0.0 then 647 X1 := 2.0 ** T'Machine_Emin; 648 649 -- Following loop generates smallest denormal 650 651 loop 652 X2 := T'Machine (X1 / 2.0); 653 exit when X2 = 0.0; 654 X1 := X2; 655 end loop; 656 657 return X1; 658 659 -- Special treatment for largest positive number 660 661 elsif X = T'Last then 662 663 -- If not generating infinities, we raise a constraint error 664 665 raise Constraint_Error with "Succ of largest positive number"; 666 667 -- Otherwise generate a positive infinity 668 669 -- For infinities, return unchanged 670 671 elsif X < T'First or else X > T'Last then 672 return X; 673 674 -- Add to the given number a number equivalent to the value 675 -- of its least significant bit. Given that the most significant bit 676 -- represents a value of 1.0 * radix ** (exp - 1), the value we want 677 -- is obtained by shifting this by (mantissa-1) bits to the right, 678 -- i.e. decreasing the exponent by that amount. 679 680 else 681 Decompose (X, X_Frac, X_Exp); 682 683 -- A special case, if the number we had was a negative power of two, 684 -- then we want to add half of what we would otherwise add, since the 685 -- exponent is going to be reduced. 686 687 -- Note that X_Frac has the same sign as X, so if X_Frac is -0.5, 688 -- then we know that we have a negative number (and hence a negative 689 -- power of 2). 690 691 if X_Frac = -0.5 then 692 return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa - 1); 693 694 -- Otherwise the exponent is unchanged 695 696 else 697 return X + Gradual_Scaling (X_Exp - T'Machine_Mantissa); 698 end if; 699 end if; 700 end Succ; 701 702 ---------------- 703 -- Truncation -- 704 ---------------- 705 706 -- The basic approach is to compute 707 708 -- T'Machine (RM1 + N) - RM1 709 710 -- where N >= 0.0 and RM1 = radix ** (mantissa - 1) 711 712 -- This works provided that the intermediate result (RM1 + N) does not 713 -- have extra precision (which is why we call Machine). When we compute 714 -- RM1 + N, the exponent of N will be normalized and the mantissa shifted 715 -- appropriately so the lower order bits, which cannot contribute to the 716 -- integer part of N, fall off on the right. When we subtract RM1 again, 717 -- the significant bits of N are shifted to the left, and what we have is 718 -- an integer, because only the first e bits are different from zero 719 -- (assuming binary radix here). 720 721 function Truncation (X : T) return T is 722 Result : T; 723 724 begin 725 Result := abs X; 726 727 if Result >= Radix_To_M_Minus_1 then 728 return T'Machine (X); 729 730 else 731 Result := 732 T'Machine (Radix_To_M_Minus_1 + Result) - Radix_To_M_Minus_1; 733 734 if Result > abs X then 735 Result := Result - 1.0; 736 end if; 737 738 if X > 0.0 then 739 return Result; 740 741 elsif X < 0.0 then 742 return -Result; 743 744 -- For zero case, make sure sign of zero is preserved 745 746 else 747 return X; 748 end if; 749 end if; 750 end Truncation; 751 752 ----------------------- 753 -- Unbiased_Rounding -- 754 ----------------------- 755 756 function Unbiased_Rounding (X : T) return T is 757 Abs_X : constant T := abs X; 758 Result : T; 759 Tail : T; 760 761 begin 762 Result := Truncation (Abs_X); 763 Tail := Abs_X - Result; 764 765 if Tail > 0.5 then 766 Result := Result + 1.0; 767 768 elsif Tail = 0.5 then 769 Result := 2.0 * Truncation ((Result / 2.0) + 0.5); 770 end if; 771 772 if X > 0.0 then 773 return Result; 774 775 elsif X < 0.0 then 776 return -Result; 777 778 -- For zero case, make sure sign of zero is preserved 779 780 else 781 return X; 782 end if; 783 end Unbiased_Rounding; 784 785 ----------- 786 -- Valid -- 787 ----------- 788 789 function Valid (X : not null access T) return Boolean is 790 IEEE_Emin : constant Integer := T'Machine_Emin - 1; 791 IEEE_Emax : constant Integer := T'Machine_Emax - 1; 792 793 IEEE_Bias : constant Integer := -(IEEE_Emin - 1); 794 795 subtype IEEE_Exponent_Range is 796 Integer range IEEE_Emin - 1 .. IEEE_Emax + 1; 797 798 -- The implementation of this floating point attribute uses a 799 -- representation type Float_Rep that allows direct access to the 800 -- exponent and mantissa parts of a floating point number. 801 802 -- The Float_Rep type is an array of Float_Word elements. This 803 -- representation is chosen to make it possible to size the type based 804 -- on a generic parameter. Since the array size is known at compile 805 -- time, efficient code can still be generated. The size of Float_Word 806 -- elements should be large enough to allow accessing the exponent in 807 -- one read, but small enough so that all floating point object sizes 808 -- are a multiple of the Float_Word'Size. 809 810 -- The following conditions must be met for all possible instantiations 811 -- of the attributes package: 812 813 -- - T'Size is an integral multiple of Float_Word'Size 814 815 -- - The exponent and sign are completely contained in a single 816 -- component of Float_Rep, named Most_Significant_Word (MSW). 817 818 -- - The sign occupies the most significant bit of the MSW and the 819 -- exponent is in the following bits. Unused bits (if any) are in 820 -- the least significant part. 821 822 type Float_Word is mod 2**Positive'Min (System.Word_Size, 32); 823 type Rep_Index is range 0 .. 7; 824 825 Rep_Words : constant Positive := 826 (T'Size + Float_Word'Size - 1) / Float_Word'Size; 827 Rep_Last : constant Rep_Index := 828 Rep_Index'Min 829 (Rep_Index (Rep_Words - 1), 830 (T'Mantissa + 16) / Float_Word'Size); 831 -- Determine the number of Float_Words needed for representing the 832 -- entire floating-point value. Do not take into account excessive 833 -- padding, as occurs on IA-64 where 80 bits floats get padded to 128 834 -- bits. In general, the exponent field cannot be larger than 15 bits, 835 -- even for 128-bit floating-point types, so the final format size 836 -- won't be larger than T'Mantissa + 16. 837 838 type Float_Rep is 839 array (Rep_Index range 0 .. Rep_Index (Rep_Words - 1)) of Float_Word; 840 841 pragma Suppress_Initialization (Float_Rep); 842 -- This pragma suppresses the generation of an initialization procedure 843 -- for type Float_Rep when operating in Initialize/Normalize_Scalars 844 -- mode. This is not just a matter of efficiency, but of functionality, 845 -- since Valid has a pragma Inline_Always, which is not permitted if 846 -- there are nested subprograms present. 847 848 Most_Significant_Word : constant Rep_Index := 849 Rep_Last * Standard'Default_Bit_Order; 850 -- Finding the location of the Exponent_Word is a bit tricky. In general 851 -- we assume Word_Order = Bit_Order. 852 853 Exponent_Factor : constant Float_Word := 854 2**(Float_Word'Size - 1) / 855 Float_Word (IEEE_Emax - IEEE_Emin + 3) * 856 Boolean'Pos (Most_Significant_Word /= 2) + 857 Boolean'Pos (Most_Significant_Word = 2); 858 -- Factor that the extracted exponent needs to be divided by to be in 859 -- range 0 .. IEEE_Emax - IEEE_Emin + 2. Special case: Exponent_Factor 860 -- is 1 for x86/IA64 double extended (GCC adds unused bits to the type). 861 862 Exponent_Mask : constant Float_Word := 863 Float_Word (IEEE_Emax - IEEE_Emin + 2) * 864 Exponent_Factor; 865 -- Value needed to mask out the exponent field. This assumes that the 866 -- range IEEE_Emin - 1 .. IEEE_Emax + contains 2**N values, for some N 867 -- in Natural. 868 869 function To_Float is new Ada.Unchecked_Conversion (Float_Rep, T); 870 871 type Float_Access is access all T; 872 function To_Address is 873 new Ada.Unchecked_Conversion (Float_Access, System.Address); 874 875 XA : constant System.Address := To_Address (Float_Access (X)); 876 877 R : Float_Rep; 878 pragma Import (Ada, R); 879 for R'Address use XA; 880 -- R is a view of the input floating-point parameter. Note that we 881 -- must avoid copying the actual bits of this parameter in float 882 -- form (since it may be a signalling NaN). 883 884 E : constant IEEE_Exponent_Range := 885 Integer ((R (Most_Significant_Word) and Exponent_Mask) / 886 Exponent_Factor) 887 - IEEE_Bias; 888 -- Mask/Shift T to only get bits from the exponent. Then convert biased 889 -- value to integer value. 890 891 SR : Float_Rep; 892 -- Float_Rep representation of significant of X.all 893 894 begin 895 if T'Denorm then 896 897 -- All denormalized numbers are valid, so the only invalid numbers 898 -- are overflows and NaNs, both with exponent = Emax + 1. 899 900 return E /= IEEE_Emax + 1; 901 902 end if; 903 904 -- All denormalized numbers except 0.0 are invalid 905 906 -- Set exponent of X to zero, so we end up with the significand, which 907 -- definitely is a valid number and can be converted back to a float. 908 909 SR := R; 910 SR (Most_Significant_Word) := 911 (SR (Most_Significant_Word) 912 and not Exponent_Mask) + Float_Word (IEEE_Bias) * Exponent_Factor; 913 914 return (E in IEEE_Emin .. IEEE_Emax) or else 915 ((E = IEEE_Emin - 1) and then abs To_Float (SR) = 1.0); 916 end Valid; 917 918end System.Fat_Gen; 919