Andrew John Walker on Fri, 27 Feb 2004 01:16:49 +0100


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Puzzling error


I've lately been getting the error
  ***   division by zero in gdiv, gdivgs or ginv

in a script I run and would like to try and figure out what is causing
it. This is in the recent windows compile of version 2.2.8

The script uses the definitions:
\p18
ep=1.0E-6

pz(x)=sum(n=1,25000,moebius(n)*log(zeta(n*x))/n)-1
qz(x)=sum(n=1,25000,moebius(n)*log(zeta(n*x))/n)+1
pzz(x,a,b)=sum(n=a,b,moebius(n)*log(zeta(n*x))/n)
qzz(x,a,b)=sum(n=a,b,moebius(n)*log(zeta(n*x))/n)

and one of the functions it uses is the following, which is
a higher order version of Newton's method which I have found works
very well (I wish I could find the name of it again!).

ha(x,zz) =
{
	local(a,b,c,da,db,delta,et,diff);
	a=zz;print1(abs(a));
                  if(real(x)>0.2,ep=1.0E-5,ep=1.0E-05);
	if(abs(a)>1.0,print("");return(100+100.0*I));
	if(abs(a)>0.1,
		if(abs(a)>0.50 && n>=6,return(100+100.0*I));
		if(abs(a)>0.15 && n>=10,return(100+100.0*I));
		print1(" * ");
		b=pzz(x+ep,1,200);print1(" *");
		c=pzz(x-ep,1,200);print("*");
		da=(b-c)/(2*ep);print("da=",da);
		db=(b+c-2*a)/(ep*ep);print("db=",db);
		delta=a/da;print("delta=",delta);
		et=a*db/(da*da);print("et=",et);
		diff=delta/sqrt(1-et);print("diff=",diff);
		return(x-diff);
	);
	a=a+pzz(x,201,25000);
	if(abs(a)>1.0,print("");return(100+100.0*I));
	print1(" * ");print1(abs(a));
	b=pzz(x+ep,1,25000);print1(" *");
	c=pzz(x-ep,1,25000);print("*");
	da=(b-c)/(2*ep);
	db=(b+c-2*a)/(ep*ep);
	delta=a/da;
	et=a*db/(da*da);
	diff=delta/sqrt(1-et);
	return(x-diff);
}

Running the following produces the error:

(10:54) gp > n=1
%2 = 1
(10:54) gp > x=0.42+1.852*I
%3 = 0.420000000000000000 + 1.85200000000000000*I
(10:54) gp > q=pzz(x,1,200)
%4 = -0.371525983207137192 - 0.824222420202395294*I
(10:55) gp > x=ha(x,q)
0.904087470415514532 *  **
da=0.232984177950373828 + 0.396420328478054736*I
db=0.909890443 + 0.1362236290*I
delta=-1.95476905912811306 - 0.211654834752745888*I
et=-2.788520878 + 2.775115125*I
diff=-0.826979100 - 0.3732423968*I
%5 = 1.246979100 + 2.225242396*I
(10:55) gp > q=pzz(x,1,200)
%6 = -0.2430868957 - 0.4396598426*I
(10:55) gp > x=ha(x,q)
0.502386322 *  **
da=0.2492102794 + 0.2378073987*I
db=-1.164153219 + 4.656612874*I
delta=-1.391679841 - 0.4362102564*I
et=-4.301330500 - 19.86198354*I
diff=-0.3020276940 + 0.1106739148*I
%7 = 1.549006794 + 2.114568482*I
(10:55) gp > q=pzz(x,1,200)
%8 = -0.1516423058 - 0.3915455483*I
(10:55) gp > x=ha(x,q)
0.4198848714 *  **
da=0.1827633242 + 0.2369284630*I
db=0.E1 + 0.E1*I
delta=-1.345613821 - 0.3979536607*I
et=0.E1 + 0.E1*I
  ***   division by zero in gdiv, gdivgs or ginv

Does the value of et and db shown indicate that not enough
precision was available? Using for instance \p27 stops
the error, but I'd prefer to stick with \p18 and somehow
trap or avoid this error due to speed issues. Most
starting points in the zero finding don't produce this
problem, it's only a very small fraction.

Thanks for any help in advance,
Andrew