with Ada.Numerics.Aux;
use type Ada.Numerics.Aux.Double;
with Ada.Numerics; use Ada.Numerics;
package body Math_Lib is
Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
Two_Pi : constant Real'Base := 2.0 * Pi;
Half_Pi : constant Real'Base := Pi / 2.0;
Fourth_Pi : constant Real'Base := Pi / 4.0;
Epsilon : constant Real'Base := Real'Base'Epsilon;
IEpsilon : constant Real'Base := 1.0 / Epsilon;
subtype Double is Aux.Double;
DEpsilon : constant Double := Double (Epsilon);
DIEpsilon : constant Double := Double (IEpsilon);
function Arctan
(Y : Real;
A : Real := 1.0)
return Real;
function Arctan
(Y : Real;
A : Real := 1.0;
Cycle : Real)
return Real;
function Exact_Remainder
(A : Real;
Y : Real)
return Real;
function Half_Log_Epsilon return Real;
function Local_Atan
(Y : Real;
A : Real := 1.0)
return Real;
function Log_Inverse_Epsilon return Real;
function Square_Root_Epsilon return Real;
function "**" (A1, A2 : Real) return Real is
begin
if A1 = 0.0
and then A2 = 0.0
then
raise Argument_Error;
elsif A1 < 0.0 then
raise Argument_Error;
elsif A2 = 0.0 then
return 1.0;
elsif A1 = 0.0 then
if A2 < 0.0 then
raise Constraint_Error;
else
return 0.0;
end if;
elsif A1 = 1.0 then
return 1.0;
elsif A2 = 1.0 then
return A1;
else
begin
if A2 = 2.0 then
return A1 * A1;
else
return
Real (Aux.pow (Double (A1), Double (A2)));
end if;
exception
when others =>
raise Constraint_Error;
end;
end if;
end "**";
function Arccos (A : Real) return Real is
Temp : Real'Base;
begin
if abs A > 1.0 then
raise Argument_Error;
elsif abs A < Square_Root_Epsilon then
return Pi / 2.0 - A;
elsif A = 1.0 then
return 0.0;
elsif A = -1.0 then
return Pi;
end if;
Temp := Real (Aux.acos (Double (A)));
if Temp < 0.0 then
Temp := Pi + Temp;
end if;
return Temp;
end Arccos;
function Arccos (A, Cycle : Real) return Real is
Temp : Real'Base;
begin
if Cycle <= 0.0 then
raise Argument_Error;
elsif abs A > 1.0 then
raise Argument_Error;
elsif abs A < Square_Root_Epsilon then
return Cycle / 4.0;
elsif A = 1.0 then
return 0.0;
elsif A = -1.0 then
return Cycle / 2.0;
end if;
Temp := Arctan (Sqrt (1.0 - A * A) / A, 1.0, Cycle);
if Temp < 0.0 then
Temp := Cycle / 2.0 + Temp;
end if;
return Temp;
end Arccos;
function Arccosh (A : Real) return Real is
begin
if A < 1.0 then
raise Argument_Error;
elsif A < 1.0 + Square_Root_Epsilon then
return A - 1.0;
elsif abs A > 1.0 / Square_Root_Epsilon then
return Log (A) + Log_Two;
else
return Log (A + Sqrt (A * A - 1.0));
end if;
end Arccosh;
function Arccot
(A : Real;
Y : Real := 1.0)
return Real
is
begin
return Arctan (Y, A);
end Arccot;
function Arccot
(A : Real;
Y : Real := 1.0;
Cycle : Real)
return Real
is
begin
return Arctan (Y, A, Cycle);
end Arccot;
function Arccoth (A : Real) return Real is
begin
if abs A = 1.0 then
raise Constraint_Error;
elsif abs A < 1.0 then
raise Argument_Error;
elsif abs A > 1.0 / Epsilon then
return 0.0;
else
return 0.5 * Log ((1.0 + A) / (A - 1.0));
end if;
end Arccoth;
function Arcsin (A : Real) return Real is
begin
if abs A > 1.0 then
raise Argument_Error;
elsif abs A < Square_Root_Epsilon then
return A;
elsif A = 1.0 then
return Pi / 2.0;
elsif A = -1.0 then
return -Pi / 2.0;
end if;
return Real (Aux.asin (Double (A)));
end Arcsin;
function Arcsin (A, Cycle : Real) return Real is
begin
if Cycle <= 0.0 then
raise Argument_Error;
elsif abs A > 1.0 then
raise Argument_Error;
elsif A = 0.0 then
return A;
elsif A = 1.0 then
return Cycle / 4.0;
elsif A = -1.0 then
return -Cycle / 4.0;
end if;
return Arctan (A / Sqrt (1.0 - A * A), 1.0, Cycle);
end Arcsin;
function Arcsinh (A : Real) return Real is
begin
if abs A < Square_Root_Epsilon then
return A;
elsif A > 1.0 / Square_Root_Epsilon then
return Log (A) + Log_Two;
elsif A < -1.0 / Square_Root_Epsilon then
return -(Log (-A) + Log_Two);
elsif A < 0.0 then
return -Log (abs A + Sqrt (A * A + 1.0));
else
return Log (A + Sqrt (A * A + 1.0));
end if;
end Arcsinh;
function Arctan
(Y : Real;
A : Real := 1.0)
return Real
is
begin
if A = 0.0
and then Y = 0.0
then
raise Argument_Error;
elsif Y = 0.0 then
if A > 0.0 then
return 0.0;
else return Pi;
end if;
elsif A = 0.0 then
if Y > 0.0 then
return Half_Pi;
else return -Half_Pi;
end if;
else
return Local_Atan (Y, A);
end if;
end Arctan;
function Arctan
(Y : Real;
A : Real := 1.0;
Cycle : Real)
return Real
is
begin
if Cycle <= 0.0 then
raise Argument_Error;
elsif A = 0.0
and then Y = 0.0
then
raise Argument_Error;
elsif Y = 0.0 then
if A > 0.0 then
return 0.0;
else return Cycle / 2.0;
end if;
elsif A = 0.0 then
if Y > 0.0 then
return Cycle / 4.0;
else return -Cycle / 4.0;
end if;
else
return Local_Atan (Y, A) * Cycle / Two_Pi;
end if;
end Arctan;
function Arctanh (A : Real) return Real is
begin
if abs A = 1.0 then
raise Constraint_Error;
elsif abs A > 1.0 then
raise Argument_Error;
elsif abs A < Square_Root_Epsilon then
return A;
else
return 0.5 * Log ((1.0 + A) / (1.0 - A));
end if;
end Arctanh;
function Cos (A : Real) return Real is
begin
if A = 0.0 then
return 1.0;
elsif abs A < Square_Root_Epsilon then
return 1.0;
end if;
return Real (Aux.Cos (Double (A)));
end Cos;
function Cos (A, Cycle : Real) return Real is
T : Real'Base;
begin
if Cycle <= 0.0 then
raise Argument_Error;
elsif A = 0.0 then
return 1.0;
end if;
T := Exact_Remainder (abs (A), Cycle) / Cycle;
if T = 0.25
or else T = 0.75
or else T = -0.25
or else T = -0.75
then
return 0.0;
elsif T = 0.5 or T = -0.5 then
return -1.0;
end if;
return Real (Aux.Cos (Double (T * Two_Pi)));
end Cos;
function Cosh (A : Real) return Real is
begin
if abs A < Square_Root_Epsilon then
return 1.0;
elsif abs A > Log_Inverse_Epsilon then
return Exp ((abs A) - Log_Two);
end if;
return Real (Aux.cosh (Double (A)));
exception
when others =>
raise Constraint_Error;
end Cosh;
function Cot (A : Real) return Real is
begin
if A = 0.0 then
raise Constraint_Error;
elsif abs A < Square_Root_Epsilon then
return 1.0 / A;
end if;
return Real (1.0 / Real'Base (Aux.tan (Double (A))));
end Cot;
function Cot (A, Cycle : Real) return Real is
T : Real'Base;
begin
if Cycle <= 0.0 then
raise Argument_Error;
elsif A = 0.0 then
raise Constraint_Error;
elsif abs A < Square_Root_Epsilon then
return 1.0 / A;
end if;
T := Exact_Remainder (A, Cycle) / Cycle;
if T = 0.0 or T = 0.5 or T = -0.5 then
raise Constraint_Error;
else
return Cos (T * Two_Pi) / Sin (T * Two_Pi);
end if;
end Cot;
function Coth (A : Real) return Real is
begin
if A = 0.0 then
raise Constraint_Error;
elsif A < Half_Log_Epsilon then
return -1.0;
elsif A > -Half_Log_Epsilon then
return 1.0;
elsif abs A < Square_Root_Epsilon then
return 1.0 / A;
end if;
return Real (1.0 / Real'Base (Aux.tanh (Double (A))));
end Coth;
function Exact_Remainder
(A : Real;
Y : Real)
return Real
is
Denominator : Real'Base := abs A;
Divisor : Real'Base := abs Y;
Reducer : Real'Base;
Sign : Real'Base := 1.0;
begin
if Y = 0.0 then
raise Constraint_Error;
elsif A = 0.0 then
return 0.0;
elsif A = Y then
return 0.0;
elsif Denominator < Divisor then
return A;
end if;
while Denominator >= Divisor loop
Reducer := Divisor;
begin
while Reducer * 1_048_576.0 < Denominator loop
Reducer := Reducer * 1_048_576.0;
end loop;
exception
when others => null;
end;
begin
while Reducer * 1_024.0 < Denominator loop
Reducer := Reducer * 1_024.0;
end loop;
exception
when others => null;
end;
begin
while Reducer * 2.0 < Denominator loop
Reducer := Reducer * 2.0;
end loop;
exception
when others => null;
end;
Denominator := Denominator - Reducer;
end loop;
if A < 0.0 then
return -Denominator;
else
return Denominator;
end if;
end Exact_Remainder;
function Exp (A : Real) return Real is
Result : Real'Base;
begin
if A = 0.0 then
return 1.0;
else
Result := Real (Aux.Exp (Double (A)));
if Result > Real'Last then
raise Constraint_Error;
else
return Result;
end if;
end if;
end Exp;
function Half_Log_Epsilon return Real is
begin
return Real (0.5 * Real'Base (Aux.Log (DEpsilon)));
end Half_Log_Epsilon;
function Local_Atan
(Y : Real;
A : Real := 1.0)
return Real
is
Z : Real'Base;
Raw_Atan : Real'Base;
begin
if abs Y > abs A then
Z := abs (A / Y);
else
Z := abs (Y / A);
end if;
if Z < Square_Root_Epsilon then
Raw_Atan := Z;
elsif Z = 1.0 then
Raw_Atan := Pi / 4.0;
elsif Z < Square_Root_Epsilon then
Raw_Atan := Z;
else
Raw_Atan := Real'Base (Aux.Atan (Double (Z)));
end if;
if abs Y > abs A then
Raw_Atan := Half_Pi - Raw_Atan;
end if;
if A > 0.0 then
if Y > 0.0 then
return Raw_Atan;
else return -Raw_Atan;
end if;
else if Y > 0.0 then
return Pi - Raw_Atan;
else return -(Pi - Raw_Atan);
end if;
end if;
end Local_Atan;
function Log (A : Real) return Real is
begin
if A < 0.0 then
raise Argument_Error;
elsif A = 0.0 then
raise Constraint_Error;
elsif A = 1.0 then
return 0.0;
end if;
return Real (Aux.Log (Double (A)));
end Log;
function Log (A, Base : Real) return Real is
begin
if A < 0.0 then
raise Argument_Error;
elsif Base <= 0.0 or else Base = 1.0 then
raise Argument_Error;
elsif A = 0.0 then
raise Constraint_Error;
elsif A = 1.0 then
return 0.0;
end if;
return Real (Aux.Log (Double (A)) / Aux.Log (Double (Base)));
end Log;
function Log_Inverse_Epsilon return Real is
begin
return Real (Aux.Log (DIEpsilon));
end Log_Inverse_Epsilon;
function Sin (A : Real) return Real is
begin
if abs A < Square_Root_Epsilon then
return A;
end if;
return Real (Aux.Sin (Double (A)));
end Sin;
function Sin (A, Cycle : Real) return Real is
T : Real'Base;
begin
if Cycle <= 0.0 then
raise Argument_Error;
elsif A = 0.0 then
return A;
end if;
T := Exact_Remainder (A, Cycle) / Cycle;
if T = 0.0 or T = 0.5 or T = -0.5 then
return 0.0;
elsif T = 0.25 or T = -0.75 then
return 1.0;
elsif T = -0.25 or T = 0.75 then
return -1.0;
end if;
return Real (Aux.Sin (Double (T * Two_Pi)));
end Sin;
function Sinh (A : Real) return Real is
begin
if abs A < Square_Root_Epsilon then
return A;
elsif A > Log_Inverse_Epsilon then
return Exp (A - Log_Two);
elsif A < -Log_Inverse_Epsilon then
return -Exp ((-A) - Log_Two);
end if;
return Real (Aux.Sinh (Double (A)));
exception
when others =>
raise Constraint_Error;
end Sinh;
function Square_Root_Epsilon return Real is
begin
return Real (Aux.Sqrt (DEpsilon));
end Square_Root_Epsilon;
function Sqrt (A : Real) return Real is
begin
if A < 0.0 then
raise Argument_Error;
elsif A = 0.0 then
return A;
elsif A = 1.0 then
return 1.0;
end if;
return Real (Aux.Sqrt (Double (A)));
end Sqrt;
function Tan (A : Real) return Real is
begin
if abs A < Square_Root_Epsilon then
return A;
elsif abs A = Pi / 2.0 then
raise Constraint_Error;
end if;
return Real (Aux.tan (Double (A)));
end Tan;
function Tan (A, Cycle : Real) return Real is
T : Real'Base;
begin
if Cycle <= 0.0 then
raise Argument_Error;
elsif A = 0.0 then
return A;
end if;
T := Exact_Remainder (A, Cycle) / Cycle;
if T = 0.25
or else T = 0.75
or else T = -0.25
or else T = -0.75
then
raise Constraint_Error;
else
return Sin (T * Two_Pi) / Cos (T * Two_Pi);
end if;
end Tan;
function Tanh (A : Real) return Real is
begin
if A < Half_Log_Epsilon then
return -1.0;
elsif A > -Half_Log_Epsilon then
return 1.0;
elsif abs A < Square_Root_Epsilon then
return A;
end if;
return Real (Aux.tanh (Double (A)));
end Tanh;
function LOG10 (A : REAL) return REAL is
begin
return Log (A, 10.0);
end LOG10;
function LOG2 (A : REAL) return REAL is
begin
return Log (A, 2.0);
end LOG2;
function ASIN (A : REAL) return REAL renames Arcsin;
function ACOS (A : REAL) return REAL renames Arccos;
function ATAN (A : REAL) return REAL is
begin
return Arctan (A, 1.0);
end ATAN;
function ATAN2 (A1, A2 : REAL) return REAL renames Arctan;
function SIND (A : REAL) return REAL is
begin
return Sin (A, 360.0);
end SIND;
function COSD (A : REAL) return REAL is
begin
return Cos (A, 360.0);
end COSD;
function TAND (A : REAL) return REAL is
begin
return Tan (A, 360.0);
end TAND;
function ASIND (A : REAL) return REAL is
begin
return Arcsin (A, 360.0);
end ASIND;
function ACOSD (A : REAL) return REAL is
begin
return Arccos (A, 360.0);
end ACOSD;
function Arctan (A : REAL) return REAL is
begin
return Arctan (A, 1.0, 360.0);
end Arctan;
function ATAND (A : REAL) return REAL is
begin
return Arctan (A, 1.0, 360.0);
end ATAND;
function ATAN2D (A1, A2 : REAL) return REAL is
begin
return Arctan (A1, A2, 360.0);
end ATAN2D;
end Math_Lib;