a-numaux-darwin.adb [plain text]
package body Ada.Numerics.Aux is
procedure Reduce (X : in out Double; Q : out Natural);
function Sine_Approx (X : Double) return Double;
function Cosine_Approx (X : Double) return Double;
pragma Inline (Reduce);
pragma Inline (Sine_Approx);
pragma Inline (Cosine_Approx);
function Cosine_Approx (X : Double) return Double is
XX : constant Double := X * X;
begin
return (((((16#8.DC57FBD05F640#E-08 * XX
- 16#4.9F7D00BF25D80#E-06) * XX
+ 16#1.A019F7FDEFCC2#E-04) * XX
- 16#5.B05B058F18B20#E-03) * XX
+ 16#A.AAAAAAAA73FA8#E-02) * XX
- 16#7.FFFFFFFFFFDE4#E-01) * XX
- 16#3.655E64869ECCE#E-14 + 1.0;
end Cosine_Approx;
function Sine_Approx (X : Double) return Double is
XX : constant Double := X * X;
begin
return (((((16#A.EA2D4ABE41808#E-09 * XX
- 16#6.B974C10F9D078#E-07) * XX
+ 16#2.E3BC673425B0E#E-05) * XX
- 16#D.00D00CCA7AF00#E-04) * XX
+ 16#2.222222221B190#E-02) * XX
- 16#2.AAAAAAAAAAA44#E-01) * (XX * X) + X;
end Sine_Approx;
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 Cos (X : Double) return Double is
Reduced_X : Double := abs X;
Quadrant : Natural range 0 .. 3;
begin
if Reduced_X > Pi / 4.0 then
Reduce (Reduced_X, Quadrant);
case Quadrant is
when 0 =>
return Cosine_Approx (Reduced_X);
when 1 =>
return Sine_Approx (-Reduced_X);
when 2 =>
return -Cosine_Approx (Reduced_X);
when 3 =>
return Sine_Approx (Reduced_X);
end case;
end if;
return Cosine_Approx (Reduced_X);
end Cos;
function Sin (X : Double) return Double is
Reduced_X : Double := X;
Quadrant : Natural range 0 .. 3;
begin
if abs X > Pi / 4.0 then
Reduce (Reduced_X, Quadrant);
case Quadrant is
when 0 =>
return Sine_Approx (Reduced_X);
when 1 =>
return Cosine_Approx (Reduced_X);
when 2 =>
return Sine_Approx (-Reduced_X);
when 3 =>
return -Cosine_Approx (Reduced_X);
end case;
end if;
return Sine_Approx (Reduced_X);
end Sin;
end Ada.Numerics.Aux;