! Trig package
! JWM August 1982
option "-nocheck-nodiag-noline"
const real PIO4 = 0.7853981633974483096156608
const real PIO2 = 1.570796326794896619231322
const real PI = 3.141592653589793238462643
const real PI2 = 6.283185307179586476925287
const real E = 2.718281828
real function NORMALISE(real x)
{brings result into range -pi <= result <= +pi}
x = x-int pt(x/pi2)*pi2
result = pi2+x if x <= -pi
result = x-pi2 if x >= pi
result = x
end {of normalise}
real function RAW SINE(real x)
{parameter must be between 0 and pi/4}
real y; y = x*x
result = (((0.0032811761*y-0.1335639326)*y+1.0)*x)/ c
((0.0004649838*y+0.0331027317)*y+1.0)
end {of raw sine}
real function RAW COSINE(real x)
{parameter must be between 0 and pi/4}
real y; y = x*x
result = ((0.0205121130*y-0.4558922221)*y+0.9999999992)/ c
((0.0008996261*y+0.0441077396)*y+1.0)
end {of raw cosine}
external real function SIN(real x)
result = -sin(-x) if x < 0
x = normalise(x)
x = pi-x if x > pio2
result = raw cosine(pio2-x) if x > pio4
result = raw sine(x)
end {of sin}
external real function COS(real x)
x = normalise(x)
x = -x if x < 0
if x > pio2 start
x = pi-x
result = -raw sine(pio2-x) if x > pio4
result = -raw cosine(x)
finish else start
result = raw sine(pio2-x) if x > pio4
result = raw cosine(x)
finish
end {of cos}
external real function TAN(real x)
{%real s
{ s = sin(x)
{ %result = s/sqrt(1.0-s*s)
result = sin(x)/cos(x)
end {of tan}
real function A TAN(real x)
! parameter must be in range 0 <= x <= 1
real y; y = x*x
result = c
x*(((0.0089472229*y+0.2870044785)*y+1.1303754276)*y+0.9999999992)/ c
(((0.0506770959*y+0.5749098994)*y+1.4637086946)*y+1.0)
end
external real function P ARC TAN(real x)
result = -p arc tan(-x) if x < 0
result = a tan(x) if 0 <= x <= 1.0
result = pio2 - a tan(1.0/x)
end {of arc tan1}
external real function ARC TAN(real x,y)
! Returns the value, in radians, of the angle whose tangent
! is specified by "y/x". This value is between +&- pi. If
! x<0 then it is in the second or third quadrants otherwise
! it is in the first or fourth quadrants.
real z
z = p arc tan(y/x) {don't have to worry about -ve case}
result = z if x >= 0
result = pi-z if z > 0
result = -pi-z
end {of arc tan}
external real function ARC SIN(real x)
if x >= 1.0 start
signal 10 if x > 1.0
result = pio2
finish
if x <= -1.0 start
signal 10 if x < 1.0
result = -pio2
finish
result = p arc tan(x/(sqrt(1.0-x*x)))
end {of arc sin}
external real function ARC COS(real x)
real w
if x > 0.999 start
signal 10 if x > 1.0
result = 0
finish
if x < -0.999 start
signal 10 if x < -1.0
result = pi
finish
result = pio2 if -0.001 < x < 0.001
w = p arc tan(sqrt(1.0/(x*x)-1.0))
{this next because arc tan is in range between +&-pi}
result = w if x >= 0
result = -w+pi
end {of arc cos}
real function LOG E(real z)
! z is such that 0.5 <= z <= 1
real x,y
x = z+z-1.0 {nasty VAX compiler bug}
y = x*x
result = ((0.0956558162*x+0.5297501385)*y+0.0677412133*x-0.6931471773) c
/((0.0286818192*x+0.45477291277)*y+1.3449644663*x+1.0)
end {of log e}
external real function LOG(real x)
! logarithm to base e
const real log of two = 0.6931471806
integer factor; factor = 0
signal 10 if x <= 0
if x > 1.0 start
cycle
factor = factor+1; x = x/2
repeat until 0.5 <= x <= 1.0
result = log e(x) + factor*log of two
finish else start
result = log e(x) if x >= 0.5
cycle
factor = factor+1; x = x+x
repeat until 0.5 <= x <= 1.0
result = log e(x) - factor*log of two
finish
end {of log}
external real function LOG TEN(real x)
! logarithm to base 10
result = log(x)/log(10.0)
end {of log10}
real function RAW EXP(real x)
! x must be in the range 0 <= x <= 1 (to give e to the power x)
real y; y = x*x
result = ((0.0106337905*x+0.1125548636)*y+0.5240642207*x+1.0)/ c
((0.0884921370-0.0065658101*x)*y+1.0000000007-0.4759358618*x)
end {of raw exp}
external real function EXP(real x)
! e to the power x
real temp
if -1.0 <= x <= 1.0 start
result = raw exp(x) if x >= 0
result = 1.0/raw exp(-x)
finish else start
temp = exp(0.5*x); result = temp*temp
finish
end {of exp}
external real function SIN H(real x)
result = 0.5*(exp(x)-exp(-x))
end {of sin h}
external real function COS H(real x)
result = 0.5*(exp(x)+exp(-x))
end {of cos h}
external real function TAN H(real x)
result = sin h(x)/cos h(x)
end {of tan h}
end of file