with System.Machine_Code; use System.Machine_Code;
package body Ada.Numerics.Aux is
NL : constant String := ASCII.LF & ASCII.HT;
function Is_Nan (X : Double) return Boolean;
function Logarithmic_Pow (X, Y : Double) return Double;
procedure Reduce (X : in out Double; Q : out Natural);
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;
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); 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
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
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;
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 := 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;
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;
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;
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(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;
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;