program mathpak;
{(C) Australian Computer Society 1984}
{based on paper A Compact Mathematical Function Package
by A.P. Clarke and W. Marwood,
The Australian Computer Journal,
Vol 16, No 3, August 1984, pp 107 to 114
entered by W F McGee 1 Dec 1984}
const e=1.0e4;
sqrt_pi=1.772453850905;
sqrt2 =1.414213562373;
eul =0.577215664901;
sqrt12 =3.464101615137;
ln2 =0.693147180559;
ln10 =2.302585092994;
ln17 =2.833213344056;
TERMMAX=500;
type available=(sif,cif,e1f,jnf,mbif,lagf,qf,phif,erff,erfcf,chgf2f,sxf,cxf,qmnf,
pmnf,hf,gmf,hgff,chgf1f,ghgff,cf,sf,expf,lf,l10f,done);
string25=string[25];
var prec:real;
numparm:array[available]of integer;
name:array[available]of string25;
sig_dig:integer;
sig_dig_flag:boolean;
sig_dig_real,Ubound,lngmb,gmb,mc_dig:real;
x:real;
{**********************************************************}
function lnn(x:real):real;
{natural log with zero exception}
begin
if x=0 then lnn:=-86 else lnn:=ln(x)
{IMPLEMENTATION DEPENDENT 38*ln10}
end;
{***********************************************************}
procedure sum_error(var1,dig1,var2,dig2:real);
var error:real;
begin
error:=abs(var1)*(exp(-dig1*ln10)+prec) +
abs(var2)*(exp(-dig2*ln10)+prec);
sig_dig:=trunc(ln(abs((var1+var2)/error))/ln10)
end;
{***********************************************************}
function ghgf(p,q:integer;a1,a2,a3,a4,c1,c2,c3,c4,x:real):real;
const pqmax=4;
type series=array[0..TERMMAX] of real;
poch=array[1..pqmax]of real;
var a,c:poch;
sum,y,abs_y,abs_last_y,sigma:real;
pqm,pq,k,i,j:integer;
convergent:boolean;
ser:series;

procedure accuracy;
function sig(terms:integer):real;
var i:integer;psum,s1,s2,weight:real;
begin
psum:=0; s1:=sum; s2:=ser[1];weight:=2*(p+q+1);
for i:=1 to terms-1 do begin
s1:=s1-ser[i-1];s2:=s2+ser[i+1];
psum:=psum+sqr(s1)+weight*sqr(s2)
end;
sig:=sqrt(psum+sqr(s1)+weight*sqr(ser[terms]))
end;

begin
sigma:=prec/sqrt12*sig(k-1);
if sigma=0 then sig_dig_real:=mc_dig
else if sum=0 then sig_dig_real:=0
else sig_dig_real:=ln(abs(sum/sigma))/ln10;
if sig_dig_real<=0.0 then sig_dig_real:=0.0;
if sig_dig_real> mc_dig then sig_dig_real:=mc_dig;
sig_dig:=trunc(sig_dig_real)
end;
begin
a[1]:=a1; a[2]:=a2; a[3]:=a3; a[4]:=a4;
c[1]:=c1; c[2]:=c2; c[3]:=c3; c[4]:=c4;
if p>q then pqm:=p else pqm:=q;
y:=1.0; sum:=1.0; k:=1; ser[0]:=1.0; pq:=p+q;
convergent:=false; abs_last_y:=1.0;
repeat
for i:=1 to pqm do y:=y*a[i]/c[i];
for i:=1 to p do a[i]:=a[i]+1;
for i:=1 to q do c[i]:=c[i]+1;
y:=y*x/k; abs_y:=abs(y);
if (abs_y < abs_last_y) then convergent:=true;
if (abs_y > abs_last_y) and convergent then y:=0;
sum:=sum+y; abs_last_y:=abs(y);ser[k]:=y; k:=k+1
until (sum=sum-y) or (k>499);
if (k>499) then writeln('not convergent after 499 iterations');
if sig_dig_flag then accuracy;
ghgf:=sum
end;

procedure init_consts;
begin
sig_dig:=0;Ubound:=10; gmb:=17; lngmb:=ln17;
mc_dig:=ln(1/prec)/ln10
end;

