gcc/ada/a-ngelfu.adb

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