gcc/ada/a-numaux-x86.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUNTIME 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-2004 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 2,  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.  See the GNU General Public License --
  18 -- for  more details.  You should have  received  a copy of the GNU General --
  19 -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  20 -- to  the Free Software Foundation,  59 Temple Place - Suite 330,  Boston, --
  21 -- MA 02111-1307, USA.                                                      --
  22 --                                                                          --
  23 -- As a special exception,  if other files  instantiate  generics from this --
  24 -- unit, or you link  this unit with other files  to produce an executable, --
  25 -- this  unit  does not  by itself cause  the resulting  executable  to  be --
  26 -- covered  by the  GNU  General  Public  License.  This exception does not --
  27 -- however invalidate  any other reasons why  the executable file  might be --
  28 -- covered by the  GNU Public License.                                      --
  29 --                                                                          --
  30 -- GNAT was originally developed  by the GNAT team at  New York University. --
  31 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  32 --                                                                          --
  33 ------------------------------------------------------------------------------
  34
  35 --  File a-numaux.adb <- 86numaux.adb
  36
  37 --  This version of Numerics.Aux is for the IEEE Double Extended floating
  38 --  point format on x86.
  39
  40 with System.Machine_Code; use System.Machine_Code;
  41
  42 package body Ada.Numerics.Aux is
  43
  44    NL : constant String := ASCII.LF & ASCII.HT;
  45
  46    -----------------------
  47    -- Local subprograms --
  48    -----------------------
  49
  50    function Is_Nan (X : Double) return Boolean;
  51    --  Return True iff X is a IEEE NaN value
  52
  53    function Logarithmic_Pow (X, Y : Double) return Double;
  54    --  Implementation of X**Y using Exp and Log functions (binary base)
  55    --  to calculate the exponentiation. This is used by Pow for values
  56    --  for values of Y in the open interval (-0.25, 0.25)
  57
  58    procedure Reduce (X : in out Double; Q : out Natural);
  59    --  Implements reduction of X by Pi/2. Q is the quadrant of the final
  60    --  result in the range 0 .. 3. The absolute value of X is at most Pi.
  61
  62    pragma Inline (Is_Nan);
  63    pragma Inline (Reduce);
  64
  65    --------------------------------
  66    -- Basic Elementary Functions --
  67    --------------------------------
  68
  69    --  This section implements a few elementary functions that are used to
  70    --  build the more complex ones. This ordering enables better inlining.
  71
  72    ----------
  73    -- Atan --
  74    ----------
  75
  76    function Atan (X : Double) return Double is
  77       Result  : Double;
  78
  79    begin
  80       Asm (Template =>
  81            "fld1" & NL
  82          & "fpatan",
  83          Outputs  => Double'Asm_Output ("=t", Result),
  84          Inputs   => Double'Asm_Input  ("0", X));
  85
  86       --  The result value is NaN iff input was invalid
  87
  88       if not (Result = Result) then
  89          raise Argument_Error;
  90       end if;
  91
  92       return Result;
  93    end Atan;
  94
  95    ---------
  96    -- Exp --
  97    ---------
  98
  99    function Exp (X : Double) return Double is
 100       Result : Double;
 101    begin
 102       Asm (Template =>
 103          "fldl2e               " & NL
 104        & "fmulp   %%st, %%st(1)" & NL -- X * log2 (E)
 105        & "fld     %%st(0)      " & NL
 106        & "frndint              " & NL -- Integer (X * Log2 (E))
 107        & "fsubr   %%st, %%st(1)" & NL -- Fraction (X * Log2 (E))
 108        & "fxch                 " & NL
 109        & "f2xm1                " & NL -- 2**(...) - 1
 110        & "fld1                 " & NL
 111        & "faddp   %%st, %%st(1)" & NL -- 2**(Fraction (X * Log2 (E)))
 112        & "fscale               " & NL -- E ** X
 113        & "fstp    %%st(1)      ",
 114          Outputs  => Double'Asm_Output ("=t", Result),
 115          Inputs   => Double'Asm_Input  ("0", X));
 116       return Result;
 117    end Exp;
 118
 119    ------------
 120    -- Is_Nan --
 121    ------------
 122
 123    function Is_Nan (X : Double) return Boolean is
 124    begin
 125       --  The IEEE NaN values are the only ones that do not equal themselves
 126
 127       return not (X = X);
 128    end Is_Nan;
 129
 130    ---------
 131    -- Log --
 132    ---------
 133
 134    function Log (X : Double) return Double is
 135       Result : Double;
 136
 137    begin
 138       Asm (Template =>
 139          "fldln2               " & NL
 140        & "fxch                 " & NL
 141        & "fyl2x                " & NL,
 142          Outputs  => Double'Asm_Output ("=t", Result),
 143          Inputs   => Double'Asm_Input  ("0", X));
 144       return Result;
 145    end Log;
 146
 147    ------------
 148    -- Reduce --
 149    ------------
 150
 151    procedure Reduce (X : in out Double; Q : out Natural) is
 152       Half_Pi     : constant := Pi / 2.0;
 153       Two_Over_Pi : constant := 2.0 / Pi;
 154
 155       HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
 156       M  : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
 157       P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
 158       P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
 159       P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
 160       P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
 161       P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
 162                                                                  - P4, HM);
 163       P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
 164       K  : Double := X * Two_Over_Pi;
 165    begin
 166       --  For X < 2.0**32, all products below are computed exactly.
 167       --  Due to cancellation effects all subtractions are exact as well.
 168       --  As no double extended floating-point number has more than 75
 169       --  zeros after the binary point, the result will be the correctly
 170       --  rounded result of X - K * (Pi / 2.0).
 171
 172       while abs K >= 2.0**HM loop
 173          K := K * M - (K * M - K);
 174          X := (((((X - K * P1) - K * P2) - K * P3)
 175                      - K * P4) - K * P5) - K * P6;
 176          K := X * Two_Over_Pi;
 177       end loop;
 178
 179       if K /= K then
 180
 181          --  K is not a number, because X was not finite
 182
 183          raise Constraint_Error;
 184       end if;
 185
 186       K := Double'Rounding (K);
 187       Q := Integer (K) mod 4;
 188       X := (((((X - K * P1) - K * P2) - K * P3)
 189                   - 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             when 1 =>
 274                Asm (Template  => "fsin",
 275                   Outputs  => Double'Asm_Output ("=t", Result),
 276                   Inputs   => Double'Asm_Input  ("0", -Reduced_X));
 277             when 2 =>
 278                Asm (Template  => "fcos ; fchs",
 279                   Outputs  => Double'Asm_Output ("=t", Result),
 280                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 281             when 3 =>
 282                Asm (Template  => "fsin",
 283                   Outputs  => Double'Asm_Output ("=t", Result),
 284                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 285          end case;
 286
 287       else
 288          Asm (Template  => "fcos",
 289               Outputs  => Double'Asm_Output ("=t", Result),
 290               Inputs   => Double'Asm_Input  ("0", Reduced_X));
 291       end if;
 292
 293       return Result;
 294    end Cos;
 295
 296    ---------------------
 297    -- Logarithmic_Pow --
 298    ---------------------
 299
 300    function Logarithmic_Pow (X, Y : Double) return Double is
 301       Result  : Double;
 302    begin
 303       Asm (Template => ""             --  X                  : Y
 304        & "fyl2x                " & NL --  Y * Log2 (X)
 305        & "fst     %%st(1)      " & NL --  Y * Log2 (X)       : Y * Log2 (X)
 306        & "frndint              " & NL --  Int (...)          : Y * Log2 (X)
 307        & "fsubr   %%st, %%st(1)" & NL --  Int (...)          : Fract (...)
 308        & "fxch                 " & NL --  Fract (...)        : Int (...)
 309        & "f2xm1                " & NL --  2**Fract (...) - 1 : Int (...)
 310        & "fld1                 " & NL --  1 : 2**Fract (...) - 1 : Int (...)
 311        & "faddp   %%st, %%st(1)" & NL --  2**Fract (...)     : Int (...)
 312        & "fscale               " & NL --  2**(Fract (...) + Int (...))
 313        & "fstp    %%st(1)      ",
 314          Outputs  => Double'Asm_Output ("=t", Result),
 315          Inputs   =>
 316            (Double'Asm_Input  ("0", X),
 317             Double'Asm_Input  ("u", Y)));
 318       return Result;
 319    end Logarithmic_Pow;
 320
 321    ---------
 322    -- Pow --
 323    ---------
 324
 325    function Pow (X, Y : Double) return Double is
 326       type Mantissa_Type is mod 2**Double'Machine_Mantissa;
 327       --  Modular type that can hold all bits of the mantissa of Double
 328
 329       --  For negative exponents, do divide at the end of the processing
 330
 331       Negative_Y : constant Boolean := Y < 0.0;
 332       Abs_Y      : constant Double := abs Y;
 333
 334       --  During this function the following invariant is kept:
 335       --  X ** (abs Y) = Base**(Exp_High + Exp_Mid + Exp_Low) * Factor
 336
 337       Base : Double := X;
 338
 339       Exp_High : Double := Double'Floor (Abs_Y);
 340       Exp_Mid  : Double;
 341       Exp_Low  : Double;
 342       Exp_Int  : Mantissa_Type;
 343
 344       Factor : Double := 1.0;
 345
 346    begin
 347       --  Select algorithm for calculating Pow (integer cases fall through)
 348
 349       if Exp_High >= 2.0**Double'Machine_Mantissa then
 350
 351          --  In case of Y that is IEEE infinity, just raise constraint error
 352
 353          if Exp_High > Double'Safe_Last then
 354             raise Constraint_Error;
 355          end if;
 356
 357          --  Large values of Y are even integers and will stay integer
 358          --  after division by two.
 359
 360          loop
 361             --  Exp_Mid and Exp_Low are zero, so
 362             --    X**(abs Y) = Base ** Exp_High = (Base**2) ** (Exp_High / 2)
 363
 364             Exp_High := Exp_High / 2.0;
 365             Base := Base * Base;
 366             exit when Exp_High < 2.0**Double'Machine_Mantissa;
 367          end loop;
 368
 369       elsif Exp_High /= Abs_Y then
 370          Exp_Low := Abs_Y - Exp_High;
 371          Factor := 1.0;
 372
 373          if Exp_Low /= 0.0 then
 374
 375             --  Exp_Low now is in interval (0.0, 1.0)
 376             --  Exp_Mid := Double'Floor (Exp_Low * 4.0) / 4.0;
 377
 378             Exp_Mid := 0.0;
 379             Exp_Low := Exp_Low - Exp_Mid;
 380
 381             if Exp_Low >= 0.5 then
 382                Factor := Sqrt (X);
 383                Exp_Low := Exp_Low - 0.5;  -- exact
 384
 385                if Exp_Low >= 0.25 then
 386                   Factor := Factor * Sqrt (Factor);
 387                   Exp_Low := Exp_Low - 0.25; --  exact
 388                end if;
 389
 390             elsif Exp_Low >= 0.25 then
 391                Factor := Sqrt (Sqrt (X));
 392                Exp_Low := Exp_Low - 0.25; --  exact
 393             end if;
 394
 395             --  Exp_Low now is in interval (0.0, 0.25)
 396
 397             --  This means it is safe to call Logarithmic_Pow
 398             --  for the remaining part.
 399
 400             Factor := Factor * Logarithmic_Pow (X, Exp_Low);
 401          end if;
 402
 403       elsif X = 0.0 then
 404          return 0.0;
 405       end if;
 406
 407       --  Exp_High is non-zero integer smaller than 2**Double'Machine_Mantissa
 408
 409       Exp_Int := Mantissa_Type (Exp_High);
 410
 411       --  Standard way for processing integer powers > 0
 412
 413       while Exp_Int > 1 loop
 414          if (Exp_Int and 1) = 1 then
 415
 416             --  Base**Y = Base**(Exp_Int - 1) * Exp_Int for Exp_Int > 0
 417
 418             Factor := Factor * Base;
 419          end if;
 420
 421          --  Exp_Int is even and Exp_Int > 0, so
 422          --    Base**Y = (Base**2)**(Exp_Int / 2)
 423
 424          Base := Base * Base;
 425          Exp_Int := Exp_Int / 2;
 426       end loop;
 427
 428       --  Exp_Int = 1 or Exp_Int = 0
 429
 430       if Exp_Int = 1 then
 431          Factor := Base * Factor;
 432       end if;
 433
 434       if Negative_Y then
 435          Factor := 1.0 / Factor;
 436       end if;
 437
 438       return Factor;
 439    end Pow;
 440
 441    ---------
 442    -- Sin --
 443    ---------
 444
 445    function Sin (X : Double) return Double is
 446       Reduced_X : Double := X;
 447       Result    : Double;
 448       Quadrant  : Natural range 0 .. 3;
 449
 450    begin
 451       if abs X > Pi / 4.0 then
 452          Reduce (Reduced_X, Quadrant);
 453
 454          case Quadrant is
 455             when 0 =>
 456                Asm (Template  => "fsin",
 457                   Outputs  => Double'Asm_Output ("=t", Result),
 458                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 459             when 1 =>
 460                Asm (Template  => "fcos",
 461                   Outputs  => Double'Asm_Output ("=t", Result),
 462                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 463             when 2 =>
 464                Asm (Template  => "fsin",
 465                   Outputs  => Double'Asm_Output ("=t", Result),
 466                   Inputs   => Double'Asm_Input  ("0", -Reduced_X));
 467             when 3 =>
 468                Asm (Template  => "fcos ; fchs",
 469                   Outputs  => Double'Asm_Output ("=t", Result),
 470                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 471          end case;
 472
 473       else
 474          Asm (Template  => "fsin",
 475             Outputs  => Double'Asm_Output ("=t", Result),
 476             Inputs   => Double'Asm_Input  ("0", Reduced_X));
 477       end if;
 478
 479       return Result;
 480    end Sin;
 481
 482    ---------
 483    -- Tan --
 484    ---------
 485
 486    function Tan (X : Double) return Double is
 487       Reduced_X : Double := X;
 488       Result    : Double;
 489       Quadrant  : Natural range 0 .. 3;
 490
 491    begin
 492       if abs X > Pi / 4.0 then
 493          Reduce (Reduced_X, Quadrant);
 494
 495          if Quadrant mod 2 = 0 then
 496             Asm (Template  => "fptan" & NL
 497                             & "ffree   %%st(0)"  & NL
 498                             & "fincstp",
 499                   Outputs  => Double'Asm_Output ("=t", Result),
 500                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 501          else
 502             Asm (Template  => "fsincos" & NL
 503                             & "fdivp   %%st, %%st(1)" & NL
 504                             & "fchs",
 505                   Outputs  => Double'Asm_Output ("=t", Result),
 506                   Inputs   => Double'Asm_Input  ("0", Reduced_X));
 507          end if;
 508
 509       else
 510          Asm (Template  =>
 511                "fptan                 " & NL
 512              & "ffree   %%st(0)      " & NL
 513              & "fincstp              ",
 514                Outputs  => Double'Asm_Output ("=t", Result),
 515                Inputs   => Double'Asm_Input  ("0", Reduced_X));
 516       end if;
 517
 518       return Result;
 519    end Tan;
 520
 521    ----------
 522    -- Sinh --
 523    ----------
 524
 525    function Sinh (X : Double) return Double is
 526    begin
 527       --  Mathematically Sinh (x) is defined to be (Exp (X) - Exp (-X)) / 2.0
 528
 529       if abs X < 25.0 then
 530          return (Exp (X) - Exp (-X)) / 2.0;
 531       else
 532          return Exp (X) / 2.0;
 533       end if;
 534    end Sinh;
 535
 536    ----------
 537    -- Cosh --
 538    ----------
 539
 540    function Cosh (X : Double) return Double is
 541    begin
 542       --  Mathematically Cosh (X) is defined to be (Exp (X) + Exp (-X)) / 2.0
 543
 544       if abs X < 22.0 then
 545          return (Exp (X) + Exp (-X)) / 2.0;
 546       else
 547          return Exp (X) / 2.0;
 548       end if;
 549    end Cosh;
 550
 551    ----------
 552    -- Tanh --
 553    ----------
 554
 555    function Tanh (X : Double) return Double is
 556    begin
 557       --  Return the Hyperbolic Tangent of x
 558
 559       --                                    x    -x
 560       --                                   e  - e        Sinh (X)
 561       --       Tanh (X) is defined to be -----------   = --------
 562       --                                    x    -x      Cosh (X)
 563       --                                   e  + e
 564
 565       if abs X > 23.0 then
 566          return Double'Copy_Sign (1.0, X);
 567       end if;
 568
 569       return 1.0 / (1.0 + Exp (-2.0 * X)) - 1.0 / (1.0 + Exp (2.0 * X));
 570    end Tanh;
 571
 572 end Ada.Numerics.Aux;