function chgf1(a,c,x:real):real;
{confluent hypergeometric function of the first kind}
begin
if x<0.0 then chgf1:=exp(x)*chgf1(c-a,c,-x)
else chgf1:=ghgf(1,1,a,1,1,1,c,1,1,1,x)
end;

function si(x:real):real;
{sine integral}
begin
si:=x*ghgf(1,2,0.5,1,1,1,1.5,1.5,1,1,-x*x/4.0)
end;

function e1(x:real):real;
{Exponential integral}
var var1,var2:real;
begin
var1:=x*ghgf(2,2,1,1,1,1,2,2,1,1,-x); var2:=-ln(x);
e1:=-eul+var1+var2;
if sig_dig_flag then begin
sum_error(-eul,mc_dig,var1,sig_dig_real);
sum_error(-eul+var1,sig_dig_real,var2,mc_dig)
end
end;

function ci(x:real):real;
{Cosine integral}
var var1,var2:real;
begin
var1:=-sqr(x/2)*ghgf(2,3,1,1,1,1,1.5,2,2,1,-sqr(x/2)); var2:=ln(x);
ci:=eul+var1+var2;
if sig_dig_flag then begin
sum_error(eul,mc_dig,var1,sig_dig_real);
sum_error(eul+var1,sig_dig_real,var2,mc_dig)
end
end;

function hgf(a1,a2,c,x:real):real;
{Gauss hypergeometric function}
begin
hgf:=ghgf(2,1,a1,a2,1,1,c,1,1,1,x)
end;

function gm(x:real):real;
{Gamma function}
var temp_flag:boolean;
begin
temp_flag:=sig_dig_flag; sig_dig_flag:=false;
if x=0.0 then gm:=1.0
else gm:=exp(x*lngmb)/x*chgf1(x,x+1,-gmb)+
exp((x-1)*lngmb-gmb)*ghgf(2,0,1-x,1,1,1,1,1,1,1,-1.0/gmb);
sig_dig_flag:=temp_flag
end;

function h(n:integer;x:real):real;
{Hermite polynomials}
var sign:integer;
begin
if odd(n div 2) then sign:=-1 else sign:=1;
if (n mod 2)=0 then h:=sign*gm(n+1)/gm(n/2+1)*chgf1(-n/2,0.5,x*x)
else h:=sign*gm(n+1)/gm((n+1)/2)*2*x*chgf1(-(n-1)/2,1.5,x*x)
end;

function pmn(m,n:integer;x:real):real;
{Associated Legendre polynomial of the first kind}
var sign:real;
begin
if abs(x)>1 then begin
writeln('usage: pmn needs abs(argument)<1');
pmn:=0.0
end
else begin
if odd(m) then sign:=-1.0 else sign:=1.0;
pmn:=gm(m+n+1)/(gm(m+1)*gm(n-m+1))*exp(ln(1-sqr(x))*m/2-ln2*m)*
sign*ghgf(2,1,m-n,m+n+1,1,1,1+m,1,1,1,(1-x)/2)
end
end;

function qmn(m,n:integer; x:real):real;
{Associated Legendre polynomial of the second kind}
var sign:real;
temp1:real;
begin
if abs(x)<1 then begin
writeln('usage:qmn needs abs(argument) >1');
qmn:=0.0
end
else begin
if odd(m) then sign:=-1.0 else sign:=1.0;
temp1:=(lnn(x*x-1.0)*m/2.0)-lnn(x)*(n+m+1.0);
qmn:=exp(temp1-ln2*(n+1))*sqrt_pi*
gm(n+m+1)/gm(n+1.5)*sign*
ghgf(2,1,(1+n+m)/2,(2+m+n)/2,1,1,n+1.5,1,1,1,1/(x*x))
end
end;

function cx(x:real):real;
{Fresnel integral of the first kind}
begin
cx:=x*ghgf(1,2,0.25,1,1,1,0.5,1.25,1,1,-sqr(pi*sqr(x/2)))
end;

function sx(x:real):real;
{Fresnel integral of the second kind}
begin
sx:=pi*x*x*x/6*ghgf(1,2,0.75,1,1,1,1.5,1.75,1,1,-sqr(pi*sqr(x/2)))
end;

