1------------------------------------------------------------------------------ 2-- -- 3-- GNAT RUNTIME COMPONENTS -- 4-- -- 5-- ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS -- 6-- -- 7-- B o d y -- 8-- -- 9-- Copyright (C) 1992-2002, 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 2, 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. See the GNU General Public License -- 17-- for more details. You should have received a copy of the GNU General -- 18-- Public License distributed with GNAT; see file COPYING. If not, write -- 19-- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, -- 20-- MA 02111-1307, USA. -- 21-- -- 22-- As a special exception, if other files instantiate generics from this -- 23-- unit, or you link this unit with other files to produce an executable, -- 24-- this unit does not by itself cause the resulting executable to be -- 25-- covered by the GNU General Public License. This exception does not -- 26-- however invalidate any other reasons why the executable file might be -- 27-- covered by the GNU Public License. -- 28-- -- 29-- GNAT was originally developed by the GNAT team at New York University. -- 30-- Extensive contributions were provided by Ada Core Technologies Inc. -- 31-- -- 32------------------------------------------------------------------------------ 33 34-- This body is specifically for using an Ada interface to C math.h to get 35-- the computation engine. Many special cases are handled locally to avoid 36-- unnecessary calls. This is not a "strict" implementation, but takes full 37-- advantage of the C functions, e.g. in providing interface to hardware 38-- provided versions of the elementary functions. 39 40-- Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan, 41-- sinh, cosh, tanh from C library via math.h 42 43with Ada.Numerics.Aux; 44 45package body Ada.Numerics.Generic_Elementary_Functions is 46 47 use type Ada.Numerics.Aux.Double; 48 49 Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696; 50 Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755; 51 Half_Log_Two : constant := Log_Two / 2; 52 53 subtype T is Float_Type'Base; 54 subtype Double is Aux.Double; 55 56 Two_Pi : constant T := 2.0 * Pi; 57 Half_Pi : constant T := Pi / 2.0; 58 59 Half_Log_Epsilon : constant T := T (1 - T'Model_Mantissa) * Half_Log_Two; 60 Log_Inverse_Epsilon : constant T := T (T'Model_Mantissa - 1) * Log_Two; 61 Sqrt_Epsilon : constant T := Sqrt_Two ** (1 - T'Model_Mantissa); 62 63 ----------------------- 64 -- Local Subprograms -- 65 ----------------------- 66 67 function Exp_Strict (X : Float_Type'Base) return Float_Type'Base; 68 -- Cody/Waite routine, supposedly more precise than the library 69 -- version. Currently only needed for Sinh/Cosh on X86 with the largest 70 -- FP type. 71 72 function Local_Atan 73 (Y : Float_Type'Base; 74 X : Float_Type'Base := 1.0) 75 return Float_Type'Base; 76 -- Common code for arc tangent after cyele reduction 77 78 ---------- 79 -- "**" -- 80 ---------- 81 82 function "**" (Left, Right : Float_Type'Base) return Float_Type'Base is 83 A_Right : Float_Type'Base; 84 Int_Part : Integer; 85 Result : Float_Type'Base; 86 R1 : Float_Type'Base; 87 Rest : Float_Type'Base; 88 89 begin 90 if Left = 0.0 91 and then Right = 0.0 92 then 93 raise Argument_Error; 94 95 elsif Left < 0.0 then 96 raise Argument_Error; 97 98 elsif Right = 0.0 then 99 return 1.0; 100 101 elsif Left = 0.0 then 102 if Right < 0.0 then 103 raise Constraint_Error; 104 else 105 return 0.0; 106 end if; 107 108 elsif Left = 1.0 then 109 return 1.0; 110 111 elsif Right = 1.0 then 112 return Left; 113 114 else 115 begin 116 if Right = 2.0 then 117 return Left * Left; 118 119 elsif Right = 0.5 then 120 return Sqrt (Left); 121 122 else 123 A_Right := abs (Right); 124 125 -- If exponent is larger than one, compute integer exponen- 126 -- tiation if possible, and evaluate fractional part with 127 -- more precision. The relative error is now proportional 128 -- to the fractional part of the exponent only. 129 130 if A_Right > 1.0 131 and then A_Right < Float_Type'Base (Integer'Last) 132 then 133 Int_Part := Integer (Float_Type'Base'Truncation (A_Right)); 134 Result := Left ** Int_Part; 135 Rest := A_Right - Float_Type'Base (Int_Part); 136 137 -- Compute with two leading bits of the mantissa using 138 -- square roots. Bound to be better than logarithms, and 139 -- easily extended to greater precision. 140 141 if Rest >= 0.5 then 142 R1 := Sqrt (Left); 143 Result := Result * R1; 144 Rest := Rest - 0.5; 145 146 if Rest >= 0.25 then 147 Result := Result * Sqrt (R1); 148 Rest := Rest - 0.25; 149 end if; 150 151 elsif Rest >= 0.25 then 152 Result := Result * Sqrt (Sqrt (Left)); 153 Rest := Rest - 0.25; 154 end if; 155 156 Result := Result * 157 Float_Type'Base (Aux.Pow (Double (Left), Double (Rest))); 158 159 if Right >= 0.0 then 160 return Result; 161 else 162 return (1.0 / Result); 163 end if; 164 else 165 return 166 Float_Type'Base (Aux.Pow (Double (Left), Double (Right))); 167 end if; 168 end if; 169 170 exception 171 when others => 172 raise Constraint_Error; 173 end; 174 end if; 175 end "**"; 176 177 ------------ 178 -- Arccos -- 179 ------------ 180 181 -- Natural cycle 182 183 function Arccos (X : Float_Type'Base) return Float_Type'Base is 184 Temp : Float_Type'Base; 185 186 begin 187 if abs X > 1.0 then 188 raise Argument_Error; 189 190 elsif abs X < Sqrt_Epsilon then 191 return Pi / 2.0 - X; 192 193 elsif X = 1.0 then 194 return 0.0; 195 196 elsif X = -1.0 then 197 return Pi; 198 end if; 199 200 Temp := Float_Type'Base (Aux.Acos (Double (X))); 201 202 if Temp < 0.0 then 203 Temp := Pi + Temp; 204 end if; 205 206 return Temp; 207 end Arccos; 208 209 -- Arbitrary cycle 210 211 function Arccos (X, Cycle : Float_Type'Base) return Float_Type'Base is 212 Temp : Float_Type'Base; 213 214 begin 215 if Cycle <= 0.0 then 216 raise Argument_Error; 217 218 elsif abs X > 1.0 then 219 raise Argument_Error; 220 221 elsif abs X < Sqrt_Epsilon then 222 return Cycle / 4.0; 223 224 elsif X = 1.0 then 225 return 0.0; 226 227 elsif X = -1.0 then 228 return Cycle / 2.0; 229 end if; 230 231 Temp := Arctan (Sqrt ((1.0 - X) * (1.0 + X)) / X, 1.0, Cycle); 232 233 if Temp < 0.0 then 234 Temp := Cycle / 2.0 + Temp; 235 end if; 236 237 return Temp; 238 end Arccos; 239 240 ------------- 241 -- Arccosh -- 242 ------------- 243 244 function Arccosh (X : Float_Type'Base) return Float_Type'Base is 245 begin 246 -- Return positive branch of Log (X - Sqrt (X * X - 1.0)), or 247 -- the proper approximation for X close to 1 or >> 1. 248 249 if X < 1.0 then 250 raise Argument_Error; 251 252 elsif X < 1.0 + Sqrt_Epsilon then 253 return Sqrt (2.0 * (X - 1.0)); 254 255 elsif X > 1.0 / Sqrt_Epsilon then 256 return Log (X) + Log_Two; 257 258 else 259 return Log (X + Sqrt ((X - 1.0) * (X + 1.0))); 260 end if; 261 end Arccosh; 262 263 ------------ 264 -- Arccot -- 265 ------------ 266 267 -- Natural cycle 268 269 function Arccot 270 (X : Float_Type'Base; 271 Y : Float_Type'Base := 1.0) 272 return Float_Type'Base 273 is 274 begin 275 -- Just reverse arguments 276 277 return Arctan (Y, X); 278 end Arccot; 279 280 -- Arbitrary cycle 281 282 function Arccot 283 (X : Float_Type'Base; 284 Y : Float_Type'Base := 1.0; 285 Cycle : Float_Type'Base) 286 return Float_Type'Base 287 is 288 begin 289 -- Just reverse arguments 290 291 return Arctan (Y, X, Cycle); 292 end Arccot; 293 294 ------------- 295 -- Arccoth -- 296 ------------- 297 298 function Arccoth (X : Float_Type'Base) return Float_Type'Base is 299 begin 300 if abs X > 2.0 then 301 return Arctanh (1.0 / X); 302 303 elsif abs X = 1.0 then 304 raise Constraint_Error; 305 306 elsif abs X < 1.0 then 307 raise Argument_Error; 308 309 else 310 -- 1.0 < abs X <= 2.0. One of X + 1.0 and X - 1.0 is exact, the 311 -- other has error 0 or Epsilon. 312 313 return 0.5 * (Log (abs (X + 1.0)) - Log (abs (X - 1.0))); 314 end if; 315 end Arccoth; 316 317 ------------ 318 -- Arcsin -- 319 ------------ 320 321 -- Natural cycle 322 323 function Arcsin (X : Float_Type'Base) return Float_Type'Base is 324 begin 325 if abs X > 1.0 then 326 raise Argument_Error; 327 328 elsif abs X < Sqrt_Epsilon then 329 return X; 330 331 elsif X = 1.0 then 332 return Pi / 2.0; 333 334 elsif X = -1.0 then 335 return -Pi / 2.0; 336 end if; 337 338 return Float_Type'Base (Aux.Asin (Double (X))); 339 end Arcsin; 340 341 -- Arbitrary cycle 342 343 function Arcsin (X, Cycle : Float_Type'Base) return Float_Type'Base is 344 begin 345 if Cycle <= 0.0 then 346 raise Argument_Error; 347 348 elsif abs X > 1.0 then 349 raise Argument_Error; 350 351 elsif X = 0.0 then 352 return X; 353 354 elsif X = 1.0 then 355 return Cycle / 4.0; 356 357 elsif X = -1.0 then 358 return -Cycle / 4.0; 359 end if; 360 361 return Arctan (X / Sqrt ((1.0 - X) * (1.0 + X)), 1.0, Cycle); 362 end Arcsin; 363 364 ------------- 365 -- Arcsinh -- 366 ------------- 367 368 function Arcsinh (X : Float_Type'Base) return Float_Type'Base is 369 begin 370 if abs X < Sqrt_Epsilon then 371 return X; 372 373 elsif X > 1.0 / Sqrt_Epsilon then 374 return Log (X) + Log_Two; 375 376 elsif X < -1.0 / Sqrt_Epsilon then 377 return -(Log (-X) + Log_Two); 378 379 elsif X < 0.0 then 380 return -Log (abs X + Sqrt (X * X + 1.0)); 381 382 else 383 return Log (X + Sqrt (X * X + 1.0)); 384 end if; 385 end Arcsinh; 386 387 ------------ 388 -- Arctan -- 389 ------------ 390 391 -- Natural cycle 392 393 function Arctan 394 (Y : Float_Type'Base; 395 X : Float_Type'Base := 1.0) 396 return Float_Type'Base 397 is 398 begin 399 if X = 0.0 400 and then Y = 0.0 401 then 402 raise Argument_Error; 403 404 elsif Y = 0.0 then 405 if X > 0.0 then 406 return 0.0; 407 else -- X < 0.0 408 return Pi * Float_Type'Copy_Sign (1.0, Y); 409 end if; 410 411 elsif X = 0.0 then 412 if Y > 0.0 then 413 return Half_Pi; 414 else -- Y < 0.0 415 return -Half_Pi; 416 end if; 417 418 else 419 return Local_Atan (Y, X); 420 end if; 421 end Arctan; 422 423 -- Arbitrary cycle 424 425 function Arctan 426 (Y : Float_Type'Base; 427 X : Float_Type'Base := 1.0; 428 Cycle : Float_Type'Base) 429 return Float_Type'Base 430 is 431 begin 432 if Cycle <= 0.0 then 433 raise Argument_Error; 434 435 elsif X = 0.0 436 and then Y = 0.0 437 then 438 raise Argument_Error; 439 440 elsif Y = 0.0 then 441 if X > 0.0 then 442 return 0.0; 443 else -- X < 0.0 444 return Cycle / 2.0 * Float_Type'Copy_Sign (1.0, Y); 445 end if; 446 447 elsif X = 0.0 then 448 if Y > 0.0 then 449 return Cycle / 4.0; 450 else -- Y < 0.0 451 return -Cycle / 4.0; 452 end if; 453 454 else 455 return Local_Atan (Y, X) * Cycle / Two_Pi; 456 end if; 457 end Arctan; 458 459 ------------- 460 -- Arctanh -- 461 ------------- 462 463 function Arctanh (X : Float_Type'Base) return Float_Type'Base is 464 A, B, D, A_Plus_1, A_From_1 : Float_Type'Base; 465 Mantissa : constant Integer := Float_Type'Base'Machine_Mantissa; 466 467 begin 468 -- The naive formula: 469 470 -- Arctanh (X) := (1/2) * Log (1 + X) / (1 - X) 471 472 -- is not well-behaved numerically when X < 0.5 and when X is close 473 -- to one. The following is accurate but probably not optimal. 474 475 if abs X = 1.0 then 476 raise Constraint_Error; 477 478 elsif abs X >= 1.0 - 2.0 ** (-Mantissa) then 479 480 if abs X >= 1.0 then 481 raise Argument_Error; 482 else 483 484 -- The one case that overflows if put through the method below: 485 -- abs X = 1.0 - Epsilon. In this case (1/2) log (2/Epsilon) is 486 -- accurate. This simplifies to: 487 488 return Float_Type'Copy_Sign ( 489 Half_Log_Two * Float_Type'Base (Mantissa + 1), X); 490 end if; 491 492 -- elsif abs X <= 0.5 then 493 -- why is above line commented out ??? 494 495 else 496 -- Use several piecewise linear approximations. 497 -- A is close to X, chosen so 1.0 + A, 1.0 - A, and X - A are exact. 498 -- The two scalings remove the low-order bits of X. 499 500 A := Float_Type'Base'Scaling ( 501 Float_Type'Base (Long_Long_Integer 502 (Float_Type'Base'Scaling (X, Mantissa - 1))), 1 - Mantissa); 503 504 B := X - A; -- This is exact; abs B <= 2**(-Mantissa). 505 A_Plus_1 := 1.0 + A; -- This is exact. 506 A_From_1 := 1.0 - A; -- Ditto. 507 D := A_Plus_1 * A_From_1; -- 1 - A*A. 508 509 -- use one term of the series expansion: 510 -- f (x + e) = f(x) + e * f'(x) + .. 511 512 -- The derivative of Arctanh at A is 1/(1-A*A). Next term is 513 -- A*(B/D)**2 (if a quadratic approximation is ever needed). 514 515 return 0.5 * (Log (A_Plus_1) - Log (A_From_1)) + B / D; 516 517 -- else 518 -- return 0.5 * Log ((X + 1.0) / (1.0 - X)); 519 -- why are above lines commented out ??? 520 end if; 521 end Arctanh; 522 523 --------- 524 -- Cos -- 525 --------- 526 527 -- Natural cycle 528 529 function Cos (X : Float_Type'Base) return Float_Type'Base is 530 begin 531 if X = 0.0 then 532 return 1.0; 533 534 elsif abs X < Sqrt_Epsilon then 535 return 1.0; 536 537 end if; 538 539 return Float_Type'Base (Aux.Cos (Double (X))); 540 end Cos; 541 542 -- Arbitrary cycle 543 544 function Cos (X, Cycle : Float_Type'Base) return Float_Type'Base is 545 begin 546 -- Just reuse the code for Sin. The potential small 547 -- loss of speed is negligible with proper (front-end) inlining. 548 549 return -Sin (abs X - Cycle * 0.25, Cycle); 550 end Cos; 551 552 ---------- 553 -- Cosh -- 554 ---------- 555 556 function Cosh (X : Float_Type'Base) return Float_Type'Base is 557 Lnv : constant Float_Type'Base := 8#0.542714#; 558 V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4; 559 Y : constant Float_Type'Base := abs X; 560 Z : Float_Type'Base; 561 562 begin 563 if Y < Sqrt_Epsilon then 564 return 1.0; 565 566 elsif Y > Log_Inverse_Epsilon then 567 Z := Exp_Strict (Y - Lnv); 568 return (Z + V2minus1 * Z); 569 570 else 571 Z := Exp_Strict (Y); 572 return 0.5 * (Z + 1.0 / Z); 573 end if; 574 575 end Cosh; 576 577 --------- 578 -- Cot -- 579 --------- 580 581 -- Natural cycle 582 583 function Cot (X : Float_Type'Base) return Float_Type'Base is 584 begin 585 if X = 0.0 then 586 raise Constraint_Error; 587 588 elsif abs X < Sqrt_Epsilon then 589 return 1.0 / X; 590 end if; 591 592 return 1.0 / Float_Type'Base (Aux.Tan (Double (X))); 593 end Cot; 594 595 -- Arbitrary cycle 596 597 function Cot (X, Cycle : Float_Type'Base) return Float_Type'Base is 598 T : Float_Type'Base; 599 600 begin 601 if Cycle <= 0.0 then 602 raise Argument_Error; 603 end if; 604 605 T := Float_Type'Base'Remainder (X, Cycle); 606 607 if T = 0.0 or abs T = 0.5 * Cycle then 608 raise Constraint_Error; 609 610 elsif abs T < Sqrt_Epsilon then 611 return 1.0 / T; 612 613 elsif abs T = 0.25 * Cycle then 614 return 0.0; 615 616 else 617 T := T / Cycle * Two_Pi; 618 return Cos (T) / Sin (T); 619 end if; 620 end Cot; 621 622 ---------- 623 -- Coth -- 624 ---------- 625 626 function Coth (X : Float_Type'Base) return Float_Type'Base is 627 begin 628 if X = 0.0 then 629 raise Constraint_Error; 630 631 elsif X < Half_Log_Epsilon then 632 return -1.0; 633 634 elsif X > -Half_Log_Epsilon then 635 return 1.0; 636 637 elsif abs X < Sqrt_Epsilon then 638 return 1.0 / X; 639 end if; 640 641 return 1.0 / Float_Type'Base (Aux.Tanh (Double (X))); 642 end Coth; 643 644 --------- 645 -- Exp -- 646 --------- 647 648 function Exp (X : Float_Type'Base) return Float_Type'Base is 649 Result : Float_Type'Base; 650 651 begin 652 if X = 0.0 then 653 return 1.0; 654 end if; 655 656 Result := Float_Type'Base (Aux.Exp (Double (X))); 657 658 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows 659 -- is False, then we can just leave it as an infinity (and indeed we 660 -- prefer to do so). But if Machine_Overflows is True, then we have 661 -- to raise a Constraint_Error exception as required by the RM. 662 663 if Float_Type'Machine_Overflows and then not Result'Valid then 664 raise Constraint_Error; 665 end if; 666 667 return Result; 668 end Exp; 669 670 ---------------- 671 -- Exp_Strict -- 672 ---------------- 673 674 function Exp_Strict (X : Float_Type'Base) return Float_Type'Base is 675 G : Float_Type'Base; 676 Z : Float_Type'Base; 677 678 P0 : constant := 0.25000_00000_00000_00000; 679 P1 : constant := 0.75753_18015_94227_76666E-2; 680 P2 : constant := 0.31555_19276_56846_46356E-4; 681 682 Q0 : constant := 0.5; 683 Q1 : constant := 0.56817_30269_85512_21787E-1; 684 Q2 : constant := 0.63121_89437_43985_02557E-3; 685 Q3 : constant := 0.75104_02839_98700_46114E-6; 686 687 C1 : constant := 8#0.543#; 688 C2 : constant := -2.1219_44400_54690_58277E-4; 689 Le : constant := 1.4426_95040_88896_34074; 690 691 XN : Float_Type'Base; 692 P, Q, R : Float_Type'Base; 693 694 begin 695 if X = 0.0 then 696 return 1.0; 697 end if; 698 699 XN := Float_Type'Base'Rounding (X * Le); 700 G := (X - XN * C1) - XN * C2; 701 Z := G * G; 702 P := G * ((P2 * Z + P1) * Z + P0); 703 Q := ((Q3 * Z + Q2) * Z + Q1) * Z + Q0; 704 R := 0.5 + P / (Q - P); 705 706 R := Float_Type'Base'Scaling (R, Integer (XN) + 1); 707 708 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows 709 -- is False, then we can just leave it as an infinity (and indeed we 710 -- prefer to do so). But if Machine_Overflows is True, then we have 711 -- to raise a Constraint_Error exception as required by the RM. 712 713 if Float_Type'Machine_Overflows and then not R'Valid then 714 raise Constraint_Error; 715 else 716 return R; 717 end if; 718 719 end Exp_Strict; 720 721 ---------------- 722 -- Local_Atan -- 723 ---------------- 724 725 function Local_Atan 726 (Y : Float_Type'Base; 727 X : Float_Type'Base := 1.0) 728 return Float_Type'Base 729 is 730 Z : Float_Type'Base; 731 Raw_Atan : Float_Type'Base; 732 733 begin 734 if abs Y > abs X then 735 Z := abs (X / Y); 736 else 737 Z := abs (Y / X); 738 end if; 739 740 if Z < Sqrt_Epsilon then 741 Raw_Atan := Z; 742 743 elsif Z = 1.0 then 744 Raw_Atan := Pi / 4.0; 745 746 else 747 Raw_Atan := Float_Type'Base (Aux.Atan (Double (Z))); 748 end if; 749 750 if abs Y > abs X then 751 Raw_Atan := Half_Pi - Raw_Atan; 752 end if; 753 754 if X > 0.0 then 755 if Y > 0.0 then 756 return Raw_Atan; 757 else -- Y < 0.0 758 return -Raw_Atan; 759 end if; 760 761 else -- X < 0.0 762 if Y > 0.0 then 763 return Pi - Raw_Atan; 764 else -- Y < 0.0 765 return -(Pi - Raw_Atan); 766 end if; 767 end if; 768 end Local_Atan; 769 770 --------- 771 -- Log -- 772 --------- 773 774 -- Natural base 775 776 function Log (X : Float_Type'Base) return Float_Type'Base is 777 begin 778 if X < 0.0 then 779 raise Argument_Error; 780 781 elsif X = 0.0 then 782 raise Constraint_Error; 783 784 elsif X = 1.0 then 785 return 0.0; 786 end if; 787 788 return Float_Type'Base (Aux.Log (Double (X))); 789 end Log; 790 791 -- Arbitrary base 792 793 function Log (X, Base : Float_Type'Base) return Float_Type'Base is 794 begin 795 if X < 0.0 then 796 raise Argument_Error; 797 798 elsif Base <= 0.0 or else Base = 1.0 then 799 raise Argument_Error; 800 801 elsif X = 0.0 then 802 raise Constraint_Error; 803 804 elsif X = 1.0 then 805 return 0.0; 806 end if; 807 808 return Float_Type'Base (Aux.Log (Double (X)) / Aux.Log (Double (Base))); 809 end Log; 810 811 --------- 812 -- Sin -- 813 --------- 814 815 -- Natural cycle 816 817 function Sin (X : Float_Type'Base) return Float_Type'Base is 818 begin 819 if abs X < Sqrt_Epsilon then 820 return X; 821 end if; 822 823 return Float_Type'Base (Aux.Sin (Double (X))); 824 end Sin; 825 826 -- Arbitrary cycle 827 828 function Sin (X, Cycle : Float_Type'Base) return Float_Type'Base is 829 T : Float_Type'Base; 830 831 begin 832 if Cycle <= 0.0 then 833 raise Argument_Error; 834 835 elsif X = 0.0 then 836 -- Is this test really needed on any machine ??? 837 return X; 838 end if; 839 840 T := Float_Type'Base'Remainder (X, Cycle); 841 842 -- The following two reductions reduce the argument 843 -- to the interval [-0.25 * Cycle, 0.25 * Cycle]. 844 -- This reduction is exact and is needed to prevent 845 -- inaccuracy that may result if the sinus function 846 -- a different (more accurate) value of Pi in its 847 -- reduction than is used in the multiplication with Two_Pi. 848 849 if abs T > 0.25 * Cycle then 850 T := 0.5 * Float_Type'Copy_Sign (Cycle, T) - T; 851 end if; 852 853 -- Could test for 12.0 * abs T = Cycle, and return 854 -- an exact value in those cases. It is not clear that 855 -- this is worth the extra test though. 856 857 return Float_Type'Base (Aux.Sin (Double (T / Cycle * Two_Pi))); 858 end Sin; 859 860 ---------- 861 -- Sinh -- 862 ---------- 863 864 function Sinh (X : Float_Type'Base) return Float_Type'Base is 865 Lnv : constant Float_Type'Base := 8#0.542714#; 866 V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4; 867 Y : constant Float_Type'Base := abs X; 868 F : constant Float_Type'Base := Y * Y; 869 Z : Float_Type'Base; 870 871 Float_Digits_1_6 : constant Boolean := Float_Type'Digits < 7; 872 873 begin 874 if Y < Sqrt_Epsilon then 875 return X; 876 877 elsif Y > Log_Inverse_Epsilon then 878 Z := Exp_Strict (Y - Lnv); 879 Z := Z + V2minus1 * Z; 880 881 elsif Y < 1.0 then 882 883 if Float_Digits_1_6 then 884 885 -- Use expansion provided by Cody and Waite, p. 226. Note that 886 -- leading term of the polynomial in Q is exactly 1.0. 887 888 declare 889 P0 : constant := -0.71379_3159E+1; 890 P1 : constant := -0.19033_3399E+0; 891 Q0 : constant := -0.42827_7109E+2; 892 893 begin 894 Z := Y + Y * F * (P1 * F + P0) / (F + Q0); 895 end; 896 897 else 898 declare 899 P0 : constant := -0.35181_28343_01771_17881E+6; 900 P1 : constant := -0.11563_52119_68517_68270E+5; 901 P2 : constant := -0.16375_79820_26307_51372E+3; 902 P3 : constant := -0.78966_12741_73570_99479E+0; 903 Q0 : constant := -0.21108_77005_81062_71242E+7; 904 Q1 : constant := 0.36162_72310_94218_36460E+5; 905 Q2 : constant := -0.27773_52311_96507_01667E+3; 906 907 begin 908 Z := Y + Y * F * (((P3 * F + P2) * F + P1) * F + P0) 909 / (((F + Q2) * F + Q1) * F + Q0); 910 end; 911 end if; 912 913 else 914 Z := Exp_Strict (Y); 915 Z := 0.5 * (Z - 1.0 / Z); 916 end if; 917 918 if X > 0.0 then 919 return Z; 920 else 921 return -Z; 922 end if; 923 end Sinh; 924 925 ---------- 926 -- Sqrt -- 927 ---------- 928 929 function Sqrt (X : Float_Type'Base) return Float_Type'Base is 930 begin 931 if X < 0.0 then 932 raise Argument_Error; 933 934 -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE 935 936 elsif X = 0.0 then 937 return X; 938 939 end if; 940 941 return Float_Type'Base (Aux.Sqrt (Double (X))); 942 end Sqrt; 943 944 --------- 945 -- Tan -- 946 --------- 947 948 -- Natural cycle 949 950 function Tan (X : Float_Type'Base) return Float_Type'Base is 951 begin 952 if abs X < Sqrt_Epsilon then 953 return X; 954 955 elsif abs X = Pi / 2.0 then 956 raise Constraint_Error; 957 end if; 958 959 return Float_Type'Base (Aux.Tan (Double (X))); 960 end Tan; 961 962 -- Arbitrary cycle 963 964 function Tan (X, Cycle : Float_Type'Base) return Float_Type'Base is 965 T : Float_Type'Base; 966 967 begin 968 if Cycle <= 0.0 then 969 raise Argument_Error; 970 971 elsif X = 0.0 then 972 return X; 973 end if; 974 975 T := Float_Type'Base'Remainder (X, Cycle); 976 977 if abs T = 0.25 * Cycle then 978 raise Constraint_Error; 979 980 elsif abs T = 0.5 * Cycle then 981 return 0.0; 982 983 else 984 T := T / Cycle * Two_Pi; 985 return Sin (T) / Cos (T); 986 end if; 987 988 end Tan; 989 990 ---------- 991 -- Tanh -- 992 ---------- 993 994 function Tanh (X : Float_Type'Base) return Float_Type'Base is 995 P0 : constant Float_Type'Base := -0.16134_11902E4; 996 P1 : constant Float_Type'Base := -0.99225_92967E2; 997 P2 : constant Float_Type'Base := -0.96437_49299E0; 998 999 Q0 : constant Float_Type'Base := 0.48402_35707E4; 1000 Q1 : constant Float_Type'Base := 0.22337_72071E4; 1001 Q2 : constant Float_Type'Base := 0.11274_47438E3; 1002 Q3 : constant Float_Type'Base := 0.10000000000E1; 1003 1004 Half_Ln3 : constant Float_Type'Base := 0.54930_61443; 1005 1006 P, Q, R : Float_Type'Base; 1007 Y : constant Float_Type'Base := abs X; 1008 G : constant Float_Type'Base := Y * Y; 1009 1010 Float_Type_Digits_15_Or_More : constant Boolean := 1011 Float_Type'Digits > 14; 1012 1013 begin 1014 if X < Half_Log_Epsilon then 1015 return -1.0; 1016 1017 elsif X > -Half_Log_Epsilon then 1018 return 1.0; 1019 1020 elsif Y < Sqrt_Epsilon then 1021 return X; 1022 1023 elsif Y < Half_Ln3 1024 and then Float_Type_Digits_15_Or_More 1025 then 1026 P := (P2 * G + P1) * G + P0; 1027 Q := ((Q3 * G + Q2) * G + Q1) * G + Q0; 1028 R := G * (P / Q); 1029 return X + X * R; 1030 1031 else 1032 return Float_Type'Base (Aux.Tanh (Double (X))); 1033 end if; 1034 end Tanh; 1035 1036end Ada.Numerics.Generic_Elementary_Functions; 1037