URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [gcc/] [ada/] [a-numaux-x86.adb] - Blame information for rev 706

Details | Compare with Previous | View Log


------------------------------------------------------------------------------
--                                                                          --
--                         GNAT RUN-TIME COMPONENTS                         --
--                                                                          --
--                     A D A . N U M E R I C S . A U X                      --
--                                                                          --
--                                 B o d y                                  --
--                        (Machine Version for x86)                         --
--                                                                          --
--          Copyright (C) 1998-2009, Free Software Foundation, Inc.         --
--                                                                          --
-- GNAT is free software;  you can  redistribute it  and/or modify it under --
-- terms of the  GNU General Public License as published  by the Free Soft- --
-- ware  Foundation;  either version 3,  or (at your option) any later ver- --
-- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
-- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
-- or FITNESS FOR A PARTICULAR PURPOSE.                                     --
--                                                                          --
-- As a special exception under Section 7 of GPL version 3, you are granted --
-- additional permissions described in the GCC Runtime Library Exception,   --
-- version 3.1, as published by the Free Software Foundation.               --
--                                                                          --
-- You should have received a copy of the GNU General Public License and    --
-- a copy of the GCC Runtime Library Exception along with this program;     --
-- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --
-- <http://www.gnu.org/licenses/>.                                          --
--                                                                          --
-- GNAT was originally developed  by the GNAT team at  New York University. --
-- Extensive contributions were provided by Ada Core Technologies Inc.      --
--                                                                          --
------------------------------------------------------------------------------
 
--  File a-numaux.adb <- 86numaux.adb
 
--  This version of Numerics.Aux is for the IEEE Double Extended floating
--  point format on x86.
 
with System.Machine_Code; use System.Machine_Code;
 
