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