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