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