gcc/ada/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-2016, 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    --  Implements 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.
  54
  55    pragma Inline (Is_Nan);
  56    pragma Inline (Reduce);
  57
  58    --------------------------------
  59    -- Basic Elementary Functions --
  60    --------------------------------
  61
  62    --  This section implements a few elementary functions that are used to
  63    --  build the more complex ones. This ordering enables better inlining.
  64
  65    ----------
  66    -- Atan --
  67    ----------
  68
  69    function Atan (X : Double) return Double is
  70       Result  : Double;
  71
  72    begin
  73       Asm (Template =>
  74            "fld1" & NL
  75          & "fpatan",
  76          Outputs  => Double'Asm_Output ("=t", Result),
  77          Inputs   => Double'Asm_Input  ("0", X));
  78
  79       --  The result value is NaN iff input was invalid
  80
  81       if not (Result = Result) then
  82          raise Argument_Error;
  83       end if;
  84
  85       return Result;
  86    end Atan;
  87
  88    ---------
  89    -- Exp --
  90    ---------
  91
  92    function Exp (X : Double) return Double is
  93       Result : Double;
  94    begin
  95       Asm (Template =>
  96          "fldl2e               " & NL
  97        & "fmulp   %%st, %%st(1)" & NL -- X * log2 (E)
  98        & "fld     %%st(0)      " & NL
  99        & "frndint              " & NL -- Integer (X * Log2 (E))
 100        & "fsubr   %%st, %%st(1)" & NL -- Fraction (X * Log2 (E))
 101        & "fxch                 " & NL
 102        & "f2xm1                " & NL -- 2**(...) - 1
 103        & "fld1                 " & NL
 104        & "faddp   %%st, %%st(1)" & NL -- 2**(Fraction (X * Log2 (E)))
 105        & "fscale               " & NL -- E ** X
 106        & "fstp    %%st(1)      ",
 107          Outputs  => Double'Asm_Output ("=t", Result),
 108          Inputs   => Double'Asm_Input  ("0", X));
 109       return Result;
 110    end Exp;
 111
 112    ------------
 113    -- Is_Nan --
 114    ------------
 115
 116    function Is_Nan (X : Double) return Boolean is
 117    begin
 118       --  The IEEE NaN values are the only ones that do not equal themselves
 119
 120       return not (X = X);
 121    end Is_Nan;
 122
 123    ---------
 124    -- Log --
 125    ---------
 126
 127    function Log (X : Double) return Double is
 128       Result : Double;
 129
 130    begin
 131       Asm (Template =>
 132          "fldln2               " & NL
 133        & "fxch                 " & NL
 134        & "fyl2x                " & NL,
 135          Outputs  => Double'Asm_Output ("=t", Result),
 136          Inputs   => Double'Asm_Input  ("0", X));
 137       return Result;
 138    end Log;
 139
 140    ------------
 141    -- Reduce --
 142    ------------
 143
 144    procedure Reduce (X : in out Double; Q : out Natural) is
 145       Half_Pi     : constant := Pi / 2.0;
 146       Two_Over_Pi : constant := 2.0 / Pi;
 147
 148       HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
 149       M  : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
 150       P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
 151       P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
 152       P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
 153       P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
 154       P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
 155                                                                  - P4, HM);
 156       P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
 157       K  : Double := X * Two_Over_Pi;
 158    begin
 159       --  For X < 2.0**32, all products below are computed exactly.
 160       --  Due to cancellation effects all subtractions are exact as well.
 161       --  As no double extended floating-point number has more than 75
 162       --  zeros after the binary point, the result will be the correctly
 163       --  rounded result of X - K * (Pi / 2.0).
 164
 165       while abs K >= 2.0**HM loop
 166          K := K * M - (K * M - K);
 167          X := (((((X - K * P1) - K * P2) - K * P3)
 168                      - K * P4) - K * P5) - K * P6;
 169          K := X * Two_Over_Pi;
 170       end loop;
 171
 172       if K /= K then
 173
 174          --  K is not a number, because X was not finite
 175
 176          raise Constraint_Error;
 177       end if;
 178
 179       K := Double'Rounding (K);
 180       Q := Integer (K) mod 4;
 181       X := (((((X - K * P1) - K * P2) - K * P3)
 182                   - K * P4) - K * P5) - K * P6;
 183    end Reduce;
 184
 185    ----------
 186    -- Sqrt --
 187    ----------
 188
 189    function Sqrt (X : Double) return Double is
 190       Result  : Double;
 191
 192    begin
 193       if X < 0.0 then
 194          raise Argument_Error;
 195       end if;
 196
 197       Asm (Template => "fsqrt",
 198            Outputs  => Double'Asm_Output ("=t", Result),
 199            Inputs   => Double'Asm_Input  ("0", X));
 200
 201       return Result;
 202    end Sqrt;
 203
 204    --------------------------------
 205    -- Other Elementary Functions --
 206    --------------------------------
 207
 208    --  These are built using the previously implemented basic functions
 209
 210    ----------
 211    -- Acos --
 212    ----------
 213
 214    function Acos (X : Double) return Double is
 215       Result  : Double;
 216
 217    begin
 218       Result := 2.0 * Atan (Sqrt ((1.0 - X) / (1.0 + X)));
 219
 220       --  The result value is NaN iff input was invalid
 221
 222       if Is_Nan (Result) then
 223          raise Argument_Error;
 224       end if;
 225
 226       return Result;
 227    end Acos;
 228
 229    ----------
 230    -- Asin --
 231    ----------
 232
 233    function Asin (X : Double) return Double is
 234       Result  : Double;
 235
 236    begin
 237       Result := Atan (X / Sqrt ((1.0 - X) * (1.0 + X)));
 238
 239       --  The result value is NaN iff input was invalid
 240
 241       if Is_Nan (Result) then
 242          raise Argument_Error;
 243       end if;
 244
 245       return Result;
 246    end Asin;
 247
 248    ---------
 249    -- Cos --
 250    ---------
 251
 252    function Cos (X : Double) return Double is
 253       Reduced_X : Double := abs X;
 254       Result    : Double;
 255       Quadrant  : Natural range 0 .. 3;
 256
 257    begin
 258       if Reduced_X > Pi / 4.0 then
 259          Reduce (Reduced_X, Quadrant);
 260
 261          case Quadrant is
 262             when 0 =>
 263                Asm (Template  => "fcos",
 264                   Outputs  => Double'Asm_Output ("=t", Result),
 265                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 266
 267             when 1 =>
 268                Asm (Template  => "fsin",
 269                   Outputs  => Double'Asm_Output ("=t", Result),
 270                   Inputs   => Double'Asm_Input  ("0", -Reduced_X));
 271
 272             when 2 =>
 273                Asm (Template  => "fcos ; fchs",
 274                   Outputs  => Double'Asm_Output ("=t", Result),
 275                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 276
 277             when 3 =>
 278                Asm (Template  => "fsin",
 279                   Outputs  => Double'Asm_Output ("=t", Result),
 280                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 281          end case;
 282
 283       else
 284          Asm (Template  => "fcos",
 285               Outputs  => Double'Asm_Output ("=t", Result),
 286               Inputs   => Double'Asm_Input  ("0", Reduced_X));
 287       end if;
 288
 289       return Result;
 290    end Cos;
 291
 292    ---------------------
 293    -- Logarithmic_Pow --
 294    ---------------------
 295
 296    function Logarithmic_Pow (X, Y : Double) return Double is
 297       Result  : Double;
 298    begin
 299       Asm (Template => ""             --  X                  : Y
 300        & "fyl2x                " & NL --  Y * Log2 (X)
 301        & "fld     %%st(0)      " & NL --  Y * Log2 (X)       : Y * Log2 (X)
 302        & "frndint              " & NL --  Int (...)          : Y * Log2 (X)
 303        & "fsubr   %%st, %%st(1)" & NL --  Int (...)          : Fract (...)
 304        & "fxch                 " & NL --  Fract (...)        : Int (...)
 305        & "f2xm1                " & NL --  2**Fract (...) - 1 : Int (...)
 306        & "fld1                 " & NL --  1 : 2**Fract (...) - 1 : Int (...)
 307        & "faddp   %%st, %%st(1)" & NL --  2**Fract (...)     : Int (...)
 308        & "fscale               ",     --  2**(Fract (...) + Int (...))
 309          Outputs  => Double'Asm_Output ("=t", Result),
 310          Inputs   =>
 311            (Double'Asm_Input  ("0", X),
 312             Double'Asm_Input  ("u", Y)));
 313       return Result;
 314    end Logarithmic_Pow;
 315
 316    ---------
 317    -- Pow --
 318    ---------
 319
 320    function Pow (X, Y : Double) return Double is
 321       type Mantissa_Type is mod 2**Double'Machine_Mantissa;
 322       --  Modular type that can hold all bits of the mantissa of Double
 323
 324       --  For negative exponents, do divide at the end of the processing
 325
 326       Negative_Y : constant Boolean := Y < 0.0;
 327       Abs_Y      : constant Double := abs Y;
 328
 329       --  During this function the following invariant is kept:
 330       --  X ** (abs Y) = Base**(Exp_High + Exp_Mid + Exp_Low) * Factor
 331
 332       Base : Double := X;
 333
 334       Exp_High : Double := Double'Floor (Abs_Y);
 335       Exp_Mid  : Double;
 336       Exp_Low  : Double;
 337       Exp_Int  : Mantissa_Type;
 338
 339       Factor : Double := 1.0;
 340
 341    begin
 342       --  Select algorithm for calculating Pow (integer cases fall through)
 343
 344       if Exp_High >= 2.0**Double'Machine_Mantissa then
 345
 346          --  In case of Y that is IEEE infinity, just raise constraint error
 347
 348          if Exp_High > Double'Safe_Last then
 349             raise Constraint_Error;
 350          end if;
 351
 352          --  Large values of Y are even integers and will stay integer
 353          --  after division by two.
 354
 355          loop
 356             --  Exp_Mid and Exp_Low are zero, so
 357             --    X**(abs Y) = Base ** Exp_High = (Base**2) ** (Exp_High / 2)
 358
 359             Exp_High := Exp_High / 2.0;
 360             Base := Base * Base;
 361             exit when Exp_High < 2.0**Double'Machine_Mantissa;
 362          end loop;
 363
 364       elsif Exp_High /= Abs_Y then
 365          Exp_Low := Abs_Y - Exp_High;
 366          Factor := 1.0;
 367
 368          if Exp_Low /= 0.0 then
 369
 370             --  Exp_Low now is in interval (0.0, 1.0)
 371             --  Exp_Mid := Double'Floor (Exp_Low * 4.0) / 4.0;
 372
 373             Exp_Mid := 0.0;
 374             Exp_Low := Exp_Low - Exp_Mid;
 375
 376             if Exp_Low >= 0.5 then
 377                Factor := Sqrt (X);
 378                Exp_Low := Exp_Low - 0.5;  -- exact
 379
 380                if Exp_Low >= 0.25 then
 381                   Factor := Factor * Sqrt (Factor);
 382                   Exp_Low := Exp_Low - 0.25; --  exact
 383                end if;
 384
 385             elsif Exp_Low >= 0.25 then
 386                Factor := Sqrt (Sqrt (X));
 387                Exp_Low := Exp_Low - 0.25; --  exact
 388             end if;
 389
 390             --  Exp_Low now is in interval (0.0, 0.25)
 391
 392             --  This means it is safe to call Logarithmic_Pow
 393             --  for the remaining part.
 394
 395             Factor := Factor * Logarithmic_Pow (X, Exp_Low);
 396          end if;
 397
 398       elsif X = 0.0 then
 399          return 0.0;
 400       end if;
 401
 402       --  Exp_High is non-zero integer smaller than 2**Double'Machine_Mantissa
 403
 404       Exp_Int := Mantissa_Type (Exp_High);
 405
 406       --  Standard way for processing integer powers > 0
 407
 408       while Exp_Int > 1 loop
 409          if (Exp_Int and 1) = 1 then
 410
 411             --  Base**Y = Base**(Exp_Int - 1) * Exp_Int for Exp_Int > 0
 412
 413             Factor := Factor * Base;
 414          end if;
 415
 416          --  Exp_Int is even and Exp_Int > 0, so
 417          --    Base**Y = (Base**2)**(Exp_Int / 2)
 418
 419          Base := Base * Base;
 420          Exp_Int := Exp_Int / 2;
 421       end loop;
 422
 423       --  Exp_Int = 1 or Exp_Int = 0
 424
 425       if Exp_Int = 1 then
 426          Factor := Base * Factor;
 427       end if;
 428
 429       if Negative_Y then
 430          Factor := 1.0 / Factor;
 431       end if;
 432
 433       return Factor;
 434    end Pow;
 435
 436    ---------
 437    -- Sin --
 438    ---------
 439
 440    function Sin (X : Double) return Double is
 441       Reduced_X : Double := X;
 442       Result    : Double;
 443       Quadrant  : Natural range 0 .. 3;
 444
 445    begin
 446       if abs X > Pi / 4.0 then
 447          Reduce (Reduced_X, Quadrant);
 448
 449          case Quadrant is
 450             when 0 =>
 451                Asm (Template  => "fsin",
 452                   Outputs  => Double'Asm_Output ("=t", Result),
 453                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 454
 455             when 1 =>
 456                Asm (Template  => "fcos",
 457                   Outputs  => Double'Asm_Output ("=t", Result),
 458                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 459
 460             when 2 =>
 461                Asm (Template  => "fsin",
 462                   Outputs  => Double'Asm_Output ("=t", Result),
 463                   Inputs   => Double'Asm_Input  ("0", -Reduced_X));
 464
 465             when 3 =>
 466                Asm (Template  => "fcos ; fchs",
 467                   Outputs  => Double'Asm_Output ("=t", Result),
 468                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 469          end case;
 470
 471       else
 472          Asm (Template  => "fsin",
 473             Outputs  => Double'Asm_Output ("=t", Result),
 474             Inputs   => Double'Asm_Input  ("0", Reduced_X));
 475       end if;
 476
 477       return Result;
 478    end Sin;
 479
 480    ---------
 481    -- Tan --
 482    ---------
 483
 484    function Tan (X : Double) return Double is
 485       Reduced_X : Double := X;
 486       Result    : Double;
 487       Quadrant  : Natural range 0 .. 3;
 488
 489    begin
 490       if abs X > Pi / 4.0 then
 491          Reduce (Reduced_X, Quadrant);
 492
 493          if Quadrant mod 2 = 0 then
 494             Asm (Template  => "fptan" & NL
 495                             & "ffree   %%st(0)"  & NL
 496                             & "fincstp",
 497                   Outputs  => Double'Asm_Output ("=t", Result),
 498                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 499          else
 500             Asm (Template  => "fsincos" & NL
 501                             & "fdivp   %%st, %%st(1)" & NL
 502                             & "fchs",
 503                   Outputs  => Double'Asm_Output ("=t", Result),
 504                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 505          end if;
 506
 507       else
 508          Asm (Template  =>
 509                "fptan                 " & NL
 510              & "ffree   %%st(0)      " & NL
 511              & "fincstp              ",
 512                Outputs  => Double'Asm_Output ("=t", Result),
 513                Inputs   => Double'Asm_Input  ("0", Reduced_X));
 514       end if;
 515
 516       return Result;
 517    end Tan;
 518
 519    ----------
 520    -- Sinh --
 521    ----------
 522
 523    function Sinh (X : Double) return Double is
 524    begin
 525       --  Mathematically Sinh (x) is defined to be (Exp (X) - Exp (-X)) / 2.0
 526
 527       if abs X < 25.0 then
 528          return (Exp (X) - Exp (-X)) / 2.0;
 529       else
 530          return Exp (X) / 2.0;
 531       end if;
 532    end Sinh;
 533
 534    ----------
 535    -- Cosh --
 536    ----------
 537
 538    function Cosh (X : Double) return Double is
 539    begin
 540       --  Mathematically Cosh (X) is defined to be (Exp (X) + Exp (-X)) / 2.0
 541
 542       if abs X < 22.0 then
 543          return (Exp (X) + Exp (-X)) / 2.0;
 544       else
 545          return Exp (X) / 2.0;
 546       end if;
 547    end Cosh;
 548
 549    ----------
 550    -- Tanh --
 551    ----------
 552
 553    function Tanh (X : Double) return Double is
 554    begin
 555       --  Return the Hyperbolic Tangent of x
 556
 557       --                                    x    -x
 558       --                                   e  - e        Sinh (X)
 559       --       Tanh (X) is defined to be -----------   = --------
 560       --                                    x    -x      Cosh (X)
 561       --                                   e  + e
 562
 563       if abs X > 23.0 then
 564          return Double'Copy_Sign (1.0, X);
 565       end if;
 566
 567       return 1.0 / (1.0 + Exp (-(2.0 * X))) - 1.0 / (1.0 + Exp (2.0 * X));
 568    end Tanh;
 569
 570 end Ada.Numerics.Aux;