\\ Solution of Problem 5

Digits(dig)=if(dig>0, default(realprecision,dig)); precision(0.0); 
cpoints(x)=[exp(I*x[1]),exp(I*x[2]),-1,exp(-I*x[2]),exp(-I*x[1])]~;
igamma(x)=if(imag(x)==0&real(x)<=0&real(x)==floor(x),0,1/gamma(x));  
p5error(z,c)=local(e,i); e=z; for(i=1,length(z), e[i]=igamma(z[i])-c[1]+z[i]*(c[2]+z[i]*(c[3]+z[i]*c[4]))); e;
residual(A,b)=b-A*matsolve(conj(A~)*A,conj(A~)*b);
{p5solve(z)=local(A,b,k,j,d); A=matrix(5,4,k,j,z[k]^(j-1)); b=vector(5,k,igamma(z[k]))~;
  d=precision(residual(A,b),Digits); for(k=1,5,d[k]=d[k]/abs(d[k]));
  matsolve(concat(A,d),b)}
optfun(x)=real(p5solve(cpoints(x))[5]);
{gauss_newton(x)=local(h,H,F,J,k,j,xk,xj);
  h=10^(-precision(x[1])/3);  H=matrix(2,2); F=matrix(3,3); J=vector(2)~; 
  for ( k=-1,1, xk=x[1]; x[1]=x[1]+k*h; 
    for (j=-1,1, xj=x[2]; x[2]=x[2]+j*h;
      F[k+2,j+2]=optfun(x); x[2]=xj);
    x[1]=xk);
  J[1]=(F[3,2]-F[1,2])/(2*h); J[2]=(F[2,3]-F[2,1])/(2*h);  
  H[1,1]=(F[3,2]-F[2,2]+(F[1,2]-F[2,2]))/h^2; 
  H[2,2]=(F[2,3]-F[2,2]+(F[2,1]-F[2,2]))/h^2;
  H[1,2]=(F[3,3]-F[1,3]+(F[1,1]-F[3,1]))/(4*h^2);
  H[2,1]=H[1,2];
  matsolve(H,J)}

\\ x is output from Matlab program
allocatemem(320 000 000);
x=[1.403199183178580,2.262377449725701]~; dig=10; gettime;
{ for (k=1,10, dig=2*dig; ndig=floor(dig); Digits(ndig);  x=precision(x,ndig); 
  d=gauss_newton(x); x=x-precision(d,ndig); print("milliseconds=",gettime); print(optfun(x)) ) }

quit