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