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