------------------------------------------------------------------------------
-- --
--- GNAT RUNTIME COMPONENTS --
+-- 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-2001 Free Software Foundation, Inc. --
+-- Copyright (C) 1998-2007, 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- --
-- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License --
-- for more details. You should have received a copy of the GNU General --
-- Public License distributed with GNAT; see file COPYING. If not, write --
--- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, --
--- MA 02111-1307, USA. --
+-- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, --
+-- Boston, MA 02110-1301, USA. --
-- --
-- As a special exception, if other files instantiate generics from this --
-- unit, or you link this unit with other files to produce an executable, --
package body Ada.Numerics.Aux is
- NL : constant String := ASCII.LF & ASCII.HT;
-
- type FPU_Stack_Pointer is range 0 .. 7;
- for FPU_Stack_Pointer'Size use 3;
-
- type FPU_Status_Word is record
- B : Boolean; -- FPU Busy (for 8087 compatibility only)
- ES : Boolean; -- Error Summary Status
- SF : Boolean; -- Stack Fault
-
- Top : FPU_Stack_Pointer;
-
- -- Condition Code Flags
-
- -- C2 is set by FPREM and FPREM1 to indicate incomplete reduction.
- -- In case of successfull recorction, C0, C3 and C1 are set to the
- -- three least significant bits of the result (resp. Q2, Q1 and Q0).
-
- -- C2 is used by FPTAN, FSIN, FCOS, and FSINCOS to indicate that
- -- that source operand is beyond the allowable range of
- -- -2.0**63 .. 2.0**63.
-
- C3 : Boolean;
- C2 : Boolean;
- C1 : Boolean;
- C0 : Boolean;
-
- -- Exception Flags
-
- PE : Boolean; -- Precision
- UE : Boolean; -- Underflow
- OE : Boolean; -- Overflow
- ZE : Boolean; -- Zero Divide
- DE : Boolean; -- Denormalized Operand
- IE : Boolean; -- Invalid Operation
- end record;
-
- for FPU_Status_Word use record
- B at 0 range 15 .. 15;
- C3 at 0 range 14 .. 14;
- Top at 0 range 11 .. 13;
- C2 at 0 range 10 .. 10;
- C1 at 0 range 9 .. 9;
- C0 at 0 range 8 .. 8;
- ES at 0 range 7 .. 7;
- SF at 0 range 6 .. 6;
- PE at 0 range 5 .. 5;
- UE at 0 range 4 .. 4;
- OE at 0 range 3 .. 3;
- ZE at 0 range 2 .. 2;
- DE at 0 range 1 .. 1;
- IE at 0 range 0 .. 0;
- end record;
-
- for FPU_Status_Word'Size use 16;
+ NL : constant String := ASCII.LF & ASCII.HT;
-----------------------
-- Local subprograms --
-- to calculate the exponentiation. This is used by Pow for values
-- for values of Y in the open interval (-0.25, 0.25)
- function Reduce (X : Double) return Double;
- -- Implement partial reduction of X by Pi in the x86.
-
- -- Note that for the Sin, Cos and Tan functions completely accurate
- -- reduction of the argument is done for arguments in the range of
- -- -2.0**63 .. 2.0**63, using a 66-bit approximation of Pi.
+ 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 --
- ---------------------------------
+ --------------------------------
+ -- Basic Elementary Functions --
+ --------------------------------
- -- This section implements a few elementary functions that are
- -- used to build the more complex ones. This ordering enables
- -- better inlining.
+ -- This section implements a few elementary functions that are used to
+ -- build the more complex ones. This ordering enables better inlining.
----------
-- Atan --
-- Reduce --
------------
- function Reduce (X : Double) return Double is
- Result : Double;
+ 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
- Asm
- (Template =>
- -- Partial argument reduction
- "fldpi " & NL
- & "fadd %%st(0), %%st" & NL
- & "fxch %%st(1) " & NL
- & "fprem1 " & NL
- & "fstp %%st(1) ",
- Outputs => Double'Asm_Output ("=t", Result),
- Inputs => Double'Asm_Input ("0", X));
- return Result;
+ -- 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;
----------
return Result;
end Sqrt;
- ---------------------------------
- -- Other Elementary Functions --
- ---------------------------------
+ --------------------------------
+ -- Other Elementary Functions --
+ --------------------------------
-- These are built using the previously implemented basic functions
function Acos (X : Double) return Double is
Result : Double;
+
begin
Result := 2.0 * Atan (Sqrt ((1.0 - X) / (1.0 + X)));
function Asin (X : Double) return Double is
Result : Double;
- begin
+ begin
Result := Atan (X / Sqrt ((1.0 - X) * (1.0 + X)));
-- The result value is NaN iff input was invalid
---------
function Cos (X : Double) return Double is
- Reduced_X : Double := X;
+ Reduced_X : Double := abs X;
Result : Double;
- Status : FPU_Status_Word;
+ 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;
- loop
- Asm
- (Template =>
- "fcos " & NL
- & "xorl %%eax, %%eax " & NL
- & "fnstsw %%ax ",
- Outputs => (Double'Asm_Output ("=t", Result),
- FPU_Status_Word'Asm_Output ("=a", Status)),
- Inputs => Double'Asm_Input ("0", Reduced_X));
-
- exit when not Status.C2;
-
- -- Original argument was not in range and the result
- -- is the unmodified argument.
-
- Reduced_X := Reduce (Result);
- end loop;
+ else
+ Asm (Template => "fcos",
+ Outputs => Double'Asm_Output ("=t", Result),
+ Inputs => Double'Asm_Input ("0", Reduced_X));
+ end if;
return Result;
end Cos;
function Logarithmic_Pow (X, Y : Double) return Double is
Result : Double;
-
begin
Asm (Template => "" -- X : Y
& "fyl2x " & NL -- Y * Log2 (X)
- & "fst %%st(1) " & NL -- Y * Log2 (X) : 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 " & NL -- 2**(Fract (...) + Int (...))
- & "fstp %%st(1) ",
+ & "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;
type Mantissa_Type is mod 2**Double'Machine_Mantissa;
-- Modular type that can hold all bits of the mantissa of Double
- -- For negative exponents, a division is done
- -- at the end of the processing.
+ -- For negative exponents, do divide at the end of the processing
Negative_Y : constant Boolean := Y < 0.0;
Abs_Y : constant Double := abs Y;
Factor : Double := 1.0;
begin
- -- Select algorithm for calculating Pow:
- -- integer cases fall through
+ -- Select algorithm for calculating Pow (integer cases fall through)
if Exp_High >= 2.0**Double'Machine_Mantissa then
elsif Exp_High /= Abs_Y then
Exp_Low := Abs_Y - Exp_High;
-
Factor := 1.0;
if Exp_Low /= 0.0 then
function Sin (X : Double) return Double is
Reduced_X : Double := X;
Result : Double;
- Status : FPU_Status_Word;
+ 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;
- loop
- Asm
- (Template =>
- "fsin " & NL
- & "xorl %%eax, %%eax " & NL
- & "fnstsw %%ax ",
- Outputs => (Double'Asm_Output ("=t", Result),
- FPU_Status_Word'Asm_Output ("=a", Status)),
- Inputs => Double'Asm_Input ("0", Reduced_X));
-
- exit when not Status.C2;
-
- -- Original argument was not in range and the result
- -- is the unmodified argument.
-
- Reduced_X := Reduce (Result);
- end loop;
+ else
+ Asm (Template => "fsin",
+ Outputs => Double'Asm_Output ("=t", Result),
+ Inputs => Double'Asm_Input ("0", Reduced_X));
+ end if;
return Result;
end Sin;
function Tan (X : Double) return Double is
Reduced_X : Double := X;
Result : Double;
- Status : FPU_Status_Word;
+ 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;
- loop
- Asm
- (Template =>
- "fptan " & NL
- & "xorl %%eax, %%eax " & NL
- & "fnstsw %%ax " & NL
- & "ffree %%st(0) " & NL
- & "fincstp ",
-
- Outputs => (Double'Asm_Output ("=t", Result),
- FPU_Status_Word'Asm_Output ("=a", Status)),
- Inputs => Double'Asm_Input ("0", Reduced_X));
-
- exit when not Status.C2;
-
- -- Original argument was not in range and the result
- -- is the unmodified argument.
-
- Reduced_X := Reduce (Result);
- end loop;
+ 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;
if abs X < 25.0 then
return (Exp (X) - Exp (-X)) / 2.0;
-
else
return Exp (X) / 2.0;
end if;
-
end Sinh;
----------
if abs X < 22.0 then
return (Exp (X) + Exp (-X)) / 2.0;
-
else
return Exp (X) / 2.0;
end if;
-
end Cosh;
----------
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 ----------- = --------
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));
-
+ return 1.0 / (1.0 + Exp (-(2.0 * X))) - 1.0 / (1.0 + Exp (2.0 * X));
end Tanh;
end Ada.Numerics.Aux;
OSDN Git Service
