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