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