Hi
I have a not too complicated code that calls a function in a loop. The
problem is that if the loop has 100 or so iterations then I get this
message:
Maxima encountered a Lisp error:
Error in PROGN [or a callee]: Bind stack overflow.
Even when the number of iterations is say 50ish if I run it two or
three times then again I get that same message. I can restart Maxima
and it will work again, but only with short enough loops.
In case you want to take a look at the code I'm attaching it here. I
know it's long and I could've perhaps come up with a shorter version
which generates the same error without irrelevant parts of the code
but I thought maybe there's an obvious problem that someone could
point out to me. Most of the definitions prior to and within the
indented code (where nrfordv is defined) are irrelevant. The problem
arises in the last line which calls the function nrford causing 200
calls of nrfordv (the long indented code). It fails at the 87th
iteration. However, if I change
nrford(.1,makelist(i,i,1,200),1);
to
nrford(1,makelist(i,i,1,200),1);
it works!
I appreciate any hint toward the resolution.
Thanks,
Mahdiyar
load(draw)$
load(basic)$
N:60$
D:0.01$
V:1$
assume(v>0)$
Z(x,d):=(1+d*x^2)^2/(1+d*x^2+6*d^2*x^2)$
y(x,d):=sqrt(1+6*d)/sqrt(d)*asinh(sqrt(d*(6*d+1))*x)-sqrt(6)*atanh((sqrt(6)*d*x)/sqrt(d*(6*d+1)*x^2+1))$
intres:-sqrt(6*d^2+d)/d*log((sqrt(6*d^2+d)*x+sqrt(1+(6*d^2+d)*x^2))/((sqrt(6*d^2+d)*x0+sqrt(1+(6*d^2+d)*x0^2))))+sqrt(6*d^2)/2/d*log((1+x*sqrt(6*d^2)/sqrt(1+(6*d^2+d)*x^2))*(1-x0*sqrt(6*d^2)/sqrt(1+(6*d^2+d)*x0^2))/(1-x*sqrt(6*d^2)/sqrt(1+(6*d^2+d)*x^2))/(1+x0*sqrt(6*d^2)/sqrt(1+(6*d^2+d)*x0^2)))$
yn(x,d,x0):=''intres$ /* should be used for d<-1/6 and 1+d*x^2<0 */
V(x,d,v):=(x^2-v^2)^2/(1+d*x^2)^2/4$
epx(x,d,v):=''((diff(V(x,d,v),x))^2/(V(x,d,v))^2/2)$
intres:factor(epx(x,d,v)*Z(x,d))$
epsilon(x,d,v):=''intres;
etx(x,d,v):=''((diff(V(x,d,v),x,2))/V(x,d,v))$
intres:factor(Z(x,d)*etx(x,d,v)+''(diff(Z(x,d),x))*(''(diff(V(x,d,v),x)))/V(x,d,v)/2)$
eta(x,d,v):=''intres;
intres:factor(expand(epsilon(sqrt(dx2/d),d,v)))$
epsilon2(dx2,d,v):=''intres$
intres:factor(expand(eta(sqrt(dx2/d),d,v)))$
eta2(dx2,d,v):=''intres$
Delta:sqrt(factor(V(x,d,v)/(24*%pi^2*Z(x,d)*epsilon(x,d,v))))$
/*assume(x2>x1,x1>0,v>x2)$
intres:abs(logcontract(-integrate(V(x,d,v)/''diff(V(x,d,v),x)/Z(x,d),
x, x1, x2)))$
Ne(x1,x2,d,v):= ''intres;*/
Ne(x1,x2,d,v):=abs(3/4*log((1+d*x2^2)/(1+d*x1^2))+(v^2*log(x2^2/x1^2)+(1+6*d)*(x1^2-x2^2))/(8*(d*v^2+1)));
intres:Ne(sqrt(dxs2/d),sqrt(dxe2/d),d,v)$
Ne2(dxs2,dxe2,d,v):=''intres;
ratprint: false$
intres:factor(subst(sqrt(dx2/d),x,diff(V(x,d,v),x)))$
DV2(dx2,d,v):=''intres$
nrfordv(d,v,print_details):=block([], /*dxs2,a,b,s,et,ep,xe,xl,xr,nrlist*/
if print_details>-1 then print("d=",d,"v=",v),
nrlist:[],
ep:makelist([0,0],i,1,2),
et:makelist([0,0],i,1,5),
infint:makelist([0,0],i,1,10),
s:ev(rectform(allroots(8*d*dx2*(d*v^2+1)^2=(6*d*dx2+dx2+1)*(d*v^2-dx2)^2,dx2)),numer),
a:[rhs(s[1]),rhs(s[2]),rhs(s[3])]*signum(d),
if print_details=2 then print("epsilon=1:",a),
for i:1 thru 3 do (
if (abs(imagpart(a[i]))>1d-5 or realpart(a[i])<0 or (d<0 and
d>-1/6 and realpart(a[i])>1/(1+6*d))) then
a[i]:0
else a[i]:realpart(a[i])
),
a:sort(a),
a:delete(0,a),
if length(a)=0 then print("error") else ep[1]:[0,a[1]],
if length(a)=2 then ep[2]:[a[2],inf],
if length(a)=3 then ep[2]:[a[2],a[3]],
ep:delete([0,0],ep),
if print_details=2 then print("ep:",ep),
s:append(ev(rectform(allroots(4*d*(d*v^2+1)*(24*d^2*dx2^2*v^2+4*d*dx2^2*v^2+3*d*dx2*v^2-d*v^2-12*d*dx2^3-2*dx2^3+12*d*dx2^2+dx2^2+3*dx2)=(6*d*dx2+dx2+1)^2*(d*v^2-dx2)^2,dx2)),numer),
ev(rectform(allroots(4*d*(d*v^2+1)*(24*d^2*dx2^2*v^2+4*d*dx2^2*v^2+3*d*dx2*v^2-d*v^2-12*d*dx2^3-2*dx2^3+12*d*dx2^2+dx2^2+3*dx2)=-(6*d*dx2+dx2+1)^2*(d*v^2-dx2)^2,dx2)),numer)),
b:[rhs(s[1]),rhs(s[2]),rhs(s[3]),rhs(s[4]),rhs(s[5]),rhs(s[6]),rhs(s[7]),rhs(s[8])]*signum(d),
if print_details=2 then print("|eta|=1:",b),
for i:1 thru 8 do (
if (abs(imagpart(b[i]))>1d-5 or realpart(b[i])<0 or (d<0 and
d>-1/6 and realpart(b[i])>last(a))) then
b[i]:0
else b[i]:realpart(b[i])
),
b:sort(b),
b:delete(0,b),
if abs(eta(0,d,v))<1 then (
if length(b)>0 then et[1]:[0,b[1]],
for i:2 step 2 thru 6 do (
if length(b)=i then et[floor(i/2)+1]:[b[i],inf],
if length(b)>i then et[floor(i/2)+1]:[b[i],b[i+1]] ),
if length(b)=8 then et[5]:[b[8],inf]
) else (
for i:1 step 2 thru 7 do (
if length(b)=i then et[floor((i+1)/2)]:[b[i],inf],
if length(b)>i then et[floor((i+1)/2)]:[b[i],b[i+1]] )
),
et:delete([0,0],et),
if print_details=2 then print("et:",et),
k:0,
for i:1 thru length(ep) do
for j:1 thru length(et) do (
k:k+1,
infint[k]:[max(ep[i][1],et[j][1]),min(ep[i][2],et[j][2])],
if infint[k][1]>=infint[k][2] then infint[k]:[0,0],
if (DV2(signum(d)*infint[k][1],d,v)<0 and
infint[k][2]=inf) then (infint[k]:[0,0],print("v'",i,j)),
if (infint[k][1]=0 and DV2(signum(d)*infint[k][2],d,v)>0)
then infint[k]:[0,0],
if (infint[k][1]#0 and infint[k][2]#inf and
Ne2(signum(d)*infint[k][1], signum(d)*infint[k][2], d, v)<2*N)
then infint[k]:[0,0]
),
infint:delete([0,0],infint),
if print_details>0 then print(infint),
for i:1 thru length(infint) do (
if infint[i][1]=0 then (dvl:0, xl:signum(d)*exp(-4*N)) else
/*v*exp(-4*N*(1+d*v^2)/v^2)*/
(dvl:DV2(signum(d)*infint[i][1],d,v), xl:signum(d)*infint[i][1]),
if infint[i][2]=inf then (dvr:0,
xr:infint[i][1]+10+10*8*abs(d)*N*(1+d*v^2)/(1+6*d)) else
(dvr:DV2(signum(d)*infint[i][2],d,v), xr:signum(d)*infint[i][2]),
if (dvl>=0 and dvr>=0) then xe:1 else
if (dvl<=0 and dvr<=0) then xe:2 else
(print("error"), return(1)),
dxs2:find_root(Ne2(xi,signum(d)*infint[i][xe],d,v)=N,xi,xl,xr),
if print_details=2 then print(dxs2),
push([1-6*epsilon2(dxs2,d,v)+2*eta2(dxs2,d,v),16*epsilon2(dxs2,d,v)],nrlist),
kill(dxs2)
),
if print_details>-1 then (
print(nrlist),
print("***")
),
return(nrlist)
)$
nrford(d,vl,pd):=block([nsrlist],
nsrlist:[],
for i:1 thru length(vl) do nsrlist:append(nrfordv(d,vl[i],pd),nsrlist),
return(nsrlist)
)$
nrford(.1,makelist(i,i,1,200),1);