\\ Meijer's G function as needed for Problem 9

{ prob9G(alfa) = 
  local(s,t,sterm,tprod,tsum,z,n); 
  s=0; t=0; n=0; z=alfa^2/16; 
  sterm=1; tprod=1; tsum=1-log(alfa/2)-psi(alfa+1)+2*psi(1);
  while (1,
    s=s+sterm; t=t+tsum*tprod; if (1.0+sterm==1, break);
    n=n+1; sterm=-sterm*z*(alfa-2*n)*(alfa-2*n+1)/(n*(n+1)*(2*n+1)^2);
    tprod=-tprod*z*(alfa-2*n+1)*(alfa-2*n+2)/(n^2*(2*n-1)*(2*n+1));
    tsum=tsum+1/(alfa-2*n+1)+1/(alfa-2*n+2)+1/n+1/(2*n-1)+1/(2*n+1);
   );
  2^(alfa)*alfa*(2*Pi*z*s+t)
}

\\ Numerical derivative of prob9G
{ dprob9G(alfa) =
  local(h,dh,d2h);
  h=10^(-precision(0.0)/5.0);
  dh=(prob9G(alfa+h)-prob9G(alfa-h))/(2*h);
  d2h=(prob9G(alfa+2*h)-prob9G(alfa-2*h))/(4*h);
  (4*dh-d2h)/3
}


{ timings() = 
default(realprecision,315);
gettime; s=psi(Pi/4); print(gettime);
s=psi(Pi/4+0.1); print(gettime);
s=prob9G(Pi/4); print(gettime);
s=dprob9G(Pi/4); print(gettime);
default(realprecision,630);
gettime; s=psi(Pi/4); print(gettime);
s=psi(Pi/4+0.1); print(gettime);
s=prob9G(Pi/4); print(gettime);
s=dprob9G(Pi/4); print(gettime);
default(realprecision,1260);
gettime; s=psi(Pi/4); print(gettime);
s=psi(Pi/4+0.1); print(gettime);
s=prob9G(Pi/4); print(gettime);
s=dprob9G(Pi/4); print(gettime);
default(realprecision,2520);
gettime; s=psi(Pi/4); print(gettime);
s=psi(Pi/4+0.1); print(gettime);
s=prob9G(Pi/4); print(gettime);
s=dprob9G(Pi/4); print(gettime);
\\ default(realprecision,5040);
\\ gettime; s=psi(Pi/4); print(gettime);
}

{ p9I(alfa) = (2+sin(10*alfa))*prob9G(alfa); }
{ dp9I(alfa) = (2+sin(10*alfa))*dprob9G(alfa)+10*cos(10*alfa)*prob9G(alfa); }

{ secant(x1, fx1, x2, fx2) = x2 - fx2*(x1-x2)/(fx1-fx2) }
global digits; digits=28;
{ Digits(n) = n=ceil(n); if (n>28, digits=n; default(realprecision,n));}
{ main(final) =
  default(realprecision,28);
  x1=Pi/4; fx1=dp9I(x1); print(fx1);
  x2=x1+fx1/100; fx2=dp9I(x2); print(fx2);
  while (1,
    x3=secant(x1,fx1,x2,fx2); acc3=floor(1.62*log(abs(x3-x2))/log(0.1));
    write(concat("prob9g.",acc3),x3);
    print("Error estimate = ",acc3);
    if (acc3>=final, break);
    if (2*acc3>4/5*digits, 
       Digits(min(2.5*acc3,1.25*final)); x3=precision(x3,digits));
    fx3=dp9I(x3); err=floor(2+log(abs(fx3))/log(0.1));
    print("Error bound = ",err,"  Heap state = ",getheap);
    Digits(1.62*err);
    x1=precision(x2,digits); fx1=precision(fx2,digits);
    x2=precision(x3,digits); fx2=precision(fx3,digits);
  );
  x3
}