function chgf2(a,b,x:real):real;
{confluent hypergeometric function of the second kind}
var temp1,temp2,sig1:real;
begin
if abs(x) temp1:=chgf1(a,b,x)/gm(1+a-b)/gm(b);sig1:=sig_dig_real;
temp2:=chgf1(1+a-b,2-b,x)*exp((1-b)*lnn(x))/gm(a)/gm(2-b);
chgf2:=(temp1-temp2)*pi/sin(pi*b);
if sig_dig_flag then sum_error(temp1,sig1,-temp2,sig_dig_real)
end
else chgf2:=exp(-a*lnn(x))*ghgf(2,0,a,1+a-b,1,1,1,1,1,1,-1.0/x)
{Asymptotic representation for large abs(x)}
end;

function erfc(x:real):real;
{Complementary error function}
begin
if x<0 then erfc:=2.0-exp(-x*x)*chgf2(0.5,0.5,x*x)/sqrt_pi
else erfc:=exp(-x*x)*chgf2(0.5,0.5,x*x)/sqrt_pi
end;

function erf(z:real):real;
{Error function}
begin
if abs(z)<3.0 then erf:=z*chgf1(0.5,1.5,-z*z)*2/sqrt_pi
else erf:=1.0-erfc(z)
end;

function phi(z:real):real;
{Normal probability integral}
var var1:real;
begin
var1:=erf(z/sqrt2);phi:=0.5+0.5*var1;
if sig_dig_flag then sum_error(0.5,mc_dig,0.5*var1,sig_dig_real)
end;

function q(x:real):real;
{Complementary Normal probability integral}
begin
q:=0.5*erfc(x/sqrt2)
end;

function lag(alpha,n,z:real):real;
{Laguerre polynomials}
begin
lag:=gm(alpha+n+1)/(gm(n+1)*gm(alpha+1))*chgf1(-n,alpha+1,z)
end;

function mbi(nn,x:real):real;
{Modified Bessel function}
var n:real;
begin
n:=nn;
if round(n)=n then n:=abs(n);
mbi:=exp(n*lnn(x/2))/gm(n+1)*ghgf(0,1,1,1,1,1,1+n,1,1,1,sqr(x/2))
end;

function jn(n,x:real):real;
{Bessel function}
begin
jn:=exp(n*lnn(x/2))/gm(n+1)*ghgf(0,1,1,1,1,1,1+n,1,1,1,-sqr(x/2))
end;

procedure initfunctions;
begin
name[sif]:='Sine-Integral';
name[cif]:='Cosine-Integral';
name[e1f]:='Exponential-Integral';
name[jnf]:='Bessel';
name[mbif]:='Modified-Bessel';
name[lagf]:='Laguerre';
name[qf]:='Complementary-Normal';
name[phif]:='Normal-Probability';
name[erff]:='Error';
name[erfcf]:='Complementary-Error';
name[chgf2f]:='Confluent-2nd-type';
name[sxf]:='2nd-Fresnel-Integral';
name[cxf]:='1st-Fresnel-Integral';
name[qmnf]:='2nd-Type-Legendre-(Assc)';
name[pmnf]:='1st-Type-Legendre-(Assc)';
name[hf]:='Hermite';
name[gmf]:='Gamma';
name[hgff]:='Gauss-Hypergeometric';
name[chgf1f]:='Confluent-1st-Type';
name[ghgff]:='Generalized-Hyper';
name[cf]:='Cosine';
name[sf]:='Sine';
name[lf]:='Ln';
name[l10f]:='Log10';
name[expf]:='Exp';
name[done]:='DONE';

numparm[sif]:=1;
numparm[cif]:=1;
numparm[e1f]:=1;
numparm[jnf]:=2;
numparm[mbif]:=2;
numparm[lagf]:=3;
numparm[qf]:=1;
numparm[phif]:=1;
numparm[erff]:=1;
numparm[erfcf]:=1;
numparm[chgf2f]:=3;
numparm[sxf]:=1;
numparm[cxf]:=1;
numparm[qmnf]:=3;
numparm[pmnf]:=3;
numparm[hf]:=2;
numparm[gmf]:=1;
numparm[hgff]:=4;
numparm[chgf1f]:=3;
numparm[ghgff]:=11;
numparm[cf]:=1;
numparm[sf]:=1;
numparm[lf]:=1;
numparm[l10f]:=1;
numparm[expf]:=1;
numparm[done]:=0;
end;

