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