Puzzling error

[Andrew John Walker <ajw01@uow.edu.au>]

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