gcc/ada/libgnat/a-numaux__x86.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUN-TIME COMPONENTS                         --
   4 --                                                                          --
   5 --                     A D A . N U M E R I C S . A U X                      --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                        (Machine Version for x86)                         --
   9 --                                                                          --
  10 --          Copyright (C) 1998-2017, Free Software Foundation, Inc.         --
  11 --                                                                          --
  12 -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  13 -- terms of the  GNU General Public License as published  by the Free Soft- --
  14 -- ware  Foundation;  either version 3,  or (at your option) any later ver- --
  15 -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  16 -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  17 -- or FITNESS FOR A PARTICULAR PURPOSE.                                     --
  18 --                                                                          --
  19 -- As a special exception under Section 7 of GPL version 3, you are granted --
  20 -- additional permissions described in the GCC Runtime Library Exception,   --
  21 -- version 3.1, as published by the Free Software Foundation.               --
  22 --                                                                          --
  23 -- You should have received a copy of the GNU General Public License and    --
  24 -- a copy of the GCC Runtime Library Exception along with this program;     --
  25 -- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --
  26 -- <http://www.gnu.org/licenses/>.                                          --
  27 --                                                                          --
  28 -- GNAT was originally developed  by the GNAT team at  New York University. --
  29 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  30 --                                                                          --
  31 ------------------------------------------------------------------------------
  32
  33 with System.Machine_Code; use System.Machine_Code;
  34
  35 package body Ada.Numerics.Aux is
  36
  37    NL : constant String := ASCII.LF & ASCII.HT;
  38
  39    -----------------------
  40    -- Local subprograms --
  41    -----------------------
  42
  43    function Is_Nan (X : Double) return Boolean;
  44    --  Return True iff X is a IEEE NaN value
  45
  46    function Logarithmic_Pow (X, Y : Double) return Double;
  47    --  Implementation of X**Y using Exp and Log functions (binary base)
  48    --  to calculate the exponentiation. This is used by Pow for values
  49    --  for values of Y in the open interval (-0.25, 0.25)
  50
  51    procedure Reduce (X : in out Double; Q : out Natural);
  52    --  Implement reduction of X by Pi/2. Q is the quadrant of the final
  53    --  result in the range 0..3. The absolute value of X is at most Pi/4.
  54    --  It is needed to avoid a loss of accuracy for sin near Pi and cos
  55    --  near Pi/2 due to the use of an insufficiently precise value of Pi
  56    --  in the range reduction.
  57
  58    pragma Inline (Is_Nan);
  59    pragma Inline (Reduce);
  60
  61    --------------------------------
  62    -- Basic Elementary Functions --
  63    --------------------------------
  64
  65    --  This section implements a few elementary functions that are used to
  66    --  build the more complex ones. This ordering enables better inlining.
  67
  68    ----------
  69    -- Atan --
  70    ----------
  71
  72    function Atan (X : Double) return Double is
  73       Result  : Double;
  74
  75    begin
  76       Asm (Template =>
  77            "fld1" & NL
  78          & "fpatan",
  79          Outputs  => Double'Asm_Output ("=t", Result),
  80          Inputs   => Double'Asm_Input  ("0", X));
  81
  82       --  The result value is NaN iff input was invalid
  83
  84       if not (Result = Result) then
  85          raise Argument_Error;
  86       end if;
  87
  88       return Result;
  89    end Atan;
  90
  91    ---------
  92    -- Exp --
  93    ---------
  94
  95    function Exp (X : Double) return Double is
  96       Result : Double;
  97    begin
  98       Asm (Template =>
  99          "fldl2e               " & NL
 100        & "fmulp   %%st, %%st(1)" & NL -- X * log2 (E)
 101        & "fld     %%st(0)      " & NL
 102        & "frndint              " & NL -- Integer (X * Log2 (E))
 103        & "fsubr   %%st, %%st(1)" & NL -- Fraction (X * Log2 (E))
 104        & "fxch                 " & NL
 105        & "f2xm1                " & NL -- 2**(...) - 1
 106        & "fld1                 " & NL
 107        & "faddp   %%st, %%st(1)" & NL -- 2**(Fraction (X * Log2 (E)))
 108        & "fscale               " & NL -- E ** X
 109        & "fstp    %%st(1)      ",
 110          Outputs  => Double'Asm_Output ("=t", Result),
 111          Inputs   => Double'Asm_Input  ("0", X));
 112       return Result;
 113    end Exp;
 114
 115    ------------
 116    -- Is_Nan --
 117    ------------
 118
 119    function Is_Nan (X : Double) return Boolean is
 120    begin
 121       --  The IEEE NaN values are the only ones that do not equal themselves
 122
 123       return X /= X;
 124    end Is_Nan;
 125
 126    ---------
 127    -- Log --
 128    ---------
 129
 130    function Log (X : Double) return Double is
 131       Result : Double;
 132
 133    begin
 134       Asm (Template =>
 135          "fldln2               " & NL
 136        & "fxch                 " & NL
 137        & "fyl2x                " & NL,
 138          Outputs  => Double'Asm_Output ("=t", Result),
 139          Inputs   => Double'Asm_Input  ("0", X));
 140       return Result;
 141    end Log;
 142
 143    ------------
 144    -- Reduce --
 145    ------------
 146
 147    procedure Reduce (X : in out Double; Q : out Natural) is
 148       Half_Pi     : constant := Pi / 2.0;
 149       Two_Over_Pi : constant := 2.0 / Pi;
 150
 151       HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
 152       M  : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
 153       P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
 154       P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
 155       P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
 156       P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
 157       P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
 158                                                                  - P4, HM);
 159       P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
 160       K  : Double;
 161       R  : Integer;
 162
 163    begin
 164       --  For X < 2.0**HM, all products below are computed exactly.
 165       --  Due to cancellation effects all subtractions are exact as well.
 166       --  As no double extended floating-point number has more than 75
 167       --  zeros after the binary point, the result will be the correctly
 168       --  rounded result of X - K * (Pi / 2.0).
 169
 170       K := X * Two_Over_Pi;
 171       while abs K >= 2.0**HM loop
 172          K := K * M - (K * M - K);
 173          X :=
 174            (((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
 175          K := X * Two_Over_Pi;
 176       end loop;
 177
 178       --  If K is not a number (because X was not finite) raise exception
 179
 180       if Is_Nan (K) then
 181          raise Constraint_Error;
 182       end if;
 183
 184       --  Go through an integer temporary so as to use machine instructions
 185
 186       R := Integer (Double'Rounding (K));
 187       Q := R mod 4;
 188       K := Double (R);
 189       X := (((((X - K * P1) - K * P2) - K * P3) - K * P4) - K * P5) - K * P6;
 190    end Reduce;
 191
 192    ----------
 193    -- Sqrt --
 194    ----------
 195
 196    function Sqrt (X : Double) return Double is
 197       Result  : Double;
 198
 199    begin
 200       if X < 0.0 then
 201          raise Argument_Error;
 202       end if;
 203
 204       Asm (Template => "fsqrt",
 205            Outputs  => Double'Asm_Output ("=t", Result),
 206            Inputs   => Double'Asm_Input  ("0", X));
 207
 208       return Result;
 209    end Sqrt;
 210
 211    --------------------------------
 212    -- Other Elementary Functions --
 213    --------------------------------
 214
 215    --  These are built using the previously implemented basic functions
 216
 217    ----------
 218    -- Acos --
 219    ----------
 220
 221    function Acos (X : Double) return Double is
 222       Result  : Double;
 223
 224    begin
 225       Result := 2.0 * Atan (Sqrt ((1.0 - X) / (1.0 + X)));
 226
 227       --  The result value is NaN iff input was invalid
 228
 229       if Is_Nan (Result) then
 230          raise Argument_Error;
 231       end if;
 232
 233       return Result;
 234    end Acos;
 235
 236    ----------
 237    -- Asin --
 238    ----------
 239
 240    function Asin (X : Double) return Double is
 241       Result  : Double;
 242
 243    begin
 244       Result := Atan (X / Sqrt ((1.0 - X) * (1.0 + X)));
 245
 246       --  The result value is NaN iff input was invalid
 247
 248       if Is_Nan (Result) then
 249          raise Argument_Error;
 250       end if;
 251
 252       return Result;
 253    end Asin;
 254
 255    ---------
 256    -- Cos --
 257    ---------
 258
 259    function Cos (X : Double) return Double is
 260       Reduced_X : Double := abs X;
 261       Result    : Double;
 262       Quadrant  : Natural range 0 .. 3;
 263
 264    begin
 265       if Reduced_X > Pi / 4.0 then
 266          Reduce (Reduced_X, Quadrant);
 267
 268          case Quadrant is
 269             when 0 =>
 270                Asm (Template  => "fcos",
 271                   Outputs  => Double'Asm_Output ("=t", Result),
 272                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 273
 274             when 1 =>
 275                Asm (Template  => "fsin",
 276                   Outputs  => Double'Asm_Output ("=t", Result),
 277                   Inputs   => Double'Asm_Input  ("0", -Reduced_X));
 278
 279             when 2 =>
 280                Asm (Template  => "fcos ; fchs",
 281                   Outputs  => Double'Asm_Output ("=t", Result),
 282                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 283
 284             when 3 =>
 285                Asm (Template  => "fsin",
 286                   Outputs  => Double'Asm_Output ("=t", Result),
 287                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 288          end case;
 289
 290       else
 291          Asm (Template  => "fcos",
 292               Outputs  => Double'Asm_Output ("=t", Result),
 293               Inputs   => Double'Asm_Input  ("0", Reduced_X));
 294       end if;
 295
 296       return Result;
 297    end Cos;
 298
 299    ---------------------
 300    -- Logarithmic_Pow --
 301    ---------------------
 302
 303    function Logarithmic_Pow (X, Y : Double) return Double is
 304       Result  : Double;
 305    begin
 306       Asm (Template => ""             --  X                  : Y
 307        & "fyl2x                " & NL --  Y * Log2 (X)
 308        & "fld     %%st(0)      " & NL --  Y * Log2 (X)       : Y * Log2 (X)
 309        & "frndint              " & NL --  Int (...)          : Y * Log2 (X)
 310        & "fsubr   %%st, %%st(1)" & NL --  Int (...)          : Fract (...)
 311        & "fxch                 " & NL --  Fract (...)        : Int (...)
 312        & "f2xm1                " & NL --  2**Fract (...) - 1 : Int (...)
 313        & "fld1                 " & NL --  1 : 2**Fract (...) - 1 : Int (...)
 314        & "faddp   %%st, %%st(1)" & NL --  2**Fract (...)     : Int (...)
 315        & "fscale               ",     --  2**(Fract (...) + Int (...))
 316          Outputs  => Double'Asm_Output ("=t", Result),
 317          Inputs   =>
 318            (Double'Asm_Input  ("0", X),
 319             Double'Asm_Input  ("u", Y)));
 320       return Result;
 321    end Logarithmic_Pow;
 322
 323    ---------
 324    -- Pow --
 325    ---------
 326
 327    function Pow (X, Y : Double) return Double is
 328       type Mantissa_Type is mod 2**Double'Machine_Mantissa;
 329       --  Modular type that can hold all bits of the mantissa of Double
 330
 331       --  For negative exponents, do divide at the end of the processing
 332
 333       Negative_Y : constant Boolean := Y < 0.0;
 334       Abs_Y      : constant Double := abs Y;
 335
 336       --  During this function the following invariant is kept:
 337       --  X ** (abs Y) = Base**(Exp_High + Exp_Mid + Exp_Low) * Factor
 338
 339       Base : Double := X;
 340
 341       Exp_High : Double := Double'Floor (Abs_Y);
 342       Exp_Mid  : Double;
 343       Exp_Low  : Double;
 344       Exp_Int  : Mantissa_Type;
 345
 346       Factor : Double := 1.0;
 347
 348    begin
 349       --  Select algorithm for calculating Pow (integer cases fall through)
 350
 351       if Exp_High >= 2.0**Double'Machine_Mantissa then
 352
 353          --  In case of Y that is IEEE infinity, just raise constraint error
 354
 355          if Exp_High > Double'Safe_Last then
 356             raise Constraint_Error;
 357          end if;
 358
 359          --  Large values of Y are even integers and will stay integer
 360          --  after division by two.
 361
 362          loop
 363             --  Exp_Mid and Exp_Low are zero, so
 364             --    X**(abs Y) = Base ** Exp_High = (Base**2) ** (Exp_High / 2)
 365
 366             Exp_High := Exp_High / 2.0;
 367             Base := Base * Base;
 368             exit when Exp_High < 2.0**Double'Machine_Mantissa;
 369          end loop;
 370
 371       elsif Exp_High /= Abs_Y then
 372          Exp_Low := Abs_Y - Exp_High;
 373          Factor := 1.0;
 374
 375          if Exp_Low /= 0.0 then
 376
 377             --  Exp_Low now is in interval (0.0, 1.0)
 378             --  Exp_Mid := Double'Floor (Exp_Low * 4.0) / 4.0;
 379
 380             Exp_Mid := 0.0;
 381             Exp_Low := Exp_Low - Exp_Mid;
 382
 383             if Exp_Low >= 0.5 then
 384                Factor := Sqrt (X);
 385                Exp_Low := Exp_Low - 0.5;  -- exact
 386
 387                if Exp_Low >= 0.25 then
 388                   Factor := Factor * Sqrt (Factor);
 389                   Exp_Low := Exp_Low - 0.25; --  exact
 390                end if;
 391
 392             elsif Exp_Low >= 0.25 then
 393                Factor := Sqrt (Sqrt (X));
 394                Exp_Low := Exp_Low - 0.25; --  exact
 395             end if;
 396
 397             --  Exp_Low now is in interval (0.0, 0.25)
 398
 399             --  This means it is safe to call Logarithmic_Pow
 400             --  for the remaining part.
 401
 402             Factor := Factor * Logarithmic_Pow (X, Exp_Low);
 403          end if;
 404
 405       elsif X = 0.0 then
 406          return 0.0;
 407       end if;
 408
 409       --  Exp_High is non-zero integer smaller than 2**Double'Machine_Mantissa
 410
 411       Exp_Int := Mantissa_Type (Exp_High);
 412
 413       --  Standard way for processing integer powers > 0
 414
 415       while Exp_Int > 1 loop
 416          if (Exp_Int and 1) = 1 then
 417
 418             --  Base**Y = Base**(Exp_Int - 1) * Exp_Int for Exp_Int > 0
 419
 420             Factor := Factor * Base;
 421          end if;
 422
 423          --  Exp_Int is even and Exp_Int > 0, so
 424          --    Base**Y = (Base**2)**(Exp_Int / 2)
 425
 426          Base := Base * Base;
 427          Exp_Int := Exp_Int / 2;
 428       end loop;
 429
 430       --  Exp_Int = 1 or Exp_Int = 0
 431
 432       if Exp_Int = 1 then
 433          Factor := Base * Factor;
 434       end if;
 435
 436       if Negative_Y then
 437          Factor := 1.0 / Factor;
 438       end if;
 439
 440       return Factor;
 441    end Pow;
 442
 443    ---------
 444    -- Sin --
 445    ---------
 446
 447    function Sin (X : Double) return Double is
 448       Reduced_X : Double := X;
 449       Result    : Double;
 450       Quadrant  : Natural range 0 .. 3;
 451
 452    begin
 453       if abs X > Pi / 4.0 then
 454          Reduce (Reduced_X, Quadrant);
 455
 456          case Quadrant is
 457             when 0 =>
 458                Asm (Template  => "fsin",
 459                   Outputs  => Double'Asm_Output ("=t", Result),
 460                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 461
 462             when 1 =>
 463                Asm (Template  => "fcos",
 464                   Outputs  => Double'Asm_Output ("=t", Result),
 465                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 466
 467             when 2 =>
 468                Asm (Template  => "fsin",
 469                   Outputs  => Double'Asm_Output ("=t", Result),
 470                   Inputs   => Double'Asm_Input  ("0", -Reduced_X));
 471
 472             when 3 =>
 473                Asm (Template  => "fcos ; fchs",
 474                   Outputs  => Double'Asm_Output ("=t", Result),
 475                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 476          end case;
 477
 478       else
 479          Asm (Template  => "fsin",
 480             Outputs  => Double'Asm_Output ("=t", Result),
 481             Inputs   => Double'Asm_Input  ("0", Reduced_X));
 482       end if;
 483
 484       return Result;
 485    end Sin;
 486
 487    ---------
 488    -- Tan --
 489    ---------
 490
 491    function Tan (X : Double) return Double is
 492       Reduced_X : Double := X;
 493       Result    : Double;
 494       Quadrant  : Natural range 0 .. 3;
 495
 496    begin
 497       if abs X > Pi / 4.0 then
 498          Reduce (Reduced_X, Quadrant);
 499
 500          if Quadrant mod 2 = 0 then
 501             Asm (Template  => "fptan" & NL
 502                             & "ffree   %%st(0)"  & NL
 503                             & "fincstp",
 504                   Outputs  => Double'Asm_Output ("=t", Result),
 505                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 506          else
 507             Asm (Template  => "fsincos" & NL
 508                             & "fdivp   %%st, %%st(1)" & NL
 509                             & "fchs",
 510                   Outputs  => Double'Asm_Output ("=t", Result),
 511                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 512          end if;
 513
 514       else
 515          Asm (Template  =>
 516                "fptan                 " & NL
 517              & "ffree   %%st(0)      " & NL
 518              & "fincstp              ",
 519                Outputs  => Double'Asm_Output ("=t", Result),
 520                Inputs   => Double'Asm_Input  ("0", Reduced_X));
 521       end if;
 522
 523       return Result;
 524    end Tan;
 525
 526    ----------
 527    -- Sinh --
 528    ----------
 529
 530    function Sinh (X : Double) return Double is
 531    begin
 532       --  Mathematically Sinh (x) is defined to be (Exp (X) - Exp (-X)) / 2.0
 533
 534       if abs X < 25.0 then
 535          return (Exp (X) - Exp (-X)) / 2.0;
 536       else
 537          return Exp (X) / 2.0;
 538       end if;
 539    end Sinh;
 540
 541    ----------
 542    -- Cosh --
 543    ----------
 544
 545    function Cosh (X : Double) return Double is
 546    begin
 547       --  Mathematically Cosh (X) is defined to be (Exp (X) + Exp (-X)) / 2.0
 548
 549       if abs X < 22.0 then
 550          return (Exp (X) + Exp (-X)) / 2.0;
 551       else
 552          return Exp (X) / 2.0;
 553       end if;
 554    end Cosh;
 555
 556    ----------
 557    -- Tanh --
 558    ----------
 559
 560    function Tanh (X : Double) return Double is
 561    begin
 562       --  Return the Hyperbolic Tangent of x
 563
 564       --                                    x    -x
 565       --                                   e  - e        Sinh (X)
 566       --       Tanh (X) is defined to be -----------   = --------
 567       --                                    x    -x      Cosh (X)
 568       --                                   e  + e
 569
 570       if abs X > 23.0 then
 571          return Double'Copy_Sign (1.0, X);
 572       end if;
 573
 574       return 1.0 / (1.0 + Exp (-(2.0 * X))) - 1.0 / (1.0 + Exp (2.0 * X));
 575    end Tanh;
 576
 577 end Ada.Numerics.Aux;