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