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