function calculate:available;
var i:integer;r:real;index:available;
getstring:string25;
vari:array[1..11]of real;
procedure getvar(num:integer);
var i:integer;
begin
writeln('Entry of parameters');
if num=0 then
else for i:=1 to num do begin
write('Parameter [ ',i,'] is ');readln(vari[i])
end
end;
begin
index:=sif;
repeat
write(name[index]:30);
index:=succ(index);
if index<>done then begin
writeln(name[index]:30);index:=succ(index)
end
until index=done;
writeln;writeln('Enter function desired');
readln(getstring);
index:=sif;
while not ((index=done)or(name[index]=getstring)) do index:=succ(index);
getvar(numparm[index]);
writeln;writeln('Working....');
if index=sif then r:=si(vari[1]);
if index=cif then r:=ci(vari[1]);
if index=e1f then r:=e1(vari[1]);
if index=jnf then r:=jn(round(vari[1]),vari[2]);
if index=mbif then r:=mbi(round(vari[1]),vari[2]);
if index=lagf then r:=lag(vari[1],vari[2],vari[3]);
if index=qf then r:=q(vari[1]);
if index=phif then r:=phi(vari[1]);
if index=erff then r:=erf(vari[1]);
if index=erfcf then r:=erfc(vari[1]);
if index=chgf2f then r:=chgf2(vari[1],vari[2],vari[3]);
if index=sxf then r:=sx(vari[1]);
if index=cxf then r:=cx(vari[1]);
if index=qmnf then r:=qmn(round(vari[1]),round(vari[2]),vari[3]);
if index=pmnf then r:=pmn(round(vari[1]),round(vari[2]),vari[3]);
if index=hf then r:=h(round(vari[1]),vari[2]);
if index=gmf then r:=gm(vari[1]);
if index=hgff then r:=hgf(vari[1],vari[2],vari[3],vari[4]);
if index=chgf1f then r:=chgf1(vari[1],vari[2],vari[3]);
if index=ghgff then r:=ghgf(round(vari[1]),round(vari[2]),vari[3],vari[4],
vari[5],vari[6],vari[7],vari[8],vari[9],
vari[10],vari[11]);
if index=cf then r:=cos(vari[1]);
if index=sf then r:=sin(vari[1]);
if index=lf then r:=ln(vari[1]);
if index=expf then r:=exp(vari[1]);
if index=l10f then r:=ln(vari[1])/ln10;
if index=done then r:=0.0;
clrscr;
write(name[index],'[');
if numparm[index]<>0 then begin
for i:=1 to numparm[index] do begin
write(vari[i]:12:6);if i<>numparm[index] then write(' , ')
end
end;
write(' ]= ');writeln(r);
writeln;writeln('significant digits = ',sig_dig);
calculate:=index
end;

function precision:real;
var x,y:real;
begin
y:=1.0;
x:=1.0;
repeat
x:=x+y;
y:=10.0*y
until x+1.0=x;
x:=y/100.0;
repeat
x:=x+y/1000.0
until x+1.0=x;
x:=x-y/1000;
precision:=1.0/x
end;
begin
clrscr;
writeln('FUNCTION PACKAGE BY CLARKE AND MARWOOD');
writeln('for TURBO PASCAL peak exponent 10E(+/- 38)');
writeln('ignore digits precision for sine,cos,exp,ln,log10');
prec:=precision;
writeln('PRECISION =',prec);
init_consts;
sig_dig_flag:=true;
initfunctions;
repeat
until calculate=done
end.
function Q(x:real):real;
var temp,factor,con,z,q1:real;
begin
con:=2.506628274;
if abs(x)>12.5 then q1:=0.0
else begin
temp:=sqr(x);
z:=exp(-temp/2.0);
factor:=abs(x/2)+sqrt(1.0+temp/4.0);
q1:=z/(con*factor+(2.0-con)*z)
end;
if x>=0.0 then Q:=q1 else Q:=1.0-q1;
end;

function arcq(qd:real):real;
var x0,x:real;
begin
if abs(qd)>=1.0 then writeln('out of range')
else begin
if qd>0.5 then x0:=0 else x0:=sqrt(-2*ln(qd))-1.0;
repeat
x:=x0;
x0:=x0-(q(x0)-qd)/(-exp(-sqr(x0)/2)/sqrt(2*pi));
until abs(x0-x)<0.00000001
end;
arcq:=x0
end;
 x:=x0;
x0:=x0-(q(x0)-qd)/