with System.Machine_Code; use System.Machine_Code;
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; ES : Boolean; SF : Boolean;
Top : FPU_Stack_Pointer;
C3 : Boolean;
C2 : Boolean;
C1 : Boolean;
C0 : Boolean;
PE : Boolean; UE : Boolean; OE : Boolean; ZE : Boolean; DE : Boolean; IE : Boolean; 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;
function Is_Nan (X : Double) return Boolean;
function Logarithmic_Pow (X, Y : Double) return Double;
function Reduce (X : Double) return Double;
pragma Inline (Is_Nan);
pragma Inline (Reduce);
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));
if not (Result = Result) then
raise Argument_Error;
end if;
return Result;
end Atan;
function Exp (X : Double) return Double is
Result : Double;
begin
Asm (Template =>
"fldl2e " & NL
& "fmulp %%st, %%st(1)" & NL & "fld %%st(0) " & NL
& "frndint " & NL & "fsubr %%st, %%st(1)" & NL & "fxch " & NL
& "f2xm1 " & NL & "fld1 " & NL
& "faddp %%st, %%st(1)" & NL & "fscale " & NL & "fstp %%st(1) ",
Outputs => Double'Asm_Output ("=t", Result),
Inputs => Double'Asm_Input ("0", X));
return Result;
end Exp;
function Is_Nan (X : Double) return Boolean is
begin
return not (X = X);
end Is_Nan;
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;
function Reduce (X : Double) return Double is
Result : Double;
begin
Asm
(Template =>
"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;
end Reduce;
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;
function Acos (X : Double) return Double is
Result : Double;
begin
Result := 2.0 * Atan (Sqrt ((1.0 - X) / (1.0 + X)));
if Is_Nan (Result) then
raise Argument_Error;
end if;
return Result;
end Acos;
function Asin (X : Double) return Double is
Result : Double;
begin
Result := Atan (X / Sqrt ((1.0 - X) * (1.0 + X)));
if Is_Nan (Result) then
raise Argument_Error;
end if;
return Result;
end Asin;
function Cos (X : Double) return Double is
Reduced_X : Double := X;
Result : Double;
Status : FPU_Status_Word;
begin
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;
Reduced_X := Reduce (Result);
end loop;
return Result;
end Cos;
function Logarithmic_Pow (X, Y : Double) return Double is
Result : Double;
begin
Asm (Template => "" & "fyl2x " & NL & "fst %%st(1) " & NL & "frndint " & NL & "fsubr %%st, %%st(1)" & NL & "fxch " & NL & "f2xm1 " & NL & "fld1 " & NL & "faddp %%st, %%st(1)" & NL & "fscale " & NL & "fstp %%st(1) ",
Outputs => Double'Asm_Output ("=t", Result),
Inputs =>
(Double'Asm_Input ("0", X),
Double'Asm_Input ("u", Y)));
return Result;
end Logarithmic_Pow;
function Pow (X, Y : Double) return Double is
type Mantissa_Type is mod 2**Double'Machine_Mantissa;
Negative_Y : constant Boolean := Y < 0.0;
Abs_Y : constant Double := abs Y;
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
if Exp_High >= 2.0**Double'Machine_Mantissa then
if Exp_High > Double'Safe_Last then
raise Constraint_Error;
end if;
loop
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_Mid := 0.0;
Exp_Low := Exp_Low - Exp_Mid;
if Exp_Low >= 0.5 then
Factor := Sqrt (X);
Exp_Low := Exp_Low - 0.5;
if Exp_Low >= 0.25 then
Factor := Factor * Sqrt (Factor);
Exp_Low := Exp_Low - 0.25; end if;
elsif Exp_Low >= 0.25 then
Factor := Sqrt (Sqrt (X));
Exp_Low := Exp_Low - 0.25; end if;
Factor := Factor * Logarithmic_Pow (X, Exp_Low);
end if;
elsif X = 0.0 then
return 0.0;
end if;
Exp_Int := Mantissa_Type (Exp_High);
while Exp_Int > 1 loop
if (Exp_Int and 1) = 1 then
Factor := Factor * Base;
end if;
Base := Base * Base;
Exp_Int := Exp_Int / 2;
end loop;
if Exp_Int = 1 then
Factor := Base * Factor;
end if;
if Negative_Y then
Factor := 1.0 / Factor;
end if;
return Factor;
end Pow;
function Sin (X : Double) return Double is
Reduced_X : Double := X;
Result : Double;
Status : FPU_Status_Word;
begin
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;
Reduced_X := Reduce (Result);
end loop;
return Result;
end Sin;
function Tan (X : Double) return Double is
Reduced_X : Double := X;
Result : Double;
Status : FPU_Status_Word;
begin
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;
Reduced_X := Reduce (Result);
end loop;
return Result;
end Tan;
function Sinh (X : Double) return Double is
begin
if abs X < 25.0 then
return (Exp (X) - Exp (-X)) / 2.0;
else
return Exp (X) / 2.0;
end if;
end Sinh;
function Cosh (X : Double) return Double is
begin
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
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;