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