{inprod(x,y) = sum(j=1,N1,real(x[j])*real(y[j])-imag(x[j])*imag(y[j]))}
{dig(eps) = if(eps==0,dec,-log(abs(eps))/log(10))}

{problem3(T,h,tol) =
\\ T: truncation point, h: step size, tol: error tolerance for power method
  N=floor(T/h); N1=N+1; s=vector(N1,k,h*(k-N1)); y=vector(N1);
  t=vector(N1,k,s[k]+s[k]^3/3); dt=vector(N1,k,h*(1+s[k]^2)); \\t=Phi(s)
  z=vector(N1,k,sqrt(1/2)*(cosh(t[k])-I*sinh(t[k]))); zbar=conj(z); \\contour
  z1=vector(N1,k, z[k]-0.5); z2=vector(N1,k, -2*z[k]+1.75); \\auxiliary vectors
  ct=vector(N1,k,if(abs(imag(z[k]))<925412,cotan(Pi*z[k]),sign(imag(z[k]))/I));
  w=vector(N1,k, -ct[k]*zbar[k]*dt[k]);  w[N1]=w[N1]/2;
  x0=vector(N1,j, 2/(z[j]^2+z[j]))~; x=x0/x0[N1];
  lam=0; lam0=1; iter=0;
  while ( abs(lam-lam0)>tol, lam0=lam;
    xw=vector(N1,j, x[j]*w[j]);
    for(k=1,N1, z1k=z1[k]; z2k=z2[k]; y[k]=sum(j=1,N1,xw[j]/((z[j]+z1k)^2+z2k))\
           + sum(j=1,N1,conj(xw[j]) / ((zbar[j]+z1k)^2 + z2k)) );
    yw=vector(N1,k, y[k]*w[k]); x0=x;
 for(j=1,N1,z1j=z1[j]-1;z2j=1.5-z2[j];x[j]=sum(k=1,N1,yw[k]/((z1j+z[k])^2+z2j))\
           + sum(k=1,N1, conj(yw[k]) /((z1j+zbar[k])^2 + z2j)) );
    lam=inprod(x,xw)/inprod(x0,xw); iter=iter+1;
\\  print(iter,"   ",lam);
  ); lam}

\\ Calling program

dec=48; default(realprecision,dec); 
L=log(10)*dec; tol=10^(2-dec);
T=L^(1/3); h=1/90; \\ h=1/(47.7*sinh(5/3*asinh(dec/51.4)));
res=problem3(T,h,tol); nrm=sqrt(res);
ex=1.62364719163297668812114803229217389051362004031757196115031496969474216802941770791914905330604705140840507431554829943395152091008273773620518163996048601062890070033956127396424095317410913181535526701701751709122805621388030494278219109072757261860509721348966967334902161;
d=dig(res-ex);