package body Ada.Numerics.Aux is
 
   NL : constant String := ASCII.LF & ASCII.HT;
 
   -----------------------
   -- Local subprograms --
   -----------------------
 
   function Is_Nan (X : Double) return Boolean;
   --  Return True iff X is a IEEE NaN value
 
   function Logarithmic_Pow (X, Y : Double) return Double;
   --  Implementation of X**Y using Exp and Log functions (binary base)
   --  to calculate the exponentiation. This is used by Pow for values
   --  for values of Y in the open interval (-0.25, 0.25)
 
   procedure Reduce (X : in out Double; Q : out Natural);
   --  Implements reduction of X by Pi/2. Q is the quadrant of the final
   --  result in the range 0 .. 3. The absolute value of X is at most Pi.
 
   pragma Inline (Is_Nan);
   pragma Inline (Reduce);
 
   --------------------------------
   -- Basic Elementary Functions --
   --------------------------------
 
   --  This section implements a few elementary functions that are used to
   --  build the more complex ones. This ordering enables better inlining.
 
   ----------
   -- Atan --
   ----------
 
   function Atan (X : Double) return Double is
      Result  : Double;
 
   begin
      Asm (Template =>
           "fld1" & NL
         & "fpatan",
         Outputs  => Double'Asm_Output ("=t", Result),
         Inputs   => Double'Asm_Input  ("0", X));
 
      --  The result value is NaN iff input was invalid
 
      if not (Result = Result) then
         raise Argument_Error;
      end if;
 
      return Result;
   end Atan;
 
   ---------
   -- Exp --
   ---------
 
   function Exp (X : Double) return Double is
      Result : Double;
   begin
      Asm (Template =>
         "fldl2e               " & NL
       & "fmulp   %%st, %%st(1)" & NL -- X * log2 (E)
       & "fld     %%st(0)      " & NL
       & "frndint              " & NL -- Integer (X * Log2 (E))
       & "fsubr   %%st, %%st(1)" & NL -- Fraction (X * Log2 (E))
       & "fxch                 " & NL
       & "f2xm1                " & NL -- 2**(...) - 1
       & "fld1                 " & NL
       & "faddp   %%st, %%st(1)" & NL -- 2**(Fraction (X * Log2 (E)))
       & "fscale               " & NL -- E ** X
       & "fstp    %%st(1)      ",
         Outputs  => Double'Asm_Output ("=t", Result),
         Inputs   => Double'Asm_Input  ("0", X));
      return Result;
   end Exp;
 
   ------------
   -- Is_Nan --
   ------------
 
   function Is_Nan (X : Double) return Boolean is
   begin
      --  The IEEE NaN values are the only ones that do not equal themselves
 
      return not (X = X);
   end Is_Nan;
 
   ---------
   -- Log --
   ---------
 
   function Log (X : Double) return Double is
      Result : Double;
 
   begin
      Asm (Template =>
         "fldln2               " & NL
       & "fxch                 " & NL
       & "fyl2x                " & NL,
         Outputs  => Double'Asm_Output ("=t", Result),
         Inputs   => Double'Asm_Input  ("0", X));
      return Result;
   end Log;
 
   ------------
   -- Reduce --
   ------------
 
   procedure Reduce (X : in out Double; Q : out Natural) is
      Half_Pi     : constant := Pi / 2.0;
      Two_Over_Pi : constant := 2.0 / Pi;
 
      HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
      M  : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
      P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
      P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
      P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
      P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
      P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
                                                                 - P4, HM);
      P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
      K  : Double := X * Two_Over_Pi;
   begin
      --  For X < 2.0**32, all products below are computed exactly.
      --  Due to cancellation effects all subtractions are exact as well.
      --  As no double extended floating-point number has more than 75
      --  zeros after the binary point, the result will be the correctly
      --  rounded result of X - K * (Pi / 2.0).
 
      while abs K >= 2.0**HM loop
         K := K * M - (K * M - K);
         X := (((((X - K * P1) - K * P2) - K * P3)
                     - K * P4) - K * P5) - K * P6;
         K := X * Two_Over_Pi;
      end loop;
 
      if K /= K then
 
         --  K is not a number, because X was not finite
 
         raise Constraint_Error;
      end if;
 
      K := Double'Rounding (K);
      Q := Integer (K) mod 4;
      X := (((((X - K * P1) - K * P2) - K * P3)
                  - K * P4) - K * P5) - K * P6;
   end Reduce;
 
   ----------
   -- Sqrt --
   ----------
 
   function Sqrt (X : Double) return Double is
      Result  : Double;
 
   begin
      if X < 0.0 then
         raise Argument_Error;
      end if;
 
      Asm (Template => "fsqrt",
           Outputs  => Double'Asm_Output ("=t", Result),
           Inputs   => Double'Asm_Input  ("0", X));
 
      return Result;
   end Sqrt;
 
   --------------------------------
   -- Other Elementary Functions --
   --------------------------------
 
   --  These are built using the previously implemented basic functions
 
   ----------
   -- Acos --
   ----------
 
   function Acos (X : Double) return Double is
      Result  : Double;
 
   begin
      Result := 2.0 * Atan (Sqrt ((1.0 - X) / (1.0 + X)));
 
      --  The result value is NaN iff input was invalid
 
      if Is_Nan (Result) then
         raise Argument_Error;
      end if;
 
      return Result;
   end Acos;
 
   ----------
   -- Asin --
   ----------
 
   function Asin (X : Double) return Double is
      Result  : Double;
 
   begin
      Result := Atan (X / Sqrt ((1.0 - X) * (1.0 + X)));
 
      --  The result value is NaN iff input was invalid
 
      if Is_Nan (Result) then
         raise Argument_Error;
      end if;
 
      return Result;
   end Asin;
 
   ---------
   -- Cos --
   ---------
 
   function Cos (X : Double) return Double is
      Reduced_X : Double := abs X;
      Result    : Double;
      Quadrant  : Natural range 0 .. 3;
 
   begin
      if Reduced_X > Pi / 4.0 then
         Reduce (Reduced_X, Quadrant);
 
         case Quadrant is
            when 0 =>
               Asm (Template  => "fcos",
                  Outputs  => Double'Asm_Output ("=t", Result),
                  Inputs   => Double'Asm_Input  ("0", Reduced_X));
            when 1 =>
               Asm (Template  => "fsin",
                  Outputs  => Double'Asm_Output ("=t", Result),
                  Inputs   => Double'Asm_Input  ("0", -Reduced_X));
            when 2 =>
               Asm (Template  => "fcos ; fchs",
                  Outputs  => Double'Asm_Output ("=t", Result),
                  Inputs   => Double'Asm_Input  ("0", Reduced_X));
            when 3 =>
               Asm (Template  => "fsin",
                  Outputs  => Double'Asm_Output ("=t", Result),
                  Inputs   => Double'Asm_Input  ("0", Reduced_X));
         end case;
 
      else
         Asm (Template  => "fcos",
              Outputs  => Double'Asm_Output ("=t", Result),
              Inputs   => Double'Asm_Input  ("0", Reduced_X));
      end if;
 
      return Result;
   end Cos;
 
   ---------------------
   -- Logarithmic_Pow --
   ---------------------
 
   function Logarithmic_Pow (X, Y : Double) return Double is
      Result  : Double;
   begin
      Asm (Template => ""             --  X                  : Y
       & "fyl2x                " & NL --  Y * Log2 (X)
       & "fld     %%st(0)      " & NL --  Y * Log2 (X)       : Y * Log2 (X)
       & "frndint              " & NL --  Int (...)          : Y * Log2 (X)
       & "fsubr   %%st, %%st(1)" & NL --  Int (...)          : Fract (...)
       & "fxch                 " & NL --  Fract (...)        : Int (...)
       & "f2xm1                " & NL --  2**Fract (...) - 1 : Int (...)
       & "fld1                 " & NL --  1 : 2**Fract (...) - 1 : Int (...)
       & "faddp   %%st, %%st(1)" & NL --  2**Fract (...)     : Int (...)
       & "fscale               ",     --  2**(Fract (...) + Int (...))
         Outputs  => Double'Asm_Output ("=t", Result),
         Inputs   =>
           (Double'Asm_Input  ("0", X),
            Double'Asm_Input  ("u", Y)));
      return Result;
   end Logarithmic_Pow;
 
   ---------
   -- Pow --
   ---------
 
   function Pow (X, Y : Double) return Double is
      type Mantissa_Type is mod 2**Double'Machine_Mantissa;
      --  Modular type that can hold all bits of the mantissa of Double
 
      --  For negative exponents, do divide at the end of the processing
 
      Negative_Y : constant Boolean := Y < 0.0;
      Abs_Y      : constant Double := abs Y;
 
      --  During this function the following invariant is kept:
      --  X ** (abs Y) = Base**(Exp_High + Exp_Mid + Exp_Low) * Factor
 
      Base : Double := X;
 
      Exp_High : Double := Double'Floor (Abs_Y);
      Exp_Mid  : Double;
      Exp_Low  : Double;
      Exp_Int  : Mantissa_Type;
 
      Factor : Double := 1.0;
 
   begin
      --  Select algorithm for calculating Pow (integer cases fall through)
 
      if Exp_High >= 2.0**Double'Machine_Mantissa then
 
         --  In case of Y that is IEEE infinity, just raise constraint error
 
         if Exp_High > Double'Safe_Last then
            raise Constraint_Error;
         end if;
 
         --  Large values of Y are even integers and will stay integer
         --  after division by two.
 
         loop
            --  Exp_Mid and Exp_Low are zero, so
            --    X**(abs Y) = Base ** Exp_High = (Base**2) ** (Exp_High / 2)
 
            Exp_High := Exp_High / 2.0;
            Base := Base * Base;
            exit when Exp_High < 2.0**Double'Machine_Mantissa;
         end loop;
 
      elsif Exp_High /= Abs_Y then
         Exp_Low := Abs_Y - Exp_High;
         Factor := 1.0;
 
         if Exp_Low /= 0.0 then
 
            --  Exp_Low now is in interval (0.0, 1.0)
            --  Exp_Mid := Double'Floor (Exp_Low * 4.0) / 4.0;
 
            Exp_Mid := 0.0;
            Exp_Low := Exp_Low - Exp_Mid;
 
            if Exp_Low >= 0.5 then
               Factor := Sqrt (X);
               Exp_Low := Exp_Low - 0.5;  -- exact
 
               if Exp_Low >= 0.25 then
                  Factor := Factor * Sqrt (Factor);
                  Exp_Low := Exp_Low - 0.25; --  exact
               end if;
 
            elsif Exp_Low >= 0.25 then
               Factor := Sqrt (Sqrt (X));
               Exp_Low := Exp_Low - 0.25; --  exact
            end if;
 
            --  Exp_Low now is in interval (0.0, 0.25)
 
            --  This means it is safe to call Logarithmic_Pow
            --  for the remaining part.
 
            Factor := Factor * Logarithmic_Pow (X, Exp_Low);
         end if;
 
      elsif X = 0.0 then
         return 0.0;
      end if;
 
      --  Exp_High is non-zero integer smaller than 2**Double'Machine_Mantissa
 
      Exp_Int := Mantissa_Type (Exp_High);
 
      --  Standard way for processing integer powers > 0
 
      while Exp_Int > 1 loop
         if (Exp_Int and 1) = 1 then
 
            --  Base**Y = Base**(Exp_Int - 1) * Exp_Int for Exp_Int > 0
 
            Factor := Factor * Base;
         end if;
 
         --  Exp_Int is even and Exp_Int > 0, so
         --    Base**Y = (Base**2)**(Exp_Int / 2)
 
         Base := Base * Base;
         Exp_Int := Exp_Int / 2;
      end loop;
 
      --  Exp_Int = 1 or Exp_Int = 0
 
      if Exp_Int = 1 then
         Factor := Base * Factor;
      end if;
 
      if Negative_Y then
         Factor := 1.0 / Factor;
      end if;
 
      return Factor;
   end Pow;
 
   ---------
   -- Sin --
   ---------
 
   function Sin (X : Double) return Double is
      Reduced_X : Double := X;
      Result    : Double;
      Quadrant  : Natural range 0 .. 3;
 
   begin
      if abs X > Pi / 4.0 then
         Reduce (Reduced_X, Quadrant);
 
         case Quadrant is
            when 0 =>
               Asm (Template  => "fsin",
                  Outputs  => Double'Asm_Output ("=t", Result),
                  Inputs   => Double'Asm_Input  ("0", Reduced_X));
            when 1 =>
               Asm (Template  => "fcos",
                  Outputs  => Double'Asm_Output ("=t", Result),
                  Inputs   => Double'Asm_Input  ("0", Reduced_X));
            when 2 =>
               Asm (Template  => "fsin",
                  Outputs  => Double'Asm_Output ("=t", Result),
                  Inputs   => Double'Asm_Input  ("0", -Reduced_X));
            when 3 =>
               Asm (Template  => "fcos ; fchs",
                  Outputs  => Double'Asm_Output ("=t", Result),
                  Inputs   => Double'Asm_Input  ("0", Reduced_X));
         end case;
 
      else
         Asm (Template  => "fsin",
            Outputs  => Double'Asm_Output ("=t", Result),
            Inputs   => Double'Asm_Input  ("0", Reduced_X));
      end if;
 
      return Result;
   end Sin;
 
   ---------
   -- Tan --
   ---------
 
   function Tan (X : Double) return Double is
      Reduced_X : Double := X;
      Result    : Double;
      Quadrant  : Natural range 0 .. 3;
 
   begin
      if abs X > Pi / 4.0 then
         Reduce (Reduced_X, Quadrant);
 
         if Quadrant mod 2 = 0 then
            Asm (Template  => "fptan" & NL
                            & "ffree   %%st(0)"  & NL
                            & "fincstp",
                  Outputs  => Double'Asm_Output ("=t", Result),
                  Inputs   => Double'Asm_Input  ("0", Reduced_X));
         else
            Asm (Template  => "fsincos" & NL
                            & "fdivp   %%st, %%st(1)" & NL
                            & "fchs",
                  Outputs  => Double'Asm_Output ("=t", Result),
                  Inputs   => Double'Asm_Input  ("0", Reduced_X));
         end if;
 
      else
         Asm (Template  =>
               "fptan                 " & NL
             & "ffree   %%st(0)      " & NL
             & "fincstp              ",
               Outputs  => Double'Asm_Output ("=t", Result),
               Inputs   => Double'Asm_Input  ("0", Reduced_X));
      end if;
 
      return Result;
   end Tan;
 
   ----------
   -- Sinh --
   ----------
 
   function Sinh (X : Double) return Double is
   begin
      --  Mathematically Sinh (x) is defined to be (Exp (X) - Exp (-X)) / 2.0
 
      if abs X < 25.0 then
         return (Exp (X) - Exp (-X)) / 2.0;
      else
         return Exp (X) / 2.0;
      end if;
   end Sinh;
 
   ----------
   -- Cosh --
   ----------
 
   function Cosh (X : Double) return Double is
   begin
      --  Mathematically Cosh (X) is defined to be (Exp (X) + Exp (-X)) / 2.0
 
      if abs X < 22.0 then
         return (Exp (X) + Exp (-X)) / 2.0;
      else
         return Exp (X) / 2.0;
      end if;
   end Cosh;
 
   ----------
   -- Tanh --
   ----------
 
   function Tanh (X : Double) return Double is
   begin
      --  Return the Hyperbolic Tangent of x
 
      --                                    x    -x
      --                                   e  - e        Sinh (X)
      --       Tanh (X) is defined to be -----------   = --------
      --                                    x    -x      Cosh (X)
      --                                   e  + e
 
      if abs X > 23.0 then
         return Double'Copy_Sign (1.0, X);
      end if;
 
      return 1.0 / (1.0 + Exp (-(2.0 * X))) - 1.0 / (1.0 + Exp (2.0 * X));
   end Tanh;
 
end Ada.Numerics.Aux;

Browse

Tools

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [gcc/] [ada/] [a-numaux-x86.adb] - Blame information for rev 706

Line No.	Rev	Author	Line
1	706	jeremybenn	`------------------------------------------------------------------------------`
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 = (Base2) (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			`-- BaseY = 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			`-- BaseY = (Base2)**(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;`