!This is an old ERCC IMP library (advisory have documentation). It contains
! a variety of useful procedures
! IMP DFLEPO IMP DGAHER IMP DGALAG IMP DGALEG IMP DIGVEC IMP DKUSCR
! IMP DSCRCH IMP FLEPO IMP GAHER IMP GALAG IMP GALEG IMP LOGFCT
! IMP DAPRO IMP DCHLIN IMP DELI3 IMP PLEX IMP DPLEX IMP VPLEX
! IMP DVPLEX IMP DEXPAN IMP DERFIN IMP CHINT IMP DCHINT IMP DKUMER
! IMP DAVDEN IMP DEXPFT IMP DIVINT IMP DCHCOL IMP DRUNGE IMP DRUMOD
! IMP DBUTCH IMP DNORSK IMP DGROTA IMP DLINPR IMP INFNT IMP DINFNT
! IMP GRAPH9 IMP DCURFT IMP DROOTS IMP DSQRTM IMP DIFSYS IMP CDEVLR
! IMP DRATNT IMP ELI3 IMP SELECT IMP FILON IMP DFILON IMP DLEQSP
! IMP QADMN IMP DETSP IMP INSORT IMP RLSORT IMP CDHESS
!*TOP*IMPDFLEPO
external routine imp dflepo(routine funct(long real name f, long real array name x, g), integer n, limit, iprin,
integer name conv, long real eps, long real name f,
long real array name x, h)
long real oldf, sg, ghg, z, step, ita, fa, fb, ga, gb, w, lambda
integer i, j, k, count, fnevalns
long real array g, s, gamma, sigma, y(1:n)
long real function dot(long real array name a, b)
integer i
long real s
s = 0
cycle i = 1, 1, n
s = s+a(i)*b(i)
repeat
result = s
end
long real function updot(integer longi, long real array name a, b)
comment MULTIPLY BY THE ITH ROW OF THE SYMMETRIC MATRIX A, c
!WHOSE UPPER TRIANGLE IS STORED BY ROWS.
integer j, k; long real s
k = i; s = 0
->nrl1 if i=1
cycle j = 1, 1, i-1
s = s+a(k)*b(j)
k = k+n-j
repeat; !STEPS TO DIAGONAL:NOW GO ALONG ROW.
nrl1: cycle j = i, 1, n
s = s+a(k+j-i)*b(j)
repeat
result = s
end; !OF FN UPDOT
nrl1:
->nrl2 unless conv=1
k = 1
cycle i = 1, 1, n
h(k) = 1;
->nrl30 if i=n
cycle j = 1, 1, n-i
h(k+j) = 0
repeat
nrl30:k = k+n-i+1
repeat; !FORMATION OF UNIT MATRIX IN H.
nrl2:conv = 1; step = 1; fnevalns = 0
funct(f, x, g); fnevalns = fnevalns+1
oldf = f+1
cycle count = 1, 1, limit
return unless oldf>f
cycle i = 1, 1, n
sigma(i) = x(i); gamma(i) = g(i); s(i) = -updot(i, h, g)
repeat
fb = f; gb = dot(g, s)
->nrl21 unless gb>=0
conv = -1; return
nrl21:oldf = f; ita = step
! ITA IS DETERMINED FROM ITS VALUE AT THE PREVIOUS ITERATION (STEP)
! AND IS INCREASED OR DECREASED BY 4 WHERE NECESSARY. IT SHOULD
! TEND TO 1 AT THE MINIMUM.
nrl4: fa = fb; ga = gb
cycle i = 1, 1, n; x(i) = x(i)+ita*s(i); repeat
funct(f, x, g); fb = f; gb = dot(g, s); fnevalns = fnevalns+1
->nrl3 unless (gb<0 and fb<fa)
ita = 4*ita; step = 4*step;
->nrl4
nrl3: if ita<=0 then start; conv = -1; return; finish
cycle i = 1, 1, n; y(i) = g(i); repeat
z = 3*(fa-fb)/ita+ga+gb
w = sqrt(z*z-ga*gb)
lambda = ita*(gb+w-z)/(gb-ga+2*w)
cycle i = 1, 1, n; x(i) = x(i)-lambda*s(i); repeat
funct(f, x, g); fnevalns = fnevalns+1
->nrl6 unless (f>fa or f>fb)
step = step/4
->nrl5 unless fb<fa
cycle i = 1, 1, n; x(i) = x(i)+lambda*s(i); g(i) = y(i)
repeat
f = fb;
->nrl6
nrl5: gb = dot(g, s)
->nrl7 unless (gb<0 and count>n and step<9.999999999999971@-7)
conv = -1; return
nrl7: fb = f; ita = ita-lambda;
->nrl3
nrl6: cycle i = 1, 1, n
sigma(i) = x(i)-sigma(i)
gamma(i) = g(i)-gamma(i)
repeat
sg = dot(sigma, gamma)
->nrl13 unless count>=n
return if sqrt(dot(s, s))<eps and sqrt(dot(sigma, sigma))<eps
comment end OFTESTFORVANISHINGDERIVATIVE.
nrl13:cycle i = 1, 1, n
s(i) = updot(i, h, gamma)
repeat
ghg = dot(s, gamma)
k = 1
->nrl8 if mod(sg)<9.999999999999987@-31 or c
mod(ghg)<9.999999999999987@-31
cycle i = 1, 1, n
cycle j = i, 1, n
h(k) = h(k)+sigma(i)*sigma(j)/sg-s(i)*s(j)/ghg
k = k+1; repeat
repeat
->nrl20 if iprin=0
->nrl20 unless count//iprin*iprin=count
nrl8: newlines(2);
printstring("
ITERATION ")
print(count, 2, 0)
printstring(":
NO. OF FN EVALNS = ")
write(fnevalns, 2)
newlines(2);
printstring(" X G
")
cycle i = 1, 1, n
printfl(x(i), 6); spaces(4); printfl(g(i), 6); newline
repeat;
printstring("
PRESENT MINIMUM VALUE OF FN =")
printfl(f, 6)
newlines(2)
nrl20:repeat; ! COUNT LOOP; CONV=-1 ;
printstring("
PROBLEM DID NOT CONVERGE IN SPECIFIED NO OF ITERATIONS
ROUTINE TERMINATED.")
newlines(2)
end; !OF ROUTINE IMP DFLEPO
!*END*IMPDFLEPO
!*TOP*IMPDGAHER
external routine imp dgaher(long real function f(long real x), long real eps,
long real name antgrl, test)
long real i4, i6, i9, i12, i15, c
i4 = .813128354472@-1*(f(1.65068012388)+f(-1.65068012388))+.804914090005* c
(f(.524647623275)+f(-.524647623275))
i6 = .453000990551@-2*(f(2.35060497367)+f(-2.35060497367))+.157067320323* c
(f(1.33584907401)+f(-1.33584907401))+.724629595224*(f(.436077411928)+f c
(-.436077411928))
c = mod(mod(i6)-mod(i4))
if (mod(i6))>1 then c = c/(mod(i6))
->nrl1 unless c<eps
antgrl = i6
test = c
return
nrl1:i9 = .396069772632@-4*(f(3.19099320178)+f(-3.19099320178))+ c
.494362427553-2*(f(2.26658058453)+f(-2.26658058453))
i9 = i9+.884745273943-1*(f(1.46855328921)+f(-1.46855328921))+ c
.432651559002*(f(.723551018752)+f(-.723551018752))+.720235215606*f(0)
c = mod(mod(i9)-mod(i6))
if (mod(i9))>1 then c = c/(mod(i9))
->nrl2 unless c<eps
antgrl = i9
test = c
return
nrl2:i12 = .265855168436@-6*(f(3.88972489787)+f(-3.88972489787))+ c
.857368704358@-4*(f(3.02063702512)+f(-3.02063702512))+.390539058463@-2* c
(f(2.27950708050)+f(-2.27950708050))
i12 = i12+.516079856158@-1*(f(1.59768263515)+f(-1.59768263515))+ c
.260492310264*(f(.947788391240)+f(-.947788391240))+.570135236262*(f( c
.314240376254)+f(-.314240376254))
c = mod(mod(i12)-mod(i9))
if (mod(i12))>1 then c = c/(mod(i12))
->nrl3 unless c<eps
antgrl = i12
test = c
return
nrl3:i15 = .152247580425@-8*(f(4.49999070731)+f(-4.49999070731))+ c
.105911554771@-5*(f(3.66995037340)+f(-3.66995037340))
i15 = i15+.100004441232@-3*(f(2.96716692791)+f(-2.96716692791))+ c
.277806884291@-2*(f(2.32573248617)+f(-2.32573248617))
i15 = i15+.307800338725@-1*(f(1.71999257519)+f(-1.71999257519))+ c
.158488915796*(f(1.13611558521)+f(-1.13611558521))
i15 = i15+.412028687499*(f(.565069583256)+f(-.565069583256))+ c
.564100308726*f(0)
c = mod(mod(i15)-mod(i12))
if (mod(i15))>1 then c = c/(mod(i15))
antgrl = i15
test = c
return if c<eps
printstring("
METHOD FAILED TO SATISFY GIVEN ACCURACY. LAST")
printstring(" ESTIMATE OF INTEGRAL IS ")
printfl(antgrl, 10)
printstring("
THE ERROR OF THIS ANSWER IS ")
printfl(test, 10)
return
end; !OF IMP DGAHER
!*END*IMPDGAHER
!*TOP*IMPDGALAG
external routine imp dgalag(long real function f(long real x), long real eps,
long real name antgrl, test)
long real s4, s6, s9, s12, s15, c
s4 = 6.03154104342@-1*f(0.322547689619)+3.57418692438@-1*f(1.74576110116) c
+3.88879085150@-2*f(4.53662029692)+5.39294705561@-4*f(9.395070912301)
s6 = 4.58964673950@-1*f(.222846604179)+4.17000830772@-1*f(1.18893210167)
s6 = s6+1.13373382074@-1*f(2.99273632606)+1.03991974531@-2*f(5 c
.77514356910)+2.61017202815@-4*f(9.83746741838)+8.98547906430@-7*f(15 c
.9828739806)
c = mod(mod(s6)-mod(s4))
if (mod(s6))>1 then c = c/(mod(s6))
->nrl1 unless c<eps
antgrl = s6
test = c
return
nrl1:s9 = 3.36126421798@-1*f(.152322227732)+4.11213980424@-1*f(.807220022742)
s9 = s9+1.99287525371@-1*f(2.00513515562)+4.74605627657@-2*f(3.78347397333)
s9 = s9+5.59962661079@-3*f(6.20495677788)+3.05249767093@-4*f(9.37298525169)
s9 = s9+6.59212302608@-6*f(13.4662369111)+4.11076933035@-8*f(18.8335977890)
s9 = s9+3.29087403035@-11*f(26.3740718909)
c = mod(mod(s9)-mod(s6))
if (mod(s9))>1 then c = c/(mod(s9))
->nrl2 unless c<eps
antgrl = s9
test = c
return
nrl2:s12 = 2.64731371055@-1*f(.115722117358)+3.77759275873@-1*f(.61175748452) c
+9.04492222117@-2*f(2.83375133774)
s12 = s12+2.44082011320@-1*f(1.51261026977)
s12 = s12+2.01023811546@-2*f(4.59922763942)+2.66397354187@-3*f(6 c
.84452545312)
s12 = s12+2.03231592663@-4*f(9.62131684246)+8.36505585682@-6*f(13 c
.0060549933)
s12 = s12+1.66849387654@-7*f(17.1168551874)+1.34239103051@-9*f(22 c
.1510903794)
s12 = s12+3.06160163503@-12*f(28.4879672510)+8.14807746743@-16*f(37 c
.0991210445)
c = mod(mod(s12)-mod(s9))
if (mod(s12))>1 then c = c/(mod(s12))
->nrl3 unless c<eps
antgrl = s12
test = c
return
nrl3:s15 = 2.18234885940@-1*f(.093307812017)+3.42210177923@-1*f(.492691740302)
s15 = s15+2.63027577942@-1*f(1.21559541207)+1.26425818106@-1*f(2 c
.26994952620)
s15 = s15+4.02068649210@-2*f(3.66762272175)+8.56387780361@-3*f(5 c
.42533662741)
s15 = s15+1.2124361472@-3*f(7.56591622661)+1.11674392344@-4*f(10.1202285680)
s15 = s15+6.45992676202@-6*f(13.1302824822)+2.22631690710@-7*f(16 c
.6544077083)
s15 = s15+4.22743038498@-9*f(20.7764788994)+3.92189726704@-11*f(25 c
.6238942267)
s15 = s15+1.45651526407@-13*f(31.4075191697)+1.48302705111@-16*f(38 c
.5306833065)
s15 = s15+1.60059490621@-20*f(48.0260855727)
c = mod(mod(s15)-mod(s12))
if (mod(s15))>1 then c = c/(mod(s15))
antgrl = s15
test = c
->nrl4 unless c<eps
return
nrl4:
printstring("
METHOD FAILED. LAST APPROXIMATION TO VALUE OF INTEGRAL IS ")
printfl(antgrl, 10)
printstring("
THE ERROR OF THIS ANSWER IS ")
printfl(test, 10)
return
end; !OF IMP DGALAG
!*END*IMPDGALAG
!*TOP*IMPDGALEG
external routine imp dgaleg(long real function f(long real x), long real xup, xlo,
long real name tin, eps, integer lim)
long real i5, i6, int5, int6, h, a, xu, xl
integer cntr, i
routine comp(long real a, long real name i5, i6,
long real function f(long real x))
long real c
c = 0.453089922969332*h
i5 = 0.118463442528095*(f(a+c)+f(a-c))
c = 0.269234655052842*h
i5 = i5+0.239314335249683*(f(a+c)+f(a-c))
i5 = h*(i5+0.284444444444445*f(a))
c = 0.466234757101576*h
i6 = 0.85662246189585@-1*(f(a+c)+f(a-c))
c = 0.330604693233133*h
i6 = i6+0.180380786524070*(f(a+c)+f(a-c))
c = 0.119309593041599*h
i6 = (i6+0.233956967286346*(f(a+c)+f(a-c)))*h
end; ! OF COMP
cntr = 0
a = .5*(xup+xlo)
h = xup-xlo
comp(a, i5, i6, f)
int5 = i5
int6 = i6
nrl1:
->nrl3 unless mod(int6-int5)<eps
tin = int6
eps = mod(int6-int5)
return
nrl3:cntr = cntr+1
->nrl2 if cntr<=lim
newlines(3)
printstring("
METHOD FAILED: REQUIRED ACCURACY NOT ATTAINED AFTER")
write(lim, 5)
printstring(" BISECTIONS OF THE INTERVAL")
newlines(3)
return
nrl2:xu = xup
int5 = 0
int6 = 0
h = h/2
cycle i = 1, 1, 2^^cntr
xl = xu-h
a = .5*(xu+xl)
comp(a, i5, i6, f)
int5 = int5+i5
int6 = int6+i6
xu = xl
repeat
->nrl1
end; ! OF ROUTINE IMP DGALEG
!*END*IMPDGALEG
!*TOP*IMPDIGVEC
external routine imp digvec(integer n, r, e, long real array (2) name a,
long real array name l, m, long real array (2) name x, y, integer array name it, vec, long real eps)
routine spec dirhesse(integer n, long real array (2) name a,
integer array name int)
routine spec qr(integer n, long real t, sum, prod,
integer array name its, long real array (2) name a, long real array name v, com)
routine spec reallu(long real array (2) name a,
integer array name perm, integer n, long real eps,
integer name d)
routine spec realsolve(long real array (2) name a, long real array name b, x,
integer array name perm, integer n)
routine spec comlu(long real array (2) name ax, ay,
integer array name perm, integer n, long real eps,
integer name d)
routine spec comsolve(long real array (2) name ax, ay, long real array name b, x,
integer array name perm, integer n)
long real s, t, p, q, sm, max, abs, norm, tol
integer i, k, d, rho, mult, sh, j, n2, ctr
long real array b, z(1:2*n), hx, hy(1:n, 1:n)
integer array inter, sort, perm(1:n)
if n<1 then ->er
if n=1 then start
l(1) = a(1, 1)
m(1) = 0
x(1, 1) = 1
y(1, 1) = 0
it(1) = 0
vec(1) = 0
return
finish
e = e//2
abs = 0
cycle i = 1, 1, n
cycle k = 1, 1, n
if abs<mod(a(i, k)) then abs = mod(a(i, k))
repeat
vec(i) = -1
it(i) = -n
sort(i) = i
repeat
if n=2 then start
s = (a(1, 1)+a(2, 2))/2
t = (a(1, 1)-a(2, 2))^2+4*a(1, 2)*a(2, 1)
p = sqrt(mod(t))/2
if t<0 then start
l(1) = s
l(2) = s
m(1) = p
m(2) = -p
finish else start
l(1) = s+p
l(2) = s-p
m(1) = 0
m(2) = 0
if l(1)<l(2) then start
t = l(1)
l(1) = l(2)
l(2) = t
finish
finish
it(1) = 0
it(2) = 0
finish else start
dirhesse(n, a, inter)
s = a(n, n)+a(n-1, n-1)
t = a(n, n)*a(n-1, n-1)-a(n, n-1)*a(n-1, n)
qr(n, eps, s, t, it, a, b, l)
cycle k = 1, 1, n
m(k) = b(k)^2+l(k)^2
repeat
cycle i = 1, 1, n-1
s = m(i)
cycle k = i+1, 1, n
if s<m(k) then start
m(i) = m(k)
m(k) = s
s = m(i)
d = sort(i)
sort(i) = sort(k)
sort(k) = d
finish
repeat
repeat
cycle k = 1, 1, n
m(k) = l(sort(k))
repeat
cycle k = 1, 1, n
l(k) = b(sort(k))
repeat
finish
if r=0 then return
if r>0 then d = 1 else d = -1
r = imod(r)
if r>n then r = n
ctr = 0
n2 = 2*n
tol = (n*eps*100)^(-2)
cycle rho = 1, 1, r
sh = 1
if d<0 then j = n+1-rho else j = rho
if it(j)<0 then ->v3
p = l(j)
q = m(j)
if ctr=1 then start
vec(j) = 0
cycle k = 1, 1, n
x(k, j) = x(k, j-d)
y(k, j) = -y(k, j-d)
repeat
ctr = 0
->v3
finish
if mod(q)<=eps*abs then ctr = 0 else ctr = 1
v1: cycle i = 1, 1, n
cycle k = 1, 1, n
hx(i, k) = a(i, k)
repeat
hx(i, i) = hx(i, i)-p
z(i+n) = 0
repeat
sm = 1/2^intpt(e/sh)
if ctr=1 then start
cycle i = 1, 1, n
cycle k = 1, 1, n
hy(i, k) = 0
repeat
hy(i, i) = -q
repeat
comlu(hx, hy, perm, n, sm, sh)
finish else reallu(hx, perm, n, sm, sh)
if sh<0 then start
sh = imod(sh)+1
->v1
finish
cycle k = 1, 1, n2
b(k) = 0
repeat
b(sort(j)) = 1
cycle i = 1, 1, 4
if ctr=1 then start
comsolve(hx, hy, b, z, perm, n)
finish else start
realsolve(hx, b, z, perm, n)
finish
max = 0
norm = 0
cycle k = 1, 1, n
p = z(k)^2+z(k+n)^2
norm = norm+p
if max<p then start
max = p
mult = k
finish
repeat
if norm>=tol then ->v2
p = z(mult)
q = z(mult+n)
cycle k = 1, 1, n
b(k) = (z(k)*p+z(k+n)*q)/max
b(k+n) = (z(k+n)*p-z(k)*q)/max
repeat
if norm<1@4 then start
mult = -1
cycle sh = 1, 1, n
k = sh
v5:
if mod(b(k))>9.999999999999989@-9 then start
b(k) = mult*b(k)
if k=sh then mult = -mult
finish else b(k) = b(k)+1000*eps
if k=sh and ctr=1 then start
k = sh+n
->v5
finish
repeat
finish
repeat
i = -4
v2: vec(j) = i
p = sqrt(norm)
cycle k = 1, 1, n2
b(k) = z(k)/p
repeat
z(1) = b(1)
z(2) = b(2)
z(n+1) = b(n+1)
z(n+2) = b(n+2)
if n=2 then ->v4
cycle sh = 3, 1, n
i = 0
v6: s = b(sh+i)
cycle k = 2, 1, sh-1
s = s+a(sh, k-1)*b(k+i)
repeat
z(sh+i) = s
if i=0 and ctr=1 then start
i = n
->v6
finish
repeat
cycle k = n-1, -1, 1
mult = inter(k)
if mult#k+1 then start
i = k+1
v7: s = z(mult)
z(mult) = z(i)
z(i) = s
if mult<=n then start
mult = mult+n
i = i+n
->v7
finish
finish
repeat
v4: norm = 0
cycle k = 1, 1, n
norm = norm+z(k)^2+z(k+n)^2
repeat
p = sqrt(norm)
cycle k = 1, 1, n
x(k, j) = z(k)/p
y(k, j) = z(k+n)/p
repeat
v3:repeat
return
er:newline
printstring("N TOO SMALL, NAMELY N = ")
write(n, 3)
printstring(" THIS ERROR OCCURRED IN ROUTINE == IMP DIGVEC ==")
newline
routine reallu(long real array (2) name a, integer array name perm,
integer n, long real eps, integer name d)
long real c
integer i, j, k, co
co = 0
cycle i = 1, 1, n-1
j = i+1
if mod(a(j, i))>mod(a(i, i)) then start
perm(i) = 1
cycle k = i, 1, n
c = a(i, k)
a(i, k) = a(j, k)
a(j, k) = c
repeat
finish else perm(i) = 0
if mod(a(i, i))<=eps then start
co = co+1
if co>d then start
d = -d
return
finish
a(i, i) = eps
finish
c = a(j, i)/a(i, i)
cycle k = j, 1, n
a(j, k) = a(j, k)-c*a(i, k)
repeat
a(j, i) = c
repeat
if mod(a(n, n))<=eps then start
co = co+1
if co>d then start
d = -d
return
finish
a(n, n) = eps
finish
end
routine realsolve(long real array (2) name a, long real array name b, x,
integer array name perm, integer n)
long real s
integer i, j, k
cycle i = 1, 1, n
x(i) = b(i)
repeat
cycle i = 1, 1, n-1
j = i+1
if perm(i)#0 then start
s = x(j)
x(j) = x(i)
x(i) = s
finish
x(j) = x(j)-a(j, i)*x(i)
repeat
x(n) = x(n)/a(n, n)
cycle i = n-1, -1, 1
s = x(i)
cycle k = i+1, 1, n
s = s-a(i, k)*x(k)
repeat
x(i) = s/a(i, i)
repeat
end
routine comlu(long real array (2) name ax, ay,
integer array name perm, integer n, long real eps,
integer name d)
long real cx, cy, mod
integer i, j, k, co
co = 0
cycle i = 1, 1, n-1
j = i+1
cx = ax(j, i)^2+ay(j, i)^2
cy = ax(i, i)^2+ay(i, i)^2
if cx>cy then start
perm(i) = 1
cycle k = i, 1, n
cx = ax(i, k)
ax(i, k) = ax(j, k)
ax(j, k) = cx
cx = ay(i, k)
ay(i, k) = ay(j, k)
ay(j, k) = cx
repeat
finish else perm(i) = 0
mod = ax(i, i)^2+ay(i, i)^2
if mod<=eps^2 then start
co = co+1
if co>d then start
d = -d
return
finish
ax(i, i) = eps
ay(i, i) = 0
mod = eps^2
finish
cx = (ax(j, i)*ax(i, i)+ay(j, i)*ay(i, i))/mod
cy = (ay(j, i)*ax(i, i)-ax(j, i)*ay(i, i))/mod
cycle k = j, 1, n
mod = ay(j, k)-cx*ay(i, k)-cy*ax(i, k)
ax(j, k) = ax(j, k)-cx*ax(i, k)+cy*ay(i, k)
ay(j, k) = mod
repeat
ax(j, i) = cx
ay(j, i) = cy
repeat
mod = ax(n, n)^2+ay(n, n)^2
if mod<=eps^2 then start
co = co+1
if co>d then start
d = -d
return
finish
ax(n, n) = eps
ay(n, n) = 0
finish
end
routine comsolve(long real array (2) name ax, ay, long real array name b, x,
integer array name perm, integer n)
long real sx, sy, mod
integer i, j, k
cycle i = 1, 1, 2*n
x(i) = b(i)
repeat
cycle i = 1, 1, n-1
j = i+1
if perm(i)#0 then start
mod = x(j)
x(j) = x(i)
x(i) = mod
mod = x(j+n)
x(j+n) = x(i+n)
x(i+n) = mod
finish
x(j) = x(j)-ax(j, i)*x(i)+ay(j, i)*x(i+n)
x(j+n) = x(j+n)-ax(j, i)*x(i+n)-ay(j, i)*x(i)
repeat
mod = ax(n, n)^2+ay(n, n)^2
sx = (x(n)*ax(n, n)+x(2*n)*ay(n, n))/mod
x(2*n) = (x(2*n)*ax(n, n)-x(n)*ay(n, n))/mod
x(n) = sx
cycle i = n-1, -1, 1
sx = x(i)
sy = x(i+n)
cycle k = i+1, 1, n
sx = sx-ax(i, k)*x(k)+ay(i, k)*x(k+n)
sy = sy-ax(i, k)*x(k+n)-ay(i, k)*x(k)
repeat
mod = ax(i, i)^2+ay(i, i)^2
x(i) = (sx*ax(i, i)+sy*ay(i, i))/mod
x(i+n) = (sy*ax(i, i)-sx*ay(i, i))/mod
repeat
end
routine dirhesse(integer n, long real array (2) name a,
integer array name int)
long real m
integer i, j, k, l
if n<2 then ->er
cycle i = 1, 1, n
int(i) = i+1
repeat
cycle k = 1, 1, n
if k=1 then ->l1
cycle i = 1, 1, n
m = a(i, k)
l = k-1
if k=n then ->l2
cycle j = k+1, 1, n
m = m+a(i, j)*a(j, l)
repeat
l2: if i<=2 then ->l4
if i>k then ->l3
cycle j = i-1, -1, 2
m = m-a(i, j-1)*a(j, k)
repeat
->l4
l3: cycle j = k, -1, 2
m = m-a(i, j-1)*a(j, k)
repeat
l4: a(i, k) = m
repeat
l1: l = k+1
i = l
if i>=n then ->l5
m = mod(a(i, k))
cycle j = i+1, 1, n
if mod(a(j, k))<=m then ->l8
m = mod(a(j, k))
i = j
l8: repeat
int(k) = i
if i=l then ->l7
cycle j = 1, 1, n
m = a(l, j)
a(l, j) = a(i, j)
a(i, j) = m
repeat
cycle j = 1, 1, n
m = a(j, l)
a(j, l) = a(j, i)
a(j, i) = m
repeat
l7: m = a(l, k)
if m=0 or k>n-2 then ->l5
cycle i = k+2, 1, n
a(i, k) = a(i, k)/m
repeat
l5: repeat
return
er: newline
printstring("N TOO SMALL, NAMELY N = ")
print(n, 3, 0)
printstring(" THIS ERROR OCCURRED IN ROUTINE == DIRHESSE ==")
newline
end
routine qr(integer n, long real t, sum, prod,
integer array name its, long real array (2) name a, long real array name v, com)
long real z, norm, s, x, p, q
integer i, j, k, l, m, u, preset
long real array aa(1:n, 1:n), c(1:6)
switch branch(1:3)
if n<3 then ->er
preset = 1
if sum*prod#0 then ->q1
preset = -1
q1: norm = 0
cycle i = 1, 1, n
cycle k = 1, 1, n
aa(i, k) = a(i, k)
if i<k+2 then ->q2
aa(i, k) = 0
q2:
norm = norm+mod(aa(i, k))
repeat
its(i) = 0
repeat
z = norm*t
l = 1
u = n
l2: cycle i = u, -1, l+1
if mod(aa(i, i-1))<z then ->l1
repeat
if its(u)>=40 then start
sum = 0
cycle i = 1, 1, u
sum = sum+mod(aa(l, i))
repeat
cycle i = l+1, 1, u
prod = 0
cycle j = i-1, 1, u
prod = prod+mod(aa(i, j))
repeat
if sum>prod then ->q3
sum = prod
q3: repeat
prod = sum*sum
sum = 1.5*sum
preset = 1
if its(u)>100*n then start
cycle i = 1, 1, u
its(i) = -u
v(i) = 0
com(i) = 0
repeat
return
finish
finish
if preset=1 then ->l9
sum = aa(u-1, u-1)+aa(u, u)
prod = aa(u-1, u-1)*aa(u, u)-aa(u-1, u)*aa(u, u-1)
l9: if u-2<l+1 then ->m3
cycle i = u-2, -1, l+1
c(1) = aa(i, i)*(aa(i, i)-sum)+prod+aa(i, i+1)*aa(i+1, i)
if c(1)#0 then start
c(2) = aa(i+1, i)*(aa(i, i)+aa(i+1, i+1)-sum)
c(3) = 0
if i+2>u then ->q4
c(3) = aa(i+2, i+1)*aa(i+1, i)
q4:
if mod(aa(i, i-1)*(mod(c(2))+mod(c(3)))/c(1))<z then start
m = i
->l8
finish
finish
repeat
m3: m = l
c(1) = aa(l, l)*(aa(l, l)-sum)+prod+aa(l, l+1)*aa(l+1, l)
c(2) = aa(l+1, l)*(aa(l, l)+aa(l+1, l+1)-sum)
c(3) = aa(l+2, l+1)*aa(l+1, l)
l8: if m>u-1 then ->m4
cycle k = m, 1, u-1
if k#m then start
c(1) = aa(k, k-1)
c(2) = aa(k+1, k-1)
c(3) = 0
if k+2>u then ->q5
c(3) = aa(k+2, k-1)
->q5
finish
its(u) = its(u)+1
q5: s = c(1)*c(1)+c(2)*c(2)+c(3)*c(3)
s = sqrt(s)
if c(1)>0 then ->q6
s = -s
q6:
if mod(s)<z then start
if k=m then start
p = aa(u-1, u-1)-aa(u, u)
s = sum
q = aa(u-1, u)*aa(u, u-1)
if p*p+4*q>0 and preset=-1 then start
q = sqrt(p*p+4*q)
x = (s+q)/2
if mod(s+q)<=mod(s-q) then ->q7
x = (s-q)/2
finish
x = s/2
q7: s = x
cycle i = m, 1, u
aa(i, i) = aa(i, i)-s
repeat
cycle k = m, 1, u-1
x = sqrt(aa(k, k)*aa(k, k)+aa(k+1, k)*aa(k+1, k))
v(k) = aa(k, k)/x
com(k) = aa(k+1, k)/x
aa(k, k) = x
aa(k+1, k) = 0
cycle j = k+1, 1, u
x = v(k)*aa(k, j)+com(k)*aa(k+1, j)
aa(k+1, j) = v(k)*aa(k+1, j)-com(k)*aa(k, j)
aa(k, j) = x
repeat
repeat
cycle k = m, 1, u-1
cycle i = k+1, -1, l
x = v(k)*aa(i, k)+com(k)*aa(i, k+1)
aa(i, k+1) = v(k)*aa(i, k+1)-com(k)*aa(i, k)
aa(i, k) = x
repeat
repeat
cycle i = m, 1, u
aa(i, i) = aa(i, i)+s
repeat
preset = -1
->l2
finish
aa(k+1, k-1) = 0
if k+2>u then ->l7
aa(k+2, k-1) = 0
->l7
finish
c(1) = c(1)+s
c(4) = c(1)/s
c(5) = c(2)/s
c(6) = c(3)/s
c(2) = c(2)/c(1)
c(3) = c(3)/c(1)
if k#m then start
aa(k, k-1) = -s
aa(k+1, k-1) = 0
if k+2>u then ->q8
aa(k+2, k-1) = 0
->q8
finish
if m=l then ->q8
aa(k, k-1) = -aa(k, k-1)
q8: cycle j = k, 1, u
x = aa(k, j)+c(2)*aa(k+1, j)
if k+2>u then ->q9
x = x+c(3)*aa(k+2, j)
aa(k+2, j) = aa(k+2, j)-x*c(6)
q9: aa(k, j) = aa(k, j)-x*c(4)
aa(k+1, j) = aa(k+1, j)-x*c(5)
repeat
j = u
if k+3>u then ->m1
j = k+3
m1: cycle i = 1, 1, j
x = c(4)*aa(i, k)+c(5)*aa(i, k+1)
if k+2>u then ->m2
x = x+c(6)*aa(i, k+2)
aa(i, k+2) = aa(i, k+2)-x*c(3)
m2: aa(i, k) = aa(i, k)-x
aa(i, k+1) = aa(i, k+1)-x*c(2)
repeat
l7: repeat
m4: preset = -1
->l2
l1: if i=u then start
branch(2):v(u) = aa(u, u)
com(u) = 0
u = u-1
->l3
finish
if i=u-1 then start
branch(3):p = aa(i, i)-aa(u, u)
x = aa(i, i)+aa(u, u)
q = aa(u, i)*aa(i, u)
s = aa(i, i)*aa(u, u)
if p*p+4*q<=0 then start
com(u) = (sqrt(-p*p-4*q))/2
com(u-1) = -com(u)
v(u) = x/2
v(u-1) = v(u)
finish else start
j = 1
if x>=0 then ->q0
j = -1
q0: v(u) = (x+j*sqrt(p*p+4*q))/2
v(u-1) = (s-q)/v(u)
com(u) = 0
com(u-1) = 0
finish
u = u-2
->l3
finish
l = i
->l9
l3: if u-l>=2 then ->l2
i = l
->branch(u-l+2)
branch(1):if l#1 then start
l = 1
if u>2 then ->l2
i = 2
if u=2 and mod(aa(2, 1))<z then ->l1
i = 1
->l1
finish
return
er: newline
printstring("N TOO SMALL, NAMELY N = ")
print(n, 3, 0)
printstring(" THIS ERROR OCCURRED IN ROUTINE == QR ==")
newline
end
end
!*END*IMPDIGVEC
!*TOP*IMPDKUSCR
external routine imp dkuscr(routine aux(long real array name f, y, long real x), long real x0, x1, e,
long real array name y, integer n, k)
routine spec int step(long real array name y, long real x, h,
integer n, long real name e, routine aux(long real array name f, y, long real x))
integer m, i, j, p, q, r
long real h, e1, e2, e3, x
long real array z(1:n)
m = 1
nrl2:
if m>=k then ->nrl1
m = 2*m
->nrl2
nrl1:h = (x1-x0)/m
p = 1; q = m; e2 = 0; e3 = 0; x = x0
nrl8:r = 1
nrl10:cycle i = p, 1, q
cycle j = 1, 1, n
z(j) = y(j)
repeat
int step(y, x, h, n, e1, aux)
if e1>e then ->nrl3
if e1<1.000000000000000@-2*e and r=1 and i#q then ->nrl4
nrl9: e2 = 0; e3 = 0
x = x+h
repeat
return
nrl3:
if e2=0 and e3=0 then ->nrl5
if e3=0 then ->nrl6
if e1>e2 and e1>e3 then ->nrl7
nrl6:e3 = e2
nrl5:e2 = e1
p = 1;
q = 2*q-2*i+2
h = h/2
cycle j = 1, 1, n
y(j) = z(j)
repeat
->nrl8
nrl7:
e = 2*e1
->nrl9
nrl4:
if fracpt((q-i)/2)>0 then ->nrl9
p = 1; q = int((q-i)/2); x = x+h
h = 2*h
->nrl8
routine int step(long real array name y, long real x, h,
integer n, long real name e, routine aux(long real array name f, y, long real x))
long real array a, b, c, d, q, r(1:n)
long real u, m; integer k, i
k = n
aux(y, a, x)
cycle i = 1, 1, k
a(i) = h*a(i)
q(i) = -5.555555559999999@-2*a(i)
r(i) = 7.916666666999999@-1*a(i)
y(i) = z(i)+2.222222221999999@-1*a(i)
repeat
u = x+2.222222221999999@-1*h
aux(y, b, u)
cycle i = 1, 1, k
b(i) = h*b(i)
r(i) = r(i)-3.375000000000000@0*b(i)
y(i) = z(i)+8.333333333999993@-2*a(i)+2.500000000000000@-1*b(i)
repeat
u = x+3.333333332999999@-1*h
aux(y, c, u)
cycle i = 1, 1, k
c(i) = h*c(i)
q(i) = q(i)+1.588235293999999@-1*c(i)
r(i) = r(i)+2.850000000000000@0*c(i)
y(i) = z(i)+5.390625000000000@-1*a(i)-1.898437499999999@0*b(i)+2 c
.109374999999999@0*c(i)
repeat
u = x+7.500000000000000@-1*h
aux(y, d, u)
cycle i = 1, 1, k
d(i) = h*d(i)
q(i) = q(i)-2.666666666999999@-1*d(i)
r(i) = r(i)-2.666666666999999@-1*d(i)
y(i) = z(i)-3.104999999999999@-1*a(i)+1.822499999999999@0*b(i)-1 c
.101600000000000@0*c(i)+4.895999999999999@-1*d(i)
repeat
u = x+9.000000000000000@-1*h
aux(y, b, u)
e = 0
cycle i = 1, 1, k
b(i) = h*b(i)
q(i) = q(i)+1.633986927999999@-1*b(i)
y(i) = z(i)+1.049382715999999@-1*a(i)+4.764705881999999@-1*c(i)+2 c
.370370370000000@-1*d(i)+1.815541030999999@-1*b(i)
if (d(i)-a(i))=0 then ->nrl1
m = -(q(i)*(i))/(d(i)-a(i))
if mod(m)<=e then ->nrl1
e = mod(m)
nrl1: repeat
end
end
!*END*IMPDKUSCR
!*TOP*IMPDSCRCH
external routine imp dscrch(long real function p(long real x), long real function f(long real x), long real x,
long real name y, long real array name a, integer n,
integer name ind)
routine spec copiv(long real array (2) name a, b, integer n, m,
long real eps, long real name det, integer name sing)
long real pi, val, th, sum
integer k, r, s, d
long real array pv(0:2*n+1), fv(0:n+1), aux(0:2*n+1, 1:1), ch(0:2*n+1,
0:2*n+1)
routine spec expand(integer d, long real array name c,
long real function fct(long real x))
if ind#0 then ->nrl1
if n<0 then ->nrl9
if n=0 then ->nrl7
if n=1 then ->nrl8
pi = 4*x-2
sum = a(n)
val = a(n-1)+pi*sum
cycle k = n-2, -1, 0
th = sum
sum = val
val = a(k)-th+pi*sum
repeat
y = (val-th)/2
return
nrl7:y = a(0)/2
return
nrl8:y = a(0)/2+(2*x-1)*a(1)
return
nrl9:ind = 1
->nrl10
nrl1:
if n<2 then ->nrl10
pi = 3.14159265359/2
d = 2*n+1
cycle r = 0, 1, d
cycle k = 0, 1, d
ch(r, k) = 0
repeat
repeat
ch(0, 0) = 1
ch(1, 1) = 2
ch(0, 1) = -1
cycle r = 2, 1, d
cycle k = 0, 1, r
ch(k, r) = -ch(k, r-2)-2*ch(k, r-1)
if k=0 then ->nrl2
ch(k, r) = ch(k, r)+4*ch(k-1, r-1)
nrl2: repeat
repeat
expand(d, pv, p)
expand(n+1, fv, f)
a(-1) = (mod(pv(d)))
a(-2) = (mod(fv(n+1)))
cycle k = 1, 1, n
val = ch(k, k)
cycle r = k-1, -1, 0
val = val*x+ch(r, k)
repeat
a(k) = val
repeat
cycle r = 1, 1, n
sum = pv(r-1)-pv(r+1)
aux(r, 1) = 2*(fv(r-1)-fv(r+1)-sum*y)
cycle s = 1, 1, n
if r=s then k = 1
if r#s then k = 0
pi = 8*k*r-2*sum*a(s)
th = pv(r+s-1)-pv(r+s+1)
if r=s then ->nrl6
k = (imod(r-s))
val = pv(k-1)-pv(k+1)
if r<s then val = -val
th = th+val
nrl6: ch(r, s) = pi+th
repeat
repeat
copiv(ch, aux, n, 1, 1@-10, th, ind)
if ind=1 then return
sum = y
cycle k = 1, 1, n
sum = sum-a(k)*aux(k, 1)
a(k) = aux(k, 1)
repeat
a(0) = 2*sum
return
nrl10:
printstring("
N TOO SMALL, NAMELY N = ")
print(n, 3, 0)
printstring(" THIS ERROR OCCURRED IN ROUTINE ")
printstring("== IMP DSCRSH ==
")
routine expand(integer d, long real array name c,
long real function fct(long real x))
long real array z(0:2*n+1)
cycle k = 0, 1, d
th = (2*k+1)*pi/(d+1)
th = cos(th/2)^2
z(k) = th
aux(k, 1) = fct(th)
repeat
cycle k = 0, 1, d
sum = 0
cycle r = 0, 1, d
val = ch(k, k)
if k=0 then ->nrl104
cycle s = k-1, -1, 0
val = val*z(r)+ch(s, k)
repeat
nrl104: sum = sum+val*aux(r, 1)
repeat
c(k) = 2*sum/(d+1)
repeat
s = -1
cycle k = 0, 1, d
if s>=0 then ->nrl102
if (mod(c(k)))<=4.999999999999999@-5 then s = 0
->nrl101
nrl102: s = k
if (mod(c(k)))<=(mod(c(k-1))) then ->nrl101
if (mod(c(k)))<=4.999999999999999@-5 then ->nrl103
s = -1
nrl101:repeat
s = -1
nrl103:if s<=0 then return
cycle k = s, 1, d
c(k) = 0
repeat
end
routine copiv(long real array (2) name a, b, integer n, m,
long real eps, long real name det, integer name sing)
long real c, sum
integer i, j, k, ri, ci, con
!========================================================================
! FILE CORRUPTION! Guessed missing lines here:
integerarray s(1:n)
cycle i = 1, 1, n
s(i) = i
repeat
cycle i = 1, 1, n
! Complete pivoting: search for maximum element in the submatrix
!========================================================================
cycle j = i+1, 1, n
c = a(j, i)/a(i, i)
cycle k = i+1, 1, n
a(j, k) = a(j, k)-c*a(i, k)
repeat
cycle k = 1, 1, m
b(j, k) = b(j, k)-c*b(i, k)
repeat
a(j, i) = c
repeat
repeat
nrl4:cycle k = 1, 1, m
b(n, k) = b(n, k)/a(n, n)
repeat
cycle i = n-1, -1, 1
cycle j = 1, 1, m
sum = 0
cycle k = i+1, 1, n
sum = sum+a(i, k)*b(k, j)
repeat
b(i, j) = (b(i, j)-sum)/a(i, i)
repeat
repeat
c = 1
cycle i = 1, 1, n
c = c*a(i, i)
cycle k = 1, 1, n
if s(k)=i then ->nrl5
repeat
nrl5: cycle j = 1, 1, m
sum = b(k, j)
b(k, j) = b(i, j)
b(i, j) = sum
repeat
con = s(k)
s(k) = s(i)
s(i) = con
repeat
sing = 0
det = (mod(c))
return
nrl10:
printstring("
N OR M TOO SMALL, NAMELY N = ")
print(n, 3, 0)
printstring(" M = ")
print(m, 3, 0)
printstring("
THIS ERROR OCCURRED IN ROUTINE ")
printstring("== COPIV ==
")
end
end
!*END*IMPDSCRCH
!*TOP*IMPFLEPO
external routine imp flepo(routine funct(real name f, real array name x, g), integer n, limit, iprin,
integer name conv, real eps, real name f, real array name x, h)
real oldf, sg, ghg, z, step, ita, fa, fb, ga, gb, w, lambda
integer i, j, k, count, fnevalns
real array g, s, gamma, sigma, y(1:n)
real function dot(real array name a, b)
integer i
real s
s = 0
cycle i = 1, 1, n
s = s+a(i)*b(i)
repeat
result = s
end
real function updot(integer i, real array name a, b)
comment MULTIPLYBYTHEITHROWOFTHESYMMETRICMATRIXA, c
!WHOSE UPPER TRIANGLE IS STORED BY ROWS.
integer j, k; real s
k = i; s = 0
->nrl1 if i=1
cycle j = 1, 1, i-1
s = s+a(k)*b(j)
k = k+n-j
repeat; !STEPS TO DIAGONAL:NOW GO ALONG ROW.
nrl1: cycle j = i, 1, n
s = s+a(k+j-i)*b(j)
repeat
result = s
end; !OF FN UPDOT
nrl1:
->nrl2 unless conv=1
k = 1
cycle i = 1, 1, n
h(k) = 1;
->nrl30 if i=n
cycle j = 1, 1, n-i
h(k+j) = 0
repeat
nrl30:k = k+n-i+1
repeat; !FORMATION OF UNIT MATRIX IN H.
nrl2:conv = 1; step = 1; fnevalns = 0
funct(f, x, g); fnevalns = fnevalns+1
oldf = f+1
cycle count = 1, 1, limit
return unless oldf>f
cycle i = 1, 1, n
sigma(i) = x(i); gamma(i) = g(i); s(i) = -updot(i, h, g)
repeat
fb = f; gb = dot(g, s)
->nrl21 unless gb>=0
conv = -1; return
nrl21:oldf = f; ita = step
! ITA IS DETERMINED FROM ITS VALUE AT THE PREVIOUS ITERATION (STEP)
! AND IS INCREASED OR DECREASED BY 4 WHERE NECESSARY. IT SHOULD
! TEND TO 1 AT THE MINIMUM.
nrl4: fa = fb; ga = gb
cycle i = 1, 1, n; x(i) = x(i)+ita*s(i); repeat
funct(f, x, g); fb = f; gb = dot(g, s); fnevalns = fnevalns+1
->nrl3 unless (gb<0 and fb<fa)
ita = 4*ita; step = 4*step;
->nrl4
nrl3: if ita<=0 then start; conv = -1; return; finish
cycle i = 1, 1, n; y(i) = g(i); repeat
z = 3*(fa-fb)/ita+ga+gb
w = sqrt(z*z-ga*gb)
lambda = ita*(gb+w-z)/(gb-ga+2*w)
cycle i = 1, 1, n; x(i) = x(i)-lambda*s(i); repeat
funct(f, x, g); fnevalns = fnevalns+1
->nrl6 unless (f>fa or f>fb)
step = step/4
->nrl5 unless fb<fa
cycle i = 1, 1, n; x(i) = x(i)+lambda*s(i); g(i) = y(i)
repeat
f = fb;
->nrl6
nrl5: gb = dot(g, s)
->nrl7 unless (gb<0 and count>n and step<9.999999999999971@-7)
conv = -1; return
nrl7: fb = f; ita = ita-lambda;
->nrl3
nrl6: cycle i = 1, 1, n
sigma(i) = x(i)-sigma(i)
gamma(i) = g(i)-gamma(i)
repeat
sg = dot(sigma, gamma)
->nrl13 unless count>=n
return if sqrt(dot(s, s))<eps and sqrt(dot(sigma, sigma))<eps
comment end OFTESTFORVANISHINGDERIVATIVE.
nrl13:cycle i = 1, 1, n
s(i) = updot(i, h, gamma)
repeat
ghg = dot(s, gamma)
k = 1
->nrl8 if mod(sg)<9.999999999999987@-31 or c
mod(ghg)<9.999999999999987@-31
cycle i = 1, 1, n
cycle j = i, 1, n
h(k) = h(k)+sigma(i)*sigma(j)/sg-s(i)*s(j)/ghg
k = k+1; repeat
repeat
->nrl20 if iprin=0
->nrl20 unless count//iprin*iprin=count
nrl8: newlines(2);
printstring("
ITERATION ")
print(count, 2, 0)
printstring(":
NO. OF FN EVALNS = ")
write(fnevalns, 2)
newlines(2);
printstring(" X G
")
cycle i = 1, 1, n
printfl(x(i), 6); spaces(4); printfl(g(i), 6); newline
repeat
printstring("
PRESENT MINIMUM VALUE OF FN =")
printfl(f, 6)
newlines(2)
nrl20:repeat; ! COUNT LOOP; CONV=-1 ;
printstring("
PROBLEM DID NOT CONVERGE IN SPECIFIED NO OF ITERATIONS
ROUTINE TERMINATED.")
newlines(2)
end; !OF ROUTINE IMP FLEPO
!
!*END*IMPFLEPO
!*TOP*IMPGAHER
external routine imp gaher(real function f(real x), real eps,
real name antgrl, test)
long real i4, i6, i9, i12, i15, c
i4 = .813128354472@-1*(f(1.65068012388)+f(-1.65068012388))+.804914090005* c
(f(.524647623275)+f(-.524647623275))
i6 = .453000990551@-2*(f(2.35060497367)+f(-2.35060497367))+.157067320323* c
(f(1.33584907401)+f(-1.33584907401))+.724629595224*(f(.436077411928)+f c
(-.436077411928))
c = mod(mod(i6)-mod(i4))
if (mod(i6))>1 then c = c/(mod(i6))
->nrl1 unless c<eps
antgrl = i6
test = c
return
nrl1:i9 = .396069772632@-4*(f(3.19099320178)+f(-3.19099320178))+ c
.494362427553-2*(f(2.26658058453)+f(-2.26658058453))
i9 = i9+.884745273943-1*(f(1.46855328921)+f(-1.46855328921))+ c
.432651559002*(f(.723551018752)+f(-.723551018752))+.720235215606*f(0)
c = mod(mod(i9)-mod(i6))
if (mod(i9))>1 then c = c/(mod(i9))
->nrl2 unless c<eps
antgrl = i9
test = c
return
nrl2:i12 = .265855168436@-6*(f(3.88972489787)+f(-3.88972489787))+ c
.857368704358@-4*(f(3.02063702512)+f(-3.02063702512))+.390539058463@-2* c
(f(2.27950708050)+f(-2.27950708050))
i12 = i12+.516079856158@-1*(f(1.59768263515)+f(-1.59768263515))+ c
.260492310264*(f(.947788391240)+f(-.947788391240))+.570135236262*(f( c
.314240376254)+f(-.314240376254))
c = mod(mod(i12)-mod(i9))
if (mod(i12))>1 then c = c/(mod(i12))
->nrl3 unless c<eps
antgrl = i12
test = c
return
nrl3:i15 = .152247580425@-8*(f(4.49999070731)+f(-4.49999070731))+ c
.105911554771@-5*(f(3.66995037340)+f(-3.66995037340))
i15 = i15+.100004441232@-3*(f(2.96716692791)+f(-2.96716692791))+ c
.277806884291@-2*(f(2.32573248617)+f(-2.32573248617))
i15 = i15+.307800338725@-1*(f(1.71999257519)+f(-1.71999257519))+ c
.158488915796*(f(1.13611558521)+f(-1.13611558521))
i15 = i15+.412028687499*(f(.565069583256)+f(-.565069583256))+ c
.564100308726*f(0)
c = mod(mod(i15)-mod(i12))
if (mod(i15))>1 then c = c/(mod(i15))
antgrl = i15
test = c
return if c<eps
printstring("
METHOD FAILED TO SATISFY GIVEN ACCURACY. LAST")
printstring(" ESTIMATE OF INTEGRAL IS ")
printfl(antgrl, 10)
printstring("
THE ERROR OF THIS ANSWER IS ")
printfl(test, 10)
return
end; !OF IMP GAHER
!*END*IMPGAHER
!*TOP*IMPGALAG
external routine imp galag(real function f(real x), real eps,
real name antgrl, test)
long real s4, s6, s9, s12, s15, c
s4 = 6.03154104342@-1*f(0.322547689619)+3.57418692438@-1*f(1.74576110116) c
+3.88879085150@-2*f(4.53662029692)+5.39294705561@-4*f(9.395070912301)
s6 = 4.58964673950@-1*f(.222846604179)+4.17000830772@-1*f(1.18893210167)
s6 = s6+1.13373382074@-1*f(2.99273632606)+1.03991974531@-2*f(5 c
.77514356910)+2.61017202815@-4*f(9.83746741838)+8.98547906430@-7*f(15 c
.9828739806)
c = mod(mod(s6)-mod(s4))
if (mod(s6))>1 then c = c/(mod(s6))
->nrl1 unless c<eps
antgrl = s6
test = c
return
nrl1:s9 = 3.36126421798@-1*f(.152322227732)+4.11213980424@-1*f(.807220022742)
s9 = s9+1.99287525371@-1*f(2.00513515562)+4.74605627657@-2*f(3.78347397333)
s9 = s9+5.59962661079@-3*f(6.20495677788)+3.05249767093@-4*f(9.37298525169)
s9 = s9+6.59212302608@-6*f(13.4662369111)+4.11076933035@-8*f(18.8335977890)
s9 = s9+3.29087403035@-11*f(26.3740718909)
c = mod(mod(s9)-mod(s6))
if (mod(s9))>1 then c = c/(mod(s9))
->nrl2 unless c<eps
antgrl = s9
test = c
return
nrl2:s12 = 2.64731371055@-1*f(.115722117358)+3.77759275873@-1*f(.61175748452) c
+9.04492222117@-2*f(2.83375133774)
s12 = s12+2.44082011320@-1*f(1.51261026977)
s12 = s12+2.01023811546@-2*f(4.59922763942)+2.66397354187@-3*f(6 c
.84452545312)
s12 = s12+2.03231592663@-4*f(9.62131684246)+8.36505585682@-6*f(13 c
.0060549933)
s12 = s12+1.66849387654@-7*f(17.1168551874)+1.34239103051@-9*f(22 c
.1510903794)
s12 = s12+3.06160163503@-12*f(28.4879672510)+8.14807746743@-16*f(37 c
.0991210445)
c = mod(mod(s12)-mod(s9))
if (mod(s12))>1 then c = c/(mod(s12))
->nrl3 unless c<eps
antgrl = s12
test = c
return
nrl3:s15 = 2.18234885940@-1*f(.093307812017)+3.42210177923@-1*f(.492691740302)
s15 = s15+2.63027577942@-1*f(1.21559541207)+1.26425818106@-1*f(2 c
.26994952620)
s15 = s15+4.02068649210@-2*f(3.66762272175)+8.56387780361@-3*f(5 c
.42533662741)
s15 = s15+1.2124361472@-3*f(7.56591622661)+1.11674392344@-4*f(10.1202285680)
s15 = s15+6.45992676202@-6*f(13.1302824822)+2.22631690710@-7*f(16 c
.6544077083)
s15 = s15+4.22743038498@-9*f(20.7764788994)+3.92189726704@-11*f(25 c
.6238942267)
s15 = s15+1.45651526407@-13*f(31.4075191697)+1.48302705111@-16*f(38 c
.5306833065)
s15 = s15+1.60059490621@-20*f(48.0260855727)
c = mod(mod(s15)-mod(s12))
if (mod(s15))>1 then c = c/(mod(s15))
antgrl = s15
test = c
->nrl4 unless c<eps
return
nrl4:
printstring("
METHOD FAILED. LAST APPROXIMATION TO VALUE OF INTEGRAL IS ")
printfl(antgrl, 10)
printstring("
THE ERROR OF THIS ANSWER IS ")
printfl(test, 10)
return
end; !OF IMP GALAG
!*END*IMPGALAG
!*TOP*IMPGALEG
external routine imp galeg(real function f(long real x), real xup, xlo,
real name tin, eps, integer lim)
long real i5, i6, int5, int6, h, a, xu, xl
integer cntr, i
routine comp(long real a, long real name i5, i6, real function f(long real x))
long real c
c = .45308992*h
i5 = .11846344*(f(a+c)+f(a-c))
c = .26923466*h
i5 = i5+0.23931434*(f(a+c)+f(a-c))
i5 = h*(i5+.28444444*f(a))
c = .46623476*h
i6 = .85662246@-1*(f(a+c)+f(a-c))
c = .33060469*h
i6 = i6+.18038079*(f(a+c)+f(a-c))
c = .11930959*h
i6 = (i6+0.23395697*(f(a+c)+f(a-c)))*h
end; !END OF COMP
cntr = 0
a = .5*(xup+xlo)
h = xup-xlo
comp(a, i5, i6, f)
int5 = i5
int6 = i6
nrl1:
->nrl3 unless mod(int6-int5)<eps
tin = int6
eps = mod(int6-int5)
return
nrl3:cntr = cntr+1
->nrl2 if cntr<=lim
newlines(3)
printstring("
METHOD FAILED: REQUIRED ACCURACY NOT ATTAINED AFTER")
write(lim, 5)
printstring(" BISECTIONS OF THE INTERVAL")
newlines(3)
return
nrl2:xu = xup
int5 = 0
int6 = 0
h = h/2
cycle i = 1, 1, 2^^cntr
xl = xu-h
a = .5*(xu+xl)
comp(a, i5, i6, f)
int5 = int5+i5
int6 = int6+i6
xu = xl
repeat
->nrl1
end; ! OF ROUTINE IMP GALEG
!*END*IMPGALEG
!*TOP*IMPLOGFCT
external real function imp logfct(integer i, j)
integer n, k
real x, p
switch s(0:20)
if j=1 then p = 0.4342944819 else p = 1.0
n = i//50; if n>20 then n = 20
if i<0 then start
printstring(" ARGUMENT I IS NEGATIVE. ")
newline
monitor
stop
finish
->s(n)
s(0):x = 0;
->nrl1
s(1):x = 148.4777669516;
->nrl1
s(2):x = 363.7393755532;
->nrl1
s(3):x = 605.0201058475;
->nrl1
s(4):x = 863.2319871918;
->nrl1
s(5):x = 1134.0452317854;
->nrl1
s(6):x = 1414.058499460;
->nrl1
s(7):x = 1704.1247472880;
->nrl1
s(8):x = 2000.5006979990;
->nrl1
s(9):x = 2303.1351597674;
->nrl1
s(10):x = 2611.3304584723;
->nrl1
s(11):x = 2924.5291022234;
->nrl1
s(12):x = 3242.2753354039;
->nrl1
s(13):x = 3564.1895888475;
->nrl1
s(14):x = 3889.9508323156;
->nrl1
s(15):x = 4219.2839911876;
->nrl1
s(16):x = 4551.9507307617;
->nrl1
s(17):x = 4887.7425519021;
->nrl1
s(18):x = 5226.4755155920;
->nrl1
s(19):x = 5567.9861427925;
->nrl1
s(20):x = 5912.1281786647
nrl1:if i=50*n then result = x*p
cycle k = (50*n)+1, 1, i
x = x+log(k)
repeat
result = x*p
end
!*TOP*IMPDAPRO
external routine imp dapro(long real function f(long real x, long real array name v), integer k0, ma, me, n,
integer name m, long real name s, long real array name x, y, v,
tol, sf)
routine spec cholesky(integer n, long real array (2) name a, long real array name b,
integer name fault)
integer i, k, l, loop, ind
long real function spec af(integer j, long real x)
long real array dy(ma:me), b(1:n), c(1:n, 1:n)
long real z, z1
if n<2 then ->nrl100
loop = 0
nrl1:loop = loop+1
cycle k = 1, 1, n
b(k) = 0
cycle l = 1, 1, n
c(k, l) = 0
repeat
repeat
cycle i = ma, 1, me
dy(i) = y(i)-f(x(i), v)
if k0=0 then ->nrl101
dy(i) = dy(i)-v(k0)
nrl101:repeat
cycle i = ma, 1, me
cycle k = 1, 1, n
z = af(k, x(i))
b(k) = b(k)+dy(i)*z
cycle l = k, 1, n
z1 = z^2
if l=k then ->nrl102
z1 = af(l, x(i))*z
nrl102: c(l, k) = c(l, k)+z1
repeat
repeat
repeat
cholesky(n, c, b, ind)
if ind=1 then ->nrl2
cycle k = 1, 1, n
v(k) = v(k)+b(k)
repeat
cycle k = 1, 1, n
if (mod(b(k)))>tol(k) then ->nrl103
repeat
m = loop
->nrl3
nrl103:
if loop<m then ->nrl1
nrl2:m = -loop-1
nrl3:s = 0
cycle i = ma, 1, me
sf(i) = f(x(i), v)
if k0=0 then ->nrl105
sf(i) = sf(i)+v(k0)
nrl105:s = s+(y(i)-sf(i))^2
repeat
s = sqrt(s/(me-ma+1))
return
nrl100:newline; spaces(6)
printstring("THIS ERROR OCCURRED IN ROUTINE ===IMP DAPRO===")
newline
long real function af(integer j, long real x)
long real zw1, zw2, h
if j=k0 then result = 1
h = tol(j); if h<5@-2 then h = 5@-2
if 0.1<h then h = 0.1
v(j) = v(j)+h/20
zw1 = f(x, v)
v(j) = v(j)-h/10
zw2 = f(x, v)
v(j) = v(j)+h/20
result = 10*(zw1-zw2)/h
end
routine cholesky(integer n, long real array (2) name a, long real array name b,
integer name fault)
long real c, sum
integer i, j, k
fault = 1
if n<2 then ->nrl10
cycle i = 1, 1, n-1
cycle j = i+1, 1, n
a(i, j) = 0
repeat
repeat
if a(1, 1)<=0 then return
c = sqrt(a(1, 1))
a(1, 1) = c
cycle k = 2, 1, n
a(1, k) = a(k, 1)/c
repeat
cycle i = 2, 1, n
sum = 0
cycle j = 1, 1, i-1
sum = sum+a(j, i)*a(j, i)
repeat
if a(i, i)<=sum then return
a(i, i) = sqrt(a(i, i)-sum)
if i=n then ->nrl1
cycle k = i+1, 1, n
sum = 0
cycle j = 1, 1, i-1
sum = sum+a(j, i)*a(j, k)
repeat
a(i, k) = (a(k, i)-sum)/a(i, i)
repeat
nrl1: repeat
b(1) = b(1)/a(1, 1)
cycle i = 2, 1, n
sum = 0
cycle k = 1, 1, i-1
sum = sum+a(k, i)*b(k)
repeat
b(i) = (b(i)-sum)/a(i, i)
repeat
b(n) = b(n)/a(n, n)
cycle i = n-1, -1, 1
sum = 0
cycle k = i+1, 1, n
sum = sum+a(i, k)*b(k)
repeat
b(i) = (b(i)-sum)/a(i, i)
repeat
fault = 0
return
nrl10:
printstring("
N TOO SMALL, NAMELY N = ")
print(n, 3, 0)
printstring(" THIS ERROR OCCURRED IN ROUTINE ")
printstring( c
"== CHOLESKY == CH
")
end
end
!*END*IMPDAPRO
!*TOP*IMPDCHLIN
external routine imp dchlin(integer m, n, integer name ifault,
long real array (2) name a, long real array name b, x, long real name eps)
long real max, rho, sum, rhold
integer first, u, v, w, k, np, i
long real array apt(1:n, 1:n+1), llt(1:n, 1:n), s(1:n+1), bt(1:m)
integer array j(1:n+2), ref(1:n+1)
routine spec solve(integer tr, long real array name in, out)
routine spec mult(integer tr, r, c, long real array (2) name a, long real array name in, out)
routine spec replace(integer in, out)
integer function spec sign(long real x)
ifault = 1
if n<2 or m<=n then ->er
first = 0
np = n+1
cycle u = 1, 1, np
ref(u) = u
repeat
c1:cycle v = 1, 1, np
w = ref(v)
cycle u = 1, 1, n
apt(u, v) = a(w, u)
repeat
bt(v) = b(w)
repeat
cycle u = 1, 1, n
cycle v = 1, 1, u
sum = 0
cycle w = 1, 1, np
sum = sum+apt(u, w)*apt(v, w)
repeat
llt(u, v) = sum
repeat
repeat
if llt(1, 1)<=0 then ->c2
rho = sqrt(llt(1, 1))
llt(1, 1) = rho
cycle w = 2, 1, n
llt(1, w) = llt(w, 1)/rho
repeat
cycle u = 2, 1, n
rho = llt(u, u)
cycle v = 1, 1, u-1
rho = rho-llt(v, u)*llt(v, u)
repeat
if rho<=0 then ->c2
rho = sqrt(rho)
llt(u, u) = rho
if u=n then ->c3
cycle w = u+1, 1, n
sum = 0
cycle v = 1, 1, u-1
sum = sum+llt(v, u)*llt(v, w)
repeat
llt(u, w) = (llt(w, u)-sum)/rho
repeat
repeat
c3:solve(1, bt, x)
mult(-1, np, n, apt, x, s)
sum = 0
max = 0
cycle u = 1, 1, np
rho = bt(u)-s(u)
s(u) = rho
sum = sum+rho*rho
max = max+mod(rho)
repeat
if max<=eps/np then rho = 0
if max>eps/np then rho = sum/max
if first=0 then ->c4
if rho>rhold then ->c4
c2:if first=-1 then return
if first=0 then ->c5
i = i+1
if j(i)=0 then ->c5
replace(j(i-1), j(i))
->c1
c5:first = -1
cycle u = 1, 1, np
ref(u) = m-n+u-1
repeat
->c1
c4:rhold = rho
if rho<eps then ->c6
cycle u = 1, 1, np
bt(u) = bt(u)-rho*sign(s(u))
repeat
solve(1, bt, x)
c6:mult(1, m, n, a, x, bt)
max = 0
k = 0
cycle u = 1, 1, m
sum = b(u)-bt(u)
if max>mod(sum) then ->c7
max = mod(sum)
k = u
w = sign(sum)
c7:repeat
if mod(max-rho)<eps then ->c8
cycle u = 1, 1, n
x(u) = w*a(k, u)
repeat
solve(0, x, x)
mult(-1, np, n, apt, x, bt)
max = 0
cycle u = 1, 1, np
j(u) = 0
if mod(s(u))<eps*mod(bt(u)) then ->c9
sum = bt(u)/s(u)
if max>=sum then ->c9
max = sum
j(1) = u
c9:repeat
j(np+1) = 0
if j(1)=0 then ->l1
i = 0
cycle u = 1, 1, np
if mod(s(u))<eps*mod(bt(u)) then ->l2
sum = bt(u)/s(u)
if mod(max-sum)>eps then ->l2
i = i+1
j(i) = ref(u)
l2:repeat
->l3
l1:j(1) = ref(1)
l3:replace(k, j(1))
i = 1
first = 1
->c1
c8:ifault = 0
eps = rho
return
er:newline
printstring("M OR N TOO SMALL, NAMELY M = ")
print(m, 3, 0)
printstring(" N = ")
print(n, 3, 0)
newline
printstring("THIS ERROR OCCURRED IN ROUTINE ")
printstring("== IMP DCHLIN ==")
newline
routine solve(integer tr, long real array name in, out)
long real s
integer i, k
if tr=1 then mult(1, n, np, apt, in, out)
out(1) = out(1)/llt(1, 1)
cycle i = 2, 1, n
s = 0
cycle k = 1, 1, i-1
s = s+llt(k, i)*out(k)
repeat
out(i) = (out(i)-s)/llt(i, i)
repeat
out(n) = out(n)/llt(n, n)
cycle i = n-1, -1, 1
s = 0
cycle k = i+1, 1, n
s = s+llt(i, k)*out(k)
repeat
out(i) = (out(i)-s)/llt(i, i)
repeat
end
routine mult(integer tr, r, c, long real array (2) name a, long real array name in, out)
long real s, m
integer i, k
cycle i = 1, 1, r
s = 0
cycle k = 1, 1, c
if tr=1 then m = a(i, k)
if tr#1 then m = a(k, i)
s = s+m*in(k)
repeat
out(i) = s
repeat
end
routine replace(integer in, out)
integer k
cycle k = 1, 1, np
if ref(k)=out then ->u1
repeat
u1: ref(k) = in
end
integer function sign(long real x)
if x>0 then result = 1
if x<0 then result = -1
result = 0
end
end
!*END*IMPDCHLIN
!*TOP*IMPDELI3
external long real function imp deli3(long real x, a, b)
long real am, ap, c, ca, cb, d, de, e, f, fa, g, h, hh, ln2, p1, pi, pm,
pz, q, r, s, t, u, v, w, y, ye, z, zd
integer k, l, m, n, bo, bk
long real array ra, rb, rr(2:14)
long real function spec arsinh(long real x)
integer function spec sign(long real x)
if x=0 then result = 0
ln2 = 0.6931471805599
pi = 3.141592653590
ca = 1@-6; cb = 1@-14
hh = x*x; f = b*hh
if a#0 then s = a else s = ca/(1+mod(x))
t = s*s; pm = t/2
e = hh*t;
z = mod(f)
r = mod(b)
h = 1+hh
if e<0.1 and z<0.1 and t<1 and r<1 then ->sl
w = 1+f
if w=0 then ->er
if b#0 then p1 = b else p1 = cb/hh
s = mod(s)
y = mod(x)
g = p1-1; if g=0 then g = cb
f = p1-t; if f=0 then f = cb*t
am = 1-t; ap = 1+e
r = p1*h; fa = g/(f*p1)
if fa>0 then bo = 1 else bo = -1
fa = mod(fa)
pz = mod(g*f)
de = sqrt(pz)
q = sqrt(mod(p1))
if pm>0.5 then pm = 0.5
pm = p1-pm
if pm>=0 then start
v = y*de; u = sqrt(r*ap)
if g<0 then v = -v
d = 1/q; c = 1
finish else start
u = sqrt(h*ap*pz)
ye = y*q; v = am*ye
q = -de/g; d = -am/de
c = 0; pz = ap-r
finish
if bo=1 then start
r = v/u
z = 1; k = 1
if pm<0 then start
h = y*sqrt(h/(ap*fa))
h = 1/h-h
z = h-2*r; r = 2+r*h
if r=0 then r = cb
if z=0 then z = h*cb
w = r/z; z = w
r = w; w = pz
finish
u = u/w; v = v/w
finish else start
t = u+mod(v)
bk = 1
if p1<0 then start
de = v/pz; ye = 2*u*ye
u = t/pz; v = (-f-g*e)/t
t = pz*mod(w)
z = (hh*r*f-g*ap+ye)/t
ye = ye/t
finish else start
de = v/w; ye = 0
u = (e+p1)/t
v = t/w; z = 1
finish
if s>1 then start
h = u; u = v
v = h
finish
finish
y = 1/y; e = s
n = 1; t = 1
l = 0; m = 0
mr:y = y-e/y
if y=0 then y = cb*sqrt(e)
f = c; c = c+d/q
g = e/q; d = 2*(f*g+d)
q = q+g; g = t
t = s+t; n = 2*n
m = 2*m
if bo=1 then start
if z<0 then m = m+k
k = sign(r); h = e/(u*u+v*v)
u = u*(1+h); v = v*(1-h)
finish else start
r = u/v; h = z*r
z = h*z; hh = e/v
if bk=1 then start
de = de/u
ye = ye*(h+1/h)+de*(1+r)
de = de*(u-hh)
if mod(ye)<1 then bk = 1 else bk = -1
finish else start
z = log(z)
k = intpt(z/ln2)+1
z = exp(z-k*ln2)
m = m+k
finish
finish
if mod(g-s)>ca*g then start
if bo=1 then start
g = (1/r-r)/2
hh = u+v*g
h = g*u-v
if hh=0 then hh = u*cb
if h=0 then h = v*cb
z = r*h
r = hh/h
finish else start
u = u+e/u
v = v+hh
finish
s = 2*sqrt(e)
e = s*t
l = 2*l
if y<0 then l = l+1
->mr
finish
if y<0 then l = l+1
e = arctan(1, t/y)+pi*l
e = e*(c*t+d)/(t*(t+q))
if bo=1 then start
h = v/(t+u)
z = 1-r*h
h = r+h
if z=0 then z = cb
if z<0 then m = m+sign(h)
s = arctan(1, h/z)+m*pi
finish else start
if bk=1 then s = arsinh(ye) else s = log(z)+m*ln2
s = s/2
finish
e = (e+sqrt(fa)*s)/n
result = e*sign(x)
sl:cycle k = 2, 1, 10
rb(k) = 0.5/k
ra(k) = 1-rb(k)
repeat
zd = 0.0454545454545
s = b+pm
cycle k = 2, 1, 10
rr(k) = s
pm = pm*t*ra(k)
s = s*b+pm
repeat
u = s*zd
s = u
bo = -1
cycle k = 10, -1, 2
u = u+(rr(k)-u)*rb(k)
bo = -bo
v = -u*bo
s = s*hh+v
repeat
if bo=1 then s = -s
u = (u+1)/2
result = (u-s*h)*sqrt(h)*x+u*arsinh(x)
er:newline
printstring("1+B*X*X = 0, NAMELY B = ")
printfl(b, 10)
printstring(" X = ")
printfl(x, 10)
newline
printstring("THIS ERROR OCCURRED IN ROUTINE == IMP DELI3 ==")
newline
result = 1@17
integer function sign(long real x)
if x>0 then result = 1
if x=0 then result = 0
result = -1
end
long real function arsinh(long real x)
long real a, y, b0, b1, b2
integer i
long real array c(0:13)
a = x*x
if a<1 then start
c(0) = 1.8717996008826; c(1) = -0.0588117611900
c(2) = 0.0047274654322; c(3) = -0.0004938363163
c(4) = 0.0000585062071; c(5) = -0.0000074669983
c(6) = 0.0000010011694; c(7) = -0.0000001390354
c(8) = 0.0000000198232; c(9) = -0.0000000028847
c(10) = 0.0000000004267; c(11) = -0.0000000000640
c(12) = 0.0000000000097; c(13) = -0.0000000000015
b1 = 0.0000000000002
y = 2*a-1
b0 = 2*y*b1+c(13)
cycle i = 12, -1, 0
b2 = b1
b1 = b0
b0 = c(i)-b2+2*y*b1
repeat
result = x*(b0-b2)/2
finish
y = mod(x)
if a>1@16 then b0 = log(y)+0.6931471805599 else b0 = log(y+sqrt(1+a))
if x<0 then b0 = -b0
result = b0
end
end
!*END*IMPDELI3
end of file