gcc/ada/a-ngelfu.adb

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