gcc/ada/a-ngelfu.adb

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