File : s-lilodo-ada.adb


   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT COMPILER COMPONENTS                         --
   4 --                                                                          --
   5 --                S Y S T E M . L I B M _ L O N G _ D O U B L E             --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --           Copyright (C) 2014-2015, 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 --                                                                          --
  19 --                                                                          --
  20 --                                                                          --
  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 is the Ada Cert Math specific version of s-lilodo.adb
  33 
  34 --  When Cody and Waite implementation is cited, it refers to the
  35 --  Software Manual for the Elementary Functions by William J. Cody, Jr.
  36 --  and William Waite, published by Prentice-Hall Series in Computational
  37 --  Mathematics. Copyright 1980. ISBN 0-13-822064-6.
  38 
  39 --  When Hart implementation is cited, it refers to
  40 --  "Computer Approximations" by John F. Hart, published by Krieger.
  41 --  Copyright 1968, Reprinted 1978 w/ corrections. ISBN 0-88275-642-7.
  42 
  43 with Ada.Numerics; use Ada.Numerics;
  44 with System.Libm; use System.Libm;
  45 with System.Libm_Long_Double.Squareroot;
  46 
  47 package body System.Libm_Long_Double is
  48    subtype LLF is Long_Long_Float;
  49 
  50    pragma Assert (LLF'Machine_Radix = 2);
  51    pragma Assert (LLF'Machine_Mantissa = 64);
  52 
  53    LLF_HM : constant Integer := Long_Long_Float'Machine_Mantissa / 2;
  54 
  55    Sqrt_Epsilon_LLF : constant Long_Long_Float :=
  56                         Sqrt_2 ** (1 - Long_Long_Float'Machine_Mantissa);
  57 
  58    type Long_Long_Float_Table is array (Positive range <>) of Long_Long_Float;
  59 
  60    --  A1 (i) = Float (2**((1-i)/16))
  61 
  62    A1_Tab_LLF : constant Long_Long_Float_Table :=
  63                  (1.0,
  64                   Long_Long_Float'Machine (Root16_Half),
  65                   Long_Long_Float'Machine (Root16_Half**2),
  66                   Long_Long_Float'Machine (Root16_Half**3),
  67                   Long_Long_Float'Machine (Root16_Half**4),
  68                   Long_Long_Float'Machine (Root16_Half**5),
  69                   Long_Long_Float'Machine (Root16_Half**6),
  70                   Long_Long_Float'Machine (Root16_Half**7),
  71                   Long_Long_Float'Machine (Root16_Half**8),
  72                   Long_Long_Float'Machine (Root16_Half**9),
  73                   Long_Long_Float'Machine (Root16_Half**10),
  74                   Long_Long_Float'Machine (Root16_Half**11),
  75                   Long_Long_Float'Machine (Root16_Half**12),
  76                   Long_Long_Float'Machine (Root16_Half**13),
  77                   Long_Long_Float'Machine (Root16_Half**14),
  78                   Long_Long_Float'Machine (Root16_Half**15),
  79                   0.5);
  80 
  81    --  A2 (i) = 2**((1-2i)/16) - A1(2i)
  82 
  83    A2_Tab_LLF : constant Long_Long_Float_Table :=
  84                  (Root16_Half     - Long_Long_Float'Machine (Root16_Half),
  85                   Root16_Half**3  - Long_Long_Float'Machine (Root16_Half**3),
  86                   Root16_Half**5  - Long_Long_Float'Machine (Root16_Half**5),
  87                   Root16_Half**7  - Long_Long_Float'Machine (Root16_Half**7),
  88                   Root16_Half**9  - Long_Long_Float'Machine (Root16_Half**9),
  89                   Root16_Half**11 - Long_Long_Float'Machine (Root16_Half**11),
  90                   Root16_Half**13 - Long_Long_Float'Machine (Root16_Half**13),
  91                   Root16_Half**15 - Long_Long_Float'Machine (Root16_Half**15));
  92 
  93    --  Intermediary functions
  94 
  95    function Reconstruct_Pow
  96      (Z : Long_Long_Float;
  97       P : Integer;
  98       M : Integer) return Long_Long_Float;
  99 
 100    procedure Reduce_044
 101      (X         : Long_Long_Float;
 102       Reduced_X : out Long_Long_Float;
 103       P         : out Integer)
 104      with Post => abs (Reduced_X) < 0.044;
 105 
 106    function Reduce_1_16 (X : Long_Long_Float) return Long_Long_Float
 107      with Post => abs (X - Reduce_1_16'Result) <= 0.0625;
 108 
 109    function Reduce_1_16 (X : Long_Long_Float) return Long_Long_Float is
 110      (LLF'Machine_Rounding (X * 16.0) * (1.0 / 16.0));
 111 
 112    package Long_Long_Float_Approximations is
 113      new Generic_Approximations (Long_Long_Float, Mantissa => 64);
 114 
 115    use Long_Long_Float_Approximations;
 116 
 117    --  Local declarations
 118 
 119    procedure Reduce_Half_Pi_Large (X : in out LLF;  N : LLF; Q : out Quadrant);
 120 
 121    procedure Split_Veltkamp
 122      (X    : Long_Long_Float;
 123       X_Hi : out Long_Long_Float;
 124       X_Lo : out Long_Long_Float)
 125      with Post => X = X_Hi + X_Lo;
 126 
 127    function Multiply_Add (X, Y, Z : LLF) return LLF is (X * Y + Z);
 128 
 129    --  The following functions reduce a positive X into the range
 130    --  -ln (2) / 2 .. ln (2) / 2
 131 
 132    --  It returns a reduced X and an integer N such that:
 133    --  X'Old = X'New + N * Log (2)
 134 
 135    --  It is used by Exp function
 136 
 137    --  The result should be correctly rounded
 138 
 139    procedure Reduce_Ln_2 (X : in out Long_Long_Float; N : out Integer)
 140      with Pre  => abs (X) <= Long_Long_Float'Ceiling
 141                     (Long_Long_Float'Pred (11356.52340_62941_439) * Inv_Ln_2);
 142    --  @llr Reduce_Ln_2 Long_Long_Float
 143    --  The following is postcondition doesn't hold. Suspicious "=" ???
 144    --  Post => abs (X) <= Ln_2 / 2.0 and
 145    --          X'Old = X + Long_Long_Float (N) * Ln_2;
 146 
 147    --  The reduction is used by the Sin, Cos and Tan functions.
 148 
 149    procedure Reduce_Half_Pi (X : in out Long_Long_Float; Q : out Quadrant)
 150      with Pre  => X >= 0.0,
 151           Post => abs (X) <= Max_Red_Trig_Arg;
 152    --  @llr Reduce_Half_Pi Long_Long_Float
 153    --  The following functions reduce a positive X into the range
 154    --  -(Pi/4 + E) .. Pi/4 + E, with E a small fraction of Pi.
 155    --
 156    --  The reason the normalization is not strict is that the computation of
 157    --  the number of times to subtract half Pi is not exact. The rounding
 158    --  error is worst for large arguments, where the number of bits behind
 159    --  the radix point reduces to half the mantissa bits.
 160 
 161    --  While it would be possible to correct for this, the polynomial
 162    --  approximations work well for values slightly outside the -Pi/4 .. Pi/4
 163    --  interval, so it is easier for both error analysis and implementation
 164    --  to leave the reduction non-strict, and assume the reduced argument is
 165    --  within -0.26 * Pi .. 0.26 * Pi rather than a quarter of pi.
 166 
 167    --  The reduction is guaranteed to be correct to within 0.501 ulp for
 168    --  values of X for which Ada's accuracy guarantees apply:
 169    --     abs X <= 2.0**(T'Machine_Mantissa / 2)
 170    --  For values outside this range, an attempt is made to have significance
 171    --  decrease only proportionally with increase of magnitued. In any case,
 172    --  for all finite arguments, the reduction will succeed, though the reduced
 173    --  value may not agree with the mathematically correct value in even its
 174    --  sign.
 175 
 176    ---------------------
 177    -- Reconstruct_Pow --
 178    ---------------------
 179 
 180    function Reconstruct_Pow
 181      (Z : Long_Long_Float;
 182       P : Integer;
 183       M : Integer) return Long_Long_Float
 184    is
 185       --  Cody and Waite implementation (in "**" function page 84)
 186 
 187       --  The following computation is carried out in two steps. First add 1 to
 188       --  Z and multiply by 2**(-P/16). Then multiply the result by 2**M.
 189 
 190       Result : Long_Long_Float;
 191 
 192    begin
 193       Result := A1_Tab_LLF (P + 1) * Z;
 194       return Long_Long_Float'Scaling (Result, M);
 195    end Reconstruct_Pow;
 196 
 197    ----------------
 198    -- Reduce_044 --
 199    ----------------
 200 
 201    procedure Reduce_044
 202      (X         : Long_Long_Float;
 203       Reduced_X : out Long_Long_Float;
 204       P         : out Integer)
 205    is
 206       --  Cody and Waite implementation (in "**" function page 84)
 207       --  The output is:
 208 
 209       --  P is the biggest odd Integer in range 1 .. 15 such that
 210       --    2^((1-P)/16) <= X.
 211 
 212       --  Reduced_X equals 2 * (X-2^(-P/16)) / (X + 2^(-P/16)).
 213       --  abs (Reduced_X) <= max (2^(2-P/16)-2^(1-P/16)) <= 0.443.
 214 
 215    begin
 216       P := 1;
 217 
 218       if X <= A1_Tab_LLF (9) then
 219          P := 9;
 220       end if;
 221 
 222       if X <= A1_Tab_LLF (P + 4) then
 223          P := P + 4;
 224       end if;
 225 
 226       if X <= A1_Tab_LLF (P + 2) then
 227          P := P + 2;
 228       end if;
 229 
 230       Reduced_X := (X - A1_Tab_LLF (P + 1)) - A2_Tab_LLF ((P + 1) / 2);
 231       Reduced_X := Reduced_X / (X + A1_Tab_LLF (P + 1));
 232       Reduced_X := Reduced_X + Reduced_X;
 233    end Reduce_044;
 234 
 235    --------------------
 236    -- Instantiations --
 237    --------------------
 238 
 239    package Instantiations is
 240       function Acos  is new Generic_Acos (LLF);
 241       function Atan2 is new Generic_Atan2 (LLF);
 242    end Instantiations;
 243 
 244    --------------------
 245    -- Split_Veltkamp --
 246    --------------------
 247 
 248    procedure Split_Veltkamp
 249      (X    : Long_Long_Float;
 250       X_Hi : out Long_Long_Float;
 251       X_Lo : out Long_Long_Float)
 252    is
 253       M : constant LLF := 0.5 + 2.0**(1 - LLF'Machine_Mantissa / 2);
 254    begin
 255       X_Hi := X * M - (X * M - X);
 256       X_Lo := X - X_Hi;
 257    end Split_Veltkamp;
 258 
 259    -----------------
 260    -- Reduce_Ln_2 --
 261    -----------------
 262 
 263    procedure Reduce_Ln_2 (X : in out Long_Long_Float; N : out Integer) is
 264       L1 : constant := Long_Long_Float'Leading_Part (Ln_2, LLF_HM);
 265       L2 : constant := Long_Long_Float'Leading_Part (Ln_2 - L1, LLF_HM);
 266       L3 : constant := Ln_2 - L2 - L1;
 267       XN : constant Long_Long_Float := Long_Long_Float'Rounding (X * Inv_Ln_2);
 268 
 269    begin
 270       --  The argument passed to the function is smaller than Ymax * 1/log(2)
 271       --  No overflow is possible for N (Ymax is the largest machine number
 272       --  less than Log (LLF'Last)).
 273 
 274       N := Integer (XN);
 275       X := ((X - XN * L1) - XN * L2) - XN * L3;
 276 
 277       if X < -Ln_2 / 2.0 then
 278          X := X + Ln_2;
 279          N := N - 1;
 280       end if;
 281 
 282       if X > Ln_2 / 2.0 then
 283          X := X - Ln_2;
 284          N := N + 1;
 285       end if;
 286    end Reduce_Ln_2;
 287 
 288    --------------------
 289    -- Reduce_Half_Pi --
 290    --------------------
 291 
 292    procedure Reduce_Half_Pi (X : in out Long_Long_Float; Q : out Quadrant) is
 293       K      : constant := Pi / 2.0;
 294       Bits_N : constant := 3;
 295       Max_N  : constant := 2.0**Bits_N - 1.0;
 296       Max_X  : constant LLF := LLF'Pred (K * Max_N); -- About 3.5 * Pi
 297 
 298       Bits_C : constant := LLF'Machine_Mantissa - Bits_N;
 299       C1     : constant LLF := LLF'Leading_Part (K, Bits_C);
 300       C2     : constant LLF := K - C1;
 301 
 302       N : constant LLF := LLF'Machine_Rounding (X * K**(-1));
 303 
 304    begin
 305       if not X'Valid then
 306          X := X - X;
 307          Q := 0;
 308 
 309       elsif abs X > Max_X then
 310          Reduce_Half_Pi_Large (X, N, Q);
 311 
 312       else
 313          pragma Assert (if X'Valid then abs N <= Max_N);
 314          X := (X - N * C1) - N * C2;
 315          Q := Integer (N) mod 4;
 316       end if;
 317    end Reduce_Half_Pi;
 318 
 319    --------------------------
 320    -- Reduce_Half_Pi_Large --
 321    --------------------------
 322 
 323    procedure Reduce_Half_Pi_Large
 324      (X : in out LLF;
 325       N : LLF;
 326       Q : out Quadrant)
 327    is
 328       type Int_64 is range -2**63 .. 2**63 - 1; -- used for conversions
 329 
 330       HM : constant Positive := LLF'Machine_Mantissa / 2;
 331       C1 : constant LLF := LLF'Leading_Part (Half_Pi, HM);
 332       C2 : constant LLF := LLF'Leading_Part (Half_Pi - C1, HM);
 333       C3 : constant LLF := LLF'Leading_Part (Half_Pi - C1 - C2, HM);
 334       C4 : constant LLF := Half_Pi - C1 - C2 - C3;
 335 
 336       K    : LLF := N;
 337       K_Hi : LLF;
 338       K_Lo : LLF;
 339 
 340    begin
 341       Q := 0;
 342 
 343       loop
 344          Split_Veltkamp (X => K, X_Hi => K_Hi, X_Lo => K_Lo);
 345 
 346          X := Multiply_Add (-K_Hi, C1, X);
 347          X := Multiply_Add (-K_Hi, C2, Multiply_Add (-K_Lo, C1, X));
 348          X := Multiply_Add (-K_Hi, C3, Multiply_Add (-K_Lo, C2, X));
 349          X := Multiply_Add (-K_Hi, C4, Multiply_Add (-K_Lo, C3, X));
 350          X := Multiply_Add (-K_Lo, C4, X);
 351 
 352          if abs K < 2.0**62 then
 353             Q := Quadrant ((Int_64 (Q) + Int_64 (N)) mod 4);
 354 
 355          elsif abs K_Lo <= 2.0**62 then
 356             Q := Quadrant ((Int_64 (Q) + Int_64 (N)) mod 4);
 357          end if;
 358 
 359          exit when X in -0.26 * Pi .. 0.26 * Pi;
 360 
 361          K := LLF'Machine_Rounding (X * Half_Pi**(-1));
 362       end loop;
 363    end Reduce_Half_Pi_Large;
 364 
 365    ----------
 366    -- Acos --
 367    ----------
 368 
 369    function Acos (X : LLF) return LLF is (Instantiations.Acos (X));
 370 
 371    -----------
 372    -- Acosh --
 373    -----------
 374 
 375    function Acosh (X : LLF) return LLF is
 376 
 377       --  Math based implementation using Log1p: x-> Log (1+x)
 378 
 379       T : constant LLF := X - 1.0;
 380 
 381    begin
 382       if X > 1.0 / Sqrt_Epsilon_LLF then
 383          return Log (X) + Ln_2;
 384       elsif X < 2.0 then
 385          return Log1p (T + Sqrt (2.0 * T + T * T));
 386       else
 387          return Log (X + Sqrt ((X - 1.0) * (X + 1.0)));
 388       end if;
 389    end Acosh;
 390 
 391    ----------
 392    -- Asin --
 393    ----------
 394 
 395    function Asin (X : LLF) return LLF is
 396      (Long_Long_Float_Approximations.Asin (X));
 397 
 398    -----------
 399    -- Asinh --
 400    -----------
 401 
 402    function Asinh (X : LLF) return LLF is
 403 
 404       --  Math based implementation using Log1p: x-> Log (1+x)
 405 
 406       Y   : constant LLF := abs X;
 407       G   : constant LLF := X * X;
 408       Res : LLF;
 409 
 410    begin
 411       if Y < Sqrt_Epsilon_LLF then
 412          Res := Y;
 413       elsif Y > 1.0 / Sqrt_Epsilon_LLF then
 414          Res := Log (Y) + Ln_2;
 415       elsif Y < 2.0 then
 416          Res := Log1p (Y + G / (1.0 + Sqrt (G + 1.0)));
 417       else
 418          Res := Log (Y + Sqrt (G + 1.0));
 419       end if;
 420 
 421       return LLF'Copy_Sign (Res, X);
 422    end Asinh;
 423 
 424    ----------
 425    -- Atan --
 426    ----------
 427 
 428    function Atan (X : LLF) return LLF is (Instantiations.Atan2 (X, 1.0));
 429 
 430    -----------
 431    -- Atan2 --
 432    -----------
 433 
 434    function Atan2 (Y, X : LLF) return LLF is (Instantiations.Atan2 (Y, X));
 435 
 436    -----------
 437    -- Atanh --
 438    -----------
 439 
 440    function Atanh (X : LLF) return LLF is
 441 
 442       --  Math based implementation using Log1p: x-> Log (1+x)
 443 
 444      (if X >= 0.0
 445       then  Log1p (2.0 * X / (1.0 - X)) / 2.0
 446       else -Log1p (-2.0 * X / (1.0 + X)) / 2.0);
 447 
 448    ---------
 449    -- Cos --
 450    ---------
 451 
 452    function Cos (X : LLF) return LLF is
 453 
 454       --  Math based implementation using Hart constants
 455 
 456       Q      : Quadrant;
 457       Result : LLF;
 458       Y      : LLF := abs (X);
 459 
 460    begin
 461       Reduce_Half_Pi (Y, Q);
 462 
 463       if Q mod 2 = 0 then
 464          Result := Approx_Cos (Y);
 465       else
 466          Result := Approx_Sin (Y);
 467       end if;
 468 
 469       return (if Q = 1 or else Q = 2 then -Result else Result);
 470    end Cos;
 471 
 472    ----------
 473    -- Cosh --
 474    ----------
 475 
 476    function Cosh (X : LLF) return LLF is
 477 
 478       --  Cody and Waite implementation (page 217)
 479 
 480       Y : constant LLF := abs (X);
 481 
 482       --  Because the overflow threshold for cosh(X) is beyond the overflow
 483       --  threshold for exp(X), it appears natural to reformulate the
 484       --  computation as:
 485 
 486       --    Cosh (X) = Exp (X - Log (2))
 487 
 488       --  But because Log (2) is not an exact machine number, the finite word
 489       --  length of the machine implies that the absolute error in X - Log (2),
 490       --  hence the transmitted error in Cosh (X) is proportional to the
 491       --  magnitude of X even when X is error-free.
 492 
 493       --  To avoid this problem, we revise the computation to
 494 
 495       --    Cosh (X) = V/2 * exp(X - Log (V))
 496 
 497       --  where Log (V) is an exact machine number slightly larger than Log (2)
 498       --  with the last few digits of its significand zero.
 499 
 500       Ln_V : constant := 8#0.542714#;
 501       --  Machine value slightly above Ln_2
 502 
 503       V_2 : constant := 0.24999_30850_04514_99336;
 504       --  V**(-2)
 505 
 506       V_2_1 : constant := 0.13830_27787_96019_02638E-4;
 507       --  V / 2 - 1
 508 
 509       --  Fix me ???
 510       Y_Bar : constant Long_Long_Float := 11356.52340_62941_439;
 511       --  Y_Bar is the last floating point for which exp (Y) does not overflow
 512       --  and exp (-Y) does not underflow
 513 
 514       W : LLF;
 515       Z : LLF;
 516 
 517    begin
 518       if Y >= Y_Bar then
 519          W := Y - Ln_V;
 520          Z := Exp (W);
 521          Z := Z + V_2 / Z;
 522 
 523          return Z + V_2_1 * Z; --  rewriting of V/2 * Z
 524 
 525       else
 526          Z := Exp (Y);
 527          return (Z + 1.0 / Z) / 2.0;
 528       end if;
 529    end Cosh;
 530 
 531    ---------
 532    -- Exp --
 533    ---------
 534 
 535    function Exp (X : LLF) return LLF is
 536 
 537       --  Cody and Waite implementation (page 60)
 538 
 539       N : Integer;
 540       R : LLF;
 541       Y : LLF := X;
 542 
 543       --  Fix me ???
 544       Ymax : constant LLF := LLF'Pred (11356.52340_62941_439);
 545       --  The largest machine number less than Log (LLF'Last)
 546 
 547    begin
 548       if abs (Y) < 2.0**(-LLF'Machine_Mantissa - 1) then
 549          return 1.0;
 550       end if;
 551 
 552       if abs Y > Ymax then
 553          return (if Y > 0.0 then Infinity else 0.0);
 554       end if;
 555 
 556       Reduce_Ln_2 (Y, N);
 557       R := Approx_Exp (Y);
 558       return Long_Long_Float'Scaling (R, N);
 559    end Exp;
 560 
 561    ----------
 562    -- Exp2 --
 563    ----------
 564 
 565    function Exp2 (X : LLF) return LLF is
 566 
 567       --  Implementation based on Cody and Waite Exp implementation (page 217)
 568       --  but using Hart constants
 569 
 570       N : Integer;
 571       Y : LLF := X;
 572       R : LLF;
 573 
 574       Result : LLF;
 575 
 576    begin
 577       if abs Y < 2.0**(-LLF'Machine_Mantissa - 1) then
 578          return 1.0;
 579       end if;
 580 
 581       if abs Y > LLF'Pred (LLF (LLF'Machine_Emax)) then
 582          return (if Y > 0.0 then Infinity else 0.0);
 583       end if;
 584 
 585       --  If X > Log(LLF'Emax) ???
 586 
 587       N := Integer (X);
 588       Y := Y - Long_Long_Float (N);
 589       R := Approx_Exp2 (Y);
 590 
 591       Result := Long_Long_Float'Scaling (R, N + 1);
 592 
 593       if Result /= Result then
 594          Result := (if X < LLF'First then 0.0 else Infinity);
 595       end if;
 596 
 597       return Result;
 598    end Exp2;
 599 
 600    ---------
 601    -- Log --
 602    ---------
 603 
 604    function Log (X : LLF) return LLF is
 605 
 606       --  Cody and Waite implementation (page 35)
 607 
 608       Exponent_X  : constant Integer := LLF'Exponent (X);
 609       XN          : LLF               := LLF (Exponent_X);
 610       Mantissa_X  : LLF               := LLF'Scaling (X, -Exponent_X);
 611       HM          : constant Integer := LLF'Machine_Mantissa / 2;
 612       L1          : constant LLF      := LLF'Leading_Part (Ln_2, HM);
 613       L2          : constant LLF      := Ln_2 - L1;
 614       Result      : LLF;
 615 
 616    begin
 617       if X <= 0.0 then
 618          if X < 0.0 then
 619             return NaN;
 620          else
 621             return -Infinity;
 622          end if;
 623 
 624          --  Making sure X is in Sqrt (0.5) .. Sqrt (2)
 625 
 626       elsif X > Long_Long_Float'Last then
 627          return X;
 628 
 629       elsif Mantissa_X <= Sqrt_Half then
 630          XN := XN - 1.0;
 631          Mantissa_X := Mantissa_X * 2.0;
 632       end if;
 633 
 634       Result := Approx_Log (Mantissa_X);
 635       Result := (XN * L2 + Result) + XN * L1;
 636 
 637       return Result;
 638    end Log;
 639 
 640    -----------
 641    -- Log1p --
 642    -----------
 643 
 644    function Log1p (X : LLF) return LLF is
 645 
 646       --  Quick implementation of Log1p not accurate to the Ada regular Log
 647       --  requirements, but accurate enough to compute inverse hyperbolic
 648       --  functions.
 649 
 650    begin
 651       if 1.0 + X = 1.0 then
 652          return X;
 653       elsif X > LLF'Last then
 654          return X;
 655       else
 656          return Log (1.0 + X) * (X / ((1.0 + X) - 1.0));
 657       end if;
 658    end Log1p;
 659 
 660    ----------
 661    -- Log2 --
 662    ----------
 663 
 664    function Log2 (X : LLF) return LLF is
 665 
 666       --  Quick implementation of Log2 not accurate to the Ada regular Log
 667       --  (base e) requirement on the whole definition interval but accurate
 668       --  enough on 0 .. 2**(-1/16).
 669 
 670      (Log (X) * (1.0 / Ln_2));
 671 
 672    ---------
 673    -- Pow --
 674    ---------
 675 
 676    function Pow (Left, Right : LLF) return LLF is
 677 
 678       --  Cody and Waite implementation (page 84)
 679 
 680       One_Over_Sixteen : constant := 0.0625;
 681 
 682       M : constant Integer := LLF'Exponent (Left);
 683       G : constant LLF      := LLF'Fraction (Left);
 684       Y : constant LLF      := Right;
 685       Z : LLF;
 686       P : Integer;
 687 
 688       U2, U1, Y1, Y2, W1, W2, W : LLF;
 689       MM, PP, IW1, I            : Integer;
 690       Is_Special                : Boolean;
 691       Special_Result            : LLF;
 692 
 693       procedure Pow_Special_Cases is new Generic_Pow_Special_Cases (LLF);
 694 
 695    begin
 696       --  Special values
 697 
 698       Pow_Special_Cases (Left, Right, Is_Special, Special_Result);
 699 
 700       if Is_Special then
 701          return Special_Result;
 702 
 703       else
 704          --  Left**Right is calculated using the formula
 705          --    2**(Right * Log2 (Left))
 706 
 707          Reduce_044 (G, Z, P);
 708 
 709          --  At this point, Z <= 0.044
 710 
 711          U2 := Approx_Power_Log (Z);
 712          U1 := LLF (M * 16 - P) * 0.0625; --  U2 + U1 = Log2 (Left)
 713 
 714          --  Forming the pseudo extended precision product of U * Right
 715 
 716          Y1  := Reduce_1_16 (Y);
 717          Y2  := Y - Y1;
 718 
 719          W   := U2 * Y + U1 * Y2;
 720          W1  := Reduce_1_16 (W);
 721          W2  := W - W1;
 722 
 723          W   := W1 + U1 * Y1;
 724          W1  := Reduce_1_16 (W);
 725          W2  := W2 + (W - W1);
 726 
 727          W   := Reduce_1_16 (W2);
 728          IW1 := Integer (16.0 * (W1 + W));
 729          W2  := W2 - W;
 730 
 731          if W2 > 0.0 then
 732             W2 := W2 - One_Over_Sixteen;
 733             IW1 := 1 + IW1;
 734          end if;
 735 
 736          if IW1 < 0 then
 737             I := 0;
 738          else
 739             I := 1;
 740          end if;
 741 
 742          MM := Integer (IW1 / 16) + I;
 743          PP := 16 * MM - IW1;
 744 
 745          Z := Approx_Exp2 (W2);
 746          return Reconstruct_Pow (Z, PP, MM);
 747       end if;
 748    end Pow;
 749 
 750    ---------
 751    -- Sin --
 752    ---------
 753 
 754    function Sin (X : LLF) return LLF is
 755 
 756       --  Math based implementation using Hart constants
 757 
 758       Y      : LLF := abs X;
 759       Q      : Quadrant;
 760       Result : LLF;
 761 
 762    begin
 763       Reduce_Half_Pi (Y, Q);
 764 
 765       if Q mod 2 = 0 then
 766          Result := Approx_Sin (Y);
 767       else
 768          Result := Approx_Cos (Y);
 769       end if;
 770 
 771       return LLF'Copy_Sign (1.0, X) * (if Q >= 2 then -Result else Result);
 772    end Sin;
 773 
 774    ----------
 775    -- Sinh --
 776    ----------
 777 
 778    function Sinh (X : LLF) return LLF is
 779 
 780       --  Cody and Waite implementation (page 217)
 781 
 782       Sign : constant LLF := LLF'Copy_Sign (1.0, X);
 783       Y    : constant LLF := abs X;
 784 
 785       --  Because the overflow threshold for sinh(X) is beyond the overflow
 786       --  threshold for exp(X), it appears natural to reformulate the
 787       --  computation as:
 788 
 789       --    Sinh (X) = Exp (X - Log (2))
 790 
 791       --  But because Log (2) is not an exact machine number, the finite word
 792       --  length of the machine implies that the absolute error in X - Log (2),
 793       --  hence the transmitted error in Sinh (X) is proportional to the
 794       --  magnitude of X even when X is error-free. To avoid this problem, we
 795       --  revise the computation to:
 796 
 797       --    Sinh (X) = V/2 * exp(X - Log (V))
 798 
 799       --  where Log (V) is an exact machine number slightly larger than Log (2)
 800       --  with the last few digits of its significand zero.
 801 
 802       Ln_V : constant := 8#0.542714#;
 803       --  Machine value slightly above Ln_2
 804 
 805       V_2 : constant := 0.24999_30850_04514_99336;
 806       --  V**(-2)
 807 
 808       V_2_1 : constant := 0.13830_27787_96019_02638E-4;
 809       --  V / 2 - 1
 810 
 811       --  Fix me ???
 812       Y_Bar : constant Long_Long_Float := 11356.52340_62941_439;
 813       --  The last floating point for which exp (X) does not overflow and
 814       --  exp (-x) does not underflow
 815 
 816       W : LLF;
 817       Z : LLF;
 818 
 819    begin
 820       if Y <= 1.0 then
 821          return Approx_Sinh (X);
 822       end if;
 823 
 824       if Y >= Y_Bar then
 825          W := Y - Ln_V;
 826          Z := Exp (W);
 827          Z := Z - V_2 / Z;
 828 
 829          return Sign * (Z + V_2_1 * Z); --  rewriting of V/2 * Z
 830 
 831       else
 832          Z := Exp (Y);
 833          return Sign * ((Z - 1.0 / Z) / 2.0);
 834       end if;
 835    end Sinh;
 836 
 837    ----------
 838    -- Sqrt --
 839    ----------
 840 
 841    function Sqrt (X : Long_Long_Float) return Long_Long_Float renames
 842      System.Libm_Long_Double.Squareroot.Sqrt;
 843 
 844    ---------
 845    -- Tan --
 846    ---------
 847 
 848    function Tan (X : LLF) return LLF is
 849 
 850       --  Math based implementation using Hart constants
 851 
 852       Y : LLF := abs X;
 853       N : Integer;
 854 
 855    begin
 856       if abs X < LLF'Last then
 857          Reduce_Half_Pi (Y, N);
 858       else
 859          return Infinity / Infinity;
 860       end if;
 861 
 862       --  The reconstruction is included in the algebraic fraction in
 863       --  Approx_Tan function.
 864 
 865       if N mod 2 = 0 then
 866          return Approx_Tan (Y) * LLF'Copy_Sign (1.0, X);
 867 
 868       else
 869          return Approx_Cot (Y) * LLF'Copy_Sign (1.0, X);
 870       end if;
 871    end Tan;
 872 
 873    ----------
 874    -- Tanh --
 875    ----------
 876 
 877    function Tanh (X : LLF) return LLF is
 878 
 879       --  Cody and Waite implementation (page 239)
 880 
 881       F      : constant LLF := abs (X);
 882       Xbig   : constant := Ln_2 * LLF (1 + LLF'Machine_Mantissa);
 883       LN_3_2 : constant := 0.54930_61443_34054_84570;
 884       Result : LLF;
 885 
 886    begin
 887       if F > Xbig then
 888          Result := 1.0;
 889       else
 890          if F > LN_3_2 then
 891             Result := 1.0 - 2.0 / (Exp (2.0 * F) + 1.0);
 892          else
 893             Result := Approx_Tanh (F);
 894          end if;
 895       end if;
 896 
 897       return LLF'Copy_Sign (Result, X);
 898    end Tanh;
 899 end System.Libm_Long_Double;