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-2012, 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 X = 0.0 then 513 return 1.0; 514 515 elsif abs X < Sqrt_Epsilon then 516 return 1.0; 517 518 end if; 519 520 return Float_Type'Base (Aux.Cos (Double (X))); 521 end Cos; 522 523 -- Arbitrary cycle 524 525 function Cos (X, Cycle : Float_Type'Base) return Float_Type'Base is 526 begin 527 -- Just reuse the code for Sin. The potential small loss of speed is 528 -- negligible with proper (front-end) inlining. 529 530 return -Sin (abs X - Cycle * 0.25, Cycle); 531 end Cos; 532 533 ---------- 534 -- Cosh -- 535 ---------- 536 537 function Cosh (X : Float_Type'Base) return Float_Type'Base is 538 Lnv : constant Float_Type'Base := 8#0.542714#; 539 V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4; 540 Y : constant Float_Type'Base := abs X; 541 Z : Float_Type'Base; 542 543 begin 544 if Y < Sqrt_Epsilon then 545 return 1.0; 546 547 elsif Y > Log_Inverse_Epsilon then 548 Z := Exp_Strict (Y - Lnv); 549 return (Z + V2minus1 * Z); 550 551 else 552 Z := Exp_Strict (Y); 553 return 0.5 * (Z + 1.0 / Z); 554 end if; 555 556 end Cosh; 557 558 --------- 559 -- Cot -- 560 --------- 561 562 -- Natural cycle 563 564 function Cot (X : Float_Type'Base) return Float_Type'Base is 565 begin 566 if X = 0.0 then 567 raise Constraint_Error; 568 569 elsif abs X < Sqrt_Epsilon then 570 return 1.0 / X; 571 end if; 572 573 return 1.0 / Float_Type'Base (Aux.Tan (Double (X))); 574 end Cot; 575 576 -- Arbitrary cycle 577 578 function Cot (X, Cycle : Float_Type'Base) return Float_Type'Base is 579 T : Float_Type'Base; 580 581 begin 582 if Cycle <= 0.0 then 583 raise Argument_Error; 584 end if; 585 586 T := Float_Type'Base'Remainder (X, Cycle); 587 588 if T = 0.0 or else abs T = 0.5 * Cycle then 589 raise Constraint_Error; 590 591 elsif abs T < Sqrt_Epsilon then 592 return 1.0 / T; 593 594 elsif abs T = 0.25 * Cycle then 595 return 0.0; 596 597 else 598 T := T / Cycle * Two_Pi; 599 return Cos (T) / Sin (T); 600 end if; 601 end Cot; 602 603 ---------- 604 -- Coth -- 605 ---------- 606 607 function Coth (X : Float_Type'Base) return Float_Type'Base is 608 begin 609 if X = 0.0 then 610 raise Constraint_Error; 611 612 elsif X < Half_Log_Epsilon then 613 return -1.0; 614 615 elsif X > -Half_Log_Epsilon then 616 return 1.0; 617 618 elsif abs X < Sqrt_Epsilon then 619 return 1.0 / X; 620 end if; 621 622 return 1.0 / Float_Type'Base (Aux.Tanh (Double (X))); 623 end Coth; 624 625 --------- 626 -- Exp -- 627 --------- 628 629 function Exp (X : Float_Type'Base) return Float_Type'Base is 630 Result : Float_Type'Base; 631 632 begin 633 if X = 0.0 then 634 return 1.0; 635 end if; 636 637 Result := Float_Type'Base (Aux.Exp (Double (X))); 638 639 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows 640 -- is False, then we can just leave it as an infinity (and indeed we 641 -- prefer to do so). But if Machine_Overflows is True, then we have 642 -- to raise a Constraint_Error exception as required by the RM. 643 644 if Float_Type'Machine_Overflows and then not Result'Valid then 645 raise Constraint_Error; 646 end if; 647 648 return Result; 649 end Exp; 650 651 ---------------- 652 -- Exp_Strict -- 653 ---------------- 654 655 function Exp_Strict (X : Float_Type'Base) return Float_Type'Base is 656 G : Float_Type'Base; 657 Z : Float_Type'Base; 658 659 P0 : constant := 0.25000_00000_00000_00000; 660 P1 : constant := 0.75753_18015_94227_76666E-2; 661 P2 : constant := 0.31555_19276_56846_46356E-4; 662 663 Q0 : constant := 0.5; 664 Q1 : constant := 0.56817_30269_85512_21787E-1; 665 Q2 : constant := 0.63121_89437_43985_02557E-3; 666 Q3 : constant := 0.75104_02839_98700_46114E-6; 667 668 C1 : constant := 8#0.543#; 669 C2 : constant := -2.1219_44400_54690_58277E-4; 670 Le : constant := 1.4426_95040_88896_34074; 671 672 XN : Float_Type'Base; 673 P, Q, R : Float_Type'Base; 674 675 begin 676 if X = 0.0 then 677 return 1.0; 678 end if; 679 680 XN := Float_Type'Base'Rounding (X * Le); 681 G := (X - XN * C1) - XN * C2; 682 Z := G * G; 683 P := G * ((P2 * Z + P1) * Z + P0); 684 Q := ((Q3 * Z + Q2) * Z + Q1) * Z + Q0; 685 R := 0.5 + P / (Q - P); 686 687 R := Float_Type'Base'Scaling (R, Integer (XN) + 1); 688 689 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows 690 -- is False, then we can just leave it as an infinity (and indeed we 691 -- prefer to do so). But if Machine_Overflows is True, then we have to 692 -- raise a Constraint_Error exception as required by the RM. 693 694 if Float_Type'Machine_Overflows and then not R'Valid then 695 raise Constraint_Error; 696 else 697 return R; 698 end if; 699 700 end Exp_Strict; 701 702 ---------------- 703 -- Local_Atan -- 704 ---------------- 705 706 function Local_Atan 707 (Y : Float_Type'Base; 708 X : Float_Type'Base := 1.0) return Float_Type'Base 709 is 710 Z : Float_Type'Base; 711 Raw_Atan : Float_Type'Base; 712 713 begin 714 Z := (if abs Y > abs X then abs (X / Y) else abs (Y / X)); 715 716 Raw_Atan := 717 (if Z < Sqrt_Epsilon then Z 718 elsif Z = 1.0 then Pi / 4.0 719 else Float_Type'Base (Aux.Atan (Double (Z)))); 720 721 if abs Y > abs X then 722 Raw_Atan := Half_Pi - Raw_Atan; 723 end if; 724 725 if X > 0.0 then 726 return Float_Type'Copy_Sign (Raw_Atan, Y); 727 else 728 return Float_Type'Copy_Sign (Pi - Raw_Atan, Y); 729 end if; 730 end Local_Atan; 731 732 --------- 733 -- Log -- 734 --------- 735 736 -- Natural base 737 738 function Log (X : Float_Type'Base) return Float_Type'Base is 739 begin 740 if X < 0.0 then 741 raise Argument_Error; 742 743 elsif X = 0.0 then 744 raise Constraint_Error; 745 746 elsif X = 1.0 then 747 return 0.0; 748 end if; 749 750 return Float_Type'Base (Aux.Log (Double (X))); 751 end Log; 752 753 -- Arbitrary base 754 755 function Log (X, Base : Float_Type'Base) return Float_Type'Base is 756 begin 757 if X < 0.0 then 758 raise Argument_Error; 759 760 elsif Base <= 0.0 or else Base = 1.0 then 761 raise Argument_Error; 762 763 elsif X = 0.0 then 764 raise Constraint_Error; 765 766 elsif X = 1.0 then 767 return 0.0; 768 end if; 769 770 return Float_Type'Base (Aux.Log (Double (X)) / Aux.Log (Double (Base))); 771 end Log; 772 773 --------- 774 -- Sin -- 775 --------- 776 777 -- Natural cycle 778 779 function Sin (X : Float_Type'Base) return Float_Type'Base is 780 begin 781 if abs X < Sqrt_Epsilon then 782 return X; 783 end if; 784 785 return Float_Type'Base (Aux.Sin (Double (X))); 786 end Sin; 787 788 -- Arbitrary cycle 789 790 function Sin (X, Cycle : Float_Type'Base) return Float_Type'Base is 791 T : Float_Type'Base; 792 793 begin 794 if Cycle <= 0.0 then 795 raise Argument_Error; 796 797 -- If X is zero, return it as the result, preserving the argument sign. 798 -- Is this test really needed on any machine ??? 799 800 elsif X = 0.0 then 801 return X; 802 end if; 803 804 T := Float_Type'Base'Remainder (X, Cycle); 805 806 -- The following two reductions reduce the argument to the interval 807 -- [-0.25 * Cycle, 0.25 * Cycle]. This reduction is exact and is needed 808 -- to prevent inaccuracy that may result if the sine function uses a 809 -- different (more accurate) value of Pi in its reduction than is used 810 -- in the multiplication with Two_Pi. 811 812 if abs T > 0.25 * Cycle then 813 T := 0.5 * Float_Type'Copy_Sign (Cycle, T) - T; 814 end if; 815 816 -- Could test for 12.0 * abs T = Cycle, and return an exact value in 817 -- those cases. It is not clear this is worth the extra test though. 818 819 return Float_Type'Base (Aux.Sin (Double (T / Cycle * Two_Pi))); 820 end Sin; 821 822 ---------- 823 -- Sinh -- 824 ---------- 825 826 function Sinh (X : Float_Type'Base) return Float_Type'Base is 827 Lnv : constant Float_Type'Base := 8#0.542714#; 828 V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4; 829 Y : constant Float_Type'Base := abs X; 830 F : constant Float_Type'Base := Y * Y; 831 Z : Float_Type'Base; 832 833 Float_Digits_1_6 : constant Boolean := Float_Type'Digits < 7; 834 835 begin 836 if Y < Sqrt_Epsilon then 837 return X; 838 839 elsif Y > Log_Inverse_Epsilon then 840 Z := Exp_Strict (Y - Lnv); 841 Z := Z + V2minus1 * Z; 842 843 elsif Y < 1.0 then 844 845 if Float_Digits_1_6 then 846 847 -- Use expansion provided by Cody and Waite, p. 226. Note that 848 -- leading term of the polynomial in Q is exactly 1.0. 849 850 declare 851 P0 : constant := -0.71379_3159E+1; 852 P1 : constant := -0.19033_3399E+0; 853 Q0 : constant := -0.42827_7109E+2; 854 855 begin 856 Z := Y + Y * F * (P1 * F + P0) / (F + Q0); 857 end; 858 859 else 860 declare 861 P0 : constant := -0.35181_28343_01771_17881E+6; 862 P1 : constant := -0.11563_52119_68517_68270E+5; 863 P2 : constant := -0.16375_79820_26307_51372E+3; 864 P3 : constant := -0.78966_12741_73570_99479E+0; 865 Q0 : constant := -0.21108_77005_81062_71242E+7; 866 Q1 : constant := 0.36162_72310_94218_36460E+5; 867 Q2 : constant := -0.27773_52311_96507_01667E+3; 868 869 begin 870 Z := Y + Y * F * (((P3 * F + P2) * F + P1) * F + P0) 871 / (((F + Q2) * F + Q1) * F + Q0); 872 end; 873 end if; 874 875 else 876 Z := Exp_Strict (Y); 877 Z := 0.5 * (Z - 1.0 / Z); 878 end if; 879 880 if X > 0.0 then 881 return Z; 882 else 883 return -Z; 884 end if; 885 end Sinh; 886 887 ---------- 888 -- Sqrt -- 889 ---------- 890 891 function Sqrt (X : Float_Type'Base) return Float_Type'Base is 892 begin 893 if X < 0.0 then 894 raise Argument_Error; 895 896 -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE 897 898 elsif X = 0.0 then 899 return X; 900 end if; 901 902 return Float_Type'Base (Aux.Sqrt (Double (X))); 903 end Sqrt; 904 905 --------- 906 -- Tan -- 907 --------- 908 909 -- Natural cycle 910 911 function Tan (X : Float_Type'Base) return Float_Type'Base is 912 begin 913 if abs X < Sqrt_Epsilon then 914 return X; 915 end if; 916 917 -- Note: if X is exactly pi/2, then we should raise an exception, since 918 -- the result would overflow. But for all floating-point formats we deal 919 -- with, it is impossible for X to be exactly pi/2, and the result is 920 -- always in range. 921 922 return Float_Type'Base (Aux.Tan (Double (X))); 923 end Tan; 924 925 -- Arbitrary cycle 926 927 function Tan (X, Cycle : Float_Type'Base) return Float_Type'Base is 928 T : Float_Type'Base; 929 930 begin 931 if Cycle <= 0.0 then 932 raise Argument_Error; 933 934 elsif X = 0.0 then 935 return X; 936 end if; 937 938 T := Float_Type'Base'Remainder (X, Cycle); 939 940 if abs T = 0.25 * Cycle then 941 raise Constraint_Error; 942 943 elsif abs T = 0.5 * Cycle then 944 return 0.0; 945 946 else 947 T := T / Cycle * Two_Pi; 948 return Sin (T) / Cos (T); 949 end if; 950 951 end Tan; 952 953 ---------- 954 -- Tanh -- 955 ---------- 956 957 function Tanh (X : Float_Type'Base) return Float_Type'Base is 958 P0 : constant Float_Type'Base := -0.16134_11902_39962_28053E+4; 959 P1 : constant Float_Type'Base := -0.99225_92967_22360_83313E+2; 960 P2 : constant Float_Type'Base := -0.96437_49277_72254_69787E+0; 961 962 Q0 : constant Float_Type'Base := 0.48402_35707_19886_88686E+4; 963 Q1 : constant Float_Type'Base := 0.22337_72071_89623_12926E+4; 964 Q2 : constant Float_Type'Base := 0.11274_47438_05349_49335E+3; 965 Q3 : constant Float_Type'Base := 0.10000_00000_00000_00000E+1; 966 967 Half_Ln3 : constant Float_Type'Base := 0.54930_61443_34054_84570; 968 969 P, Q, R : Float_Type'Base; 970 Y : constant Float_Type'Base := abs X; 971 G : constant Float_Type'Base := Y * Y; 972 973 Float_Type_Digits_15_Or_More : constant Boolean := 974 Float_Type'Digits > 14; 975 976 begin 977 if X < Half_Log_Epsilon then 978 return -1.0; 979 980 elsif X > -Half_Log_Epsilon then 981 return 1.0; 982 983 elsif Y < Sqrt_Epsilon then 984 return X; 985 986 elsif Y < Half_Ln3 987 and then Float_Type_Digits_15_Or_More 988 then 989 P := (P2 * G + P1) * G + P0; 990 Q := ((Q3 * G + Q2) * G + Q1) * G + Q0; 991 R := G * (P / Q); 992 return X + X * R; 993 994 else 995 return Float_Type'Base (Aux.Tanh (Double (X))); 996 end if; 997 end Tanh; 998 999end Ada.Numerics.Generic_Elementary_Functions; 1000