package body MATHS is -- based on Abramowitz & STEGUN 1965
-- Accuracy is limitted to 4 or 5 significant decimal figures.
	--MATH_ERROR,MISSING_FUNCTION:exception;
	--function SQRT( X:FLOAT ) return FLOAT;--X>=0. SQRT(X)*SQRT(X)=X.
	--function LOG(BASE:FLOAT; X:FLOAT) return FLOAT;
	--function LN(X:FLOAT)return FLOAT;--LOG(BASE=>e, X)
	--function LG(X:FLOAT)return FLOAT;--LOG(BASE=>2.0,X)
	--function EXP(X:FLOAT)return FLOAT;--e**X
	--function COS(X:RADIANS)return FLOAT;--cosine of X in radians
	--function SIN(X:RADIANS)return FLOAT;--sine of X in radians
	--type RADIANS, DEGREES are new FLOAT, "+" converts between them
 
function LN(X:FLOAT)return FLOAT is --Natural logariethm - base=>e
	EXPONENT:FLOAT:=0.0;
	X1:FLOAT:=X;
begin
	if X <0.0 then raise MATH_ERROR;
	end if;
	while X1 > 2.0 loop
		X1:=X1/E; EXPONENT:=EXPONENT+1.0;
	end loop; -- X1<2.0 and LN(X)=EXPONENT+LN(X1)
 
	X1:=X1-1.0;
 
	return EXPONENT+X1*(0.99949_566+
			  X1*(-0.49190_896+
			   X1*(0.28947_478+
			    X1*(-0.13606_275+
			      X1*0.03215_845
                           ))));
 
end LN;
 
function LOG(X:FLOAT) return FLOAT is
begin
	return LOG(10.0, X);
end LOG;
 
function LOG(BASE:FLOAT; X:FLOAT) return FLOAT is 
begin
	if X <0.0 then raise MATH_ERROR;
	end if;
	return LN(X)/LN(BASE);
 
end LOG;
 
function SQRT( X:FLOAT ) return FLOAT is --X>=0. SQRT(X)*SQRT(X)=X.
begin
	if X <0.0 then raise MATH_ERROR;
	end if;
	return EXP(0.5*LN(X));
 
end SQRT;
 
function EXP(X:FLOAT)return FLOAT is --e**X,inverse of ln(X)
	POWER:FLOAT; X1:FLOAT:=X;
begin	POWER:=E**INTEGER(X1); X1:=X1-FLOAT(INTEGER(X1));
	return POWER* (1.0+X1*
			(0.99986_84+X1*
			 (0.49829_26+X1*
			  (0.15953_32+X1*0.0293641
		       ))));
end EXP;
 
function LG(X:FLOAT)return FLOAT is
begin
	if X <0.0 then raise MATH_ERROR;
	end if;
	return LOG(2.0, X);
 
end LG;
 
function "+"(X:RADIANS)return DEGREES is
begin
	return DEGREES(180.0*FLOAT(X)/PI);
end "+";
 
function "+"(X:DEGREES)return RADIANS is
begin
	return RADIANS(PI*FLOAT(X)/180.0);
end "+";
 
 
function COS(X:RADIANS)return FLOAT is --cosine of X in radians
begin	
	return SIN(PI*0.5-X);
 
end COS;
 
function SIN(X:RADIANS)return FLOAT is --sine of X in radians 
	T:FLOAT:=FLOAT(X);
	ERROR_TERM:constant:=1.0E-7;
	S:FLOAT:=0.0;
	I:INTEGER:=2;
begin	if X = 0.0 or X=PI then return 0.0;
	elsif X=PI*0.5 then return 1.0;
	elsif X<0.0 then return -SIN(-X);
	elsif X>PI2 then 
		while T >PI2 loop T:=T-PI2;
		end loop;
		return SIN(RADIANS(T));
	elsif X>PI then return -SIN(X-PI);
	elsif X>PI*0.5 then return SIN(PI-X);
		while abs(T)>ERROR_TERM loop
			S:=S+FLOAT(T); 
			T:=-T*FLOAT(X*X)/FLOAT(I*(I+1));
			I:=I+2;
		end loop;
		return S;
	else raise MISSING_FUNCTION;
	end if;
end SIN;
 
end MATHS;