Solutions to Trefethen's problems from 
"Hundred-dollar, Hundred-digit challenge",
SIAM News volume 35, number 1.

Keith Briggs 2002 Mar 25

Notes on high-precision arithmetic:

    http://keithbriggs.info/xrc.html
    www.nersc.gov/~dhbailey/mpdist/mpdist.html
    www.shoup.net/ntl/


Problem 1 
---------

  integrate_0^1 cos(log(x)/x)/x dx

  Solution:
  
    0.3233674316777787...

  Method:

     integrate_0^1 cos(log(x)/x)/x dx
   = integrate_0^infty W'(y) cos(y) dy
   = sum_{n=1; n odd}^infty integrate_{-pi}^{pi} W'(u+n*pi) cos(u) du ...(1)

  where W=Lambert's W function: W(x)exp(W(x))=x; W'(x)=W(x)/(x(1+W(x))).

  The advantage of this representation is that each integral is a product
  of cos with slowly varying function (W(x)~log(x)).   Thus, it is very well
  approximated by an expansion

  integrate_{-pi}^{pi} W'(u+n*pi) cos(u) du ~=
    a_2/2! W^{(3)}(n*pi) + a_4/4! W^{(5)}(n*pi) + ...                 ...(2)

  where a_k = integrate_{-pi}^{pi} x^k cos(x) dx.
  
  Both the a_k and the derivatives W^{(k+1)} of W can be computed from 
  simple recurrences.    Already at n=13, two terms of this expansion
  give sufficient accuracy.   Thus, we do the first few terms in the sum (1)
  by numerical integration, then switch to (2) for the rest of the sum.
  
  We may estimate the truncation error as follows: W^{(3)}(x) is asymptotically
  2/(x*log(x))^3 for large x.   Bounding this by 2/x^3, we can sum the tail 
  as being < 4*pi/x^2.   From this we estimate that 10^6 terms
  are sufficient for 10 dp in the final result.

  Doing a more accurate tail estimation:

     sum_{k=n}^infty W^{(3)}(k*pi) 
  ~= 2*pi sum_{k=n}^infty 1/(x*log(x))^3 
  ~  2*pi integral_n^infty 1/(x*log(x))^3 dx
  =  2*pi [ -2Ei(-2log(n)) + (1-2log(n))/(2n^2 log(n)^2) ]

  This suggests that the tail is:
  <8e-9   at n=1000 
  <4e-11  at n=10000
  <2e-13  at n=100000
  <2e-15  at n=1000000

Problem 2
---------

  Circular mirror of radius 1/3 at each point of Z^2.
  Photon starts at (1/2,1/10) with velocity (1,0).
  What is its position at time 10?

  Method:

  1. Run a quick-and-dirty python program (trefethen2.py) using
     IEEE double arithmetic to get a sequence of bounces, which we
     will afterwards confirm with precise arithmetic.
     This program will also output data for feeding to "graph" to make
     a postscript figure (trefethen2-fig.ps).
     The sequence of the centres of circular mirrors computed is
     [[1, 0], [-1, 1], [0,  2], [ 0,  1], [-1, 2], [ 0,  0], [-1, 0], [0, 1],
      [0, 0], [-1, 0], [0, -1], [-1, -1], [-1, 0], [-1, -1], [ 1, 0]];

  2. Run a C++ program (trefethen2-xr.cc) which will recompute
     everything in such a way that the final time and velocity vector
     are accurate to any requested precision (e.g. 10 decimals).
     If the sequence of bounces above was wrong, this program will crash
     with a "sqrt of negative argument" error.

  Solution:

     x(10)    = -0.736292698609618310775737171866
     y(10)    = -0.669642696363571375194376070966
     distance =  0.99526291944335416089031180942

  More precisely:
     x(10)   = -0.736292698609618310775737171866391720451575491031024118015570
                32135829096701522929400945112122000980773908774827834875197238
                77531369315078617517331539949084939398476219513951166249120103
                0045569289350602e0
     y(10)   = -0.669642696363571375194376070967098466862227372928775421597411
                09381282909191783522403449125396431985636370847144665160568810
                35768639394200390550524425747382398387223992432483782229188924
                2655423072136579e0 
     distance = 0.9952629194433541608903118094267216210294669227341543498032088
                580729861796228306320991749818976188760306062973261905117220872
                365706973813827798622298601793010235848291712247521533893738542
                1202854314982e0

Problem 3
---------

  A=
  1/1  1/2  1/4  1/7  1/11 ...
  1/3  1/5  1/8  1/12 1/17 ...
  1/6  1/9  1/13 1/18 1/24 ...
  1/10 1/14 1/19 1/25 1/32 ...
  1/15 1/20 1/26 1/33 1/41 ...
  ....

  What is the 2-norm of A considered as a operator on l^2?
  
  Solution:

    1.27422415274...

  Method:
  
    We have  1/A_{ij}=(i+1)(i+2(j+1))/2+j(j-1)/2 (indexing from 0)
    and therefore must maximize over x_j: sqrt(sum_i (sum_j A_{ij}*x_j)^2)
    subject to sum_i x_i^2=1
    
    I truncated these sums at n=4000 and performed the constrained
    optimization.


Problem 4
---------

  minimize  exp(sin(50x))
           +sin(60exp(y))
           +sin(70sin(x))
           +sin(sin(80y))
           -sin(10(x+y))
           +(x^2+y^2)/4

  Solution:

    -3.30686864747523728...
    at 
    x=-0.0244030796943785 y=2.10612427155358

  Method:
  
     Random search followed by gradient optimization.

     Check with gp:
     43 100
     x=-0.0244030796962585; y=0.210612427157718
     objective:
     exp(sin(50*x))+sin(60*exp(y))+sin(70*sin(x))+sin(sin(80*y))-sin(10*(x+y))+(x^2+y^2)/4
     gradient: 
     x/2+50*exp(sin(50*x))*cos(50*x)-10*cos(10*(x+y))+70*cos(x)*cos(70*sin(x)
     y/2+60*exp(y)*cos(60*exp(y))-10*cos(10*(x+y))+80*cos(80*y)*cos(sin(80*y))

Problem 5
---------

  Find the cubic minimizing sup-norm error for 1/gamma(z) on unit disk.

  Solution:
    
    min sup-norm error = 0.2143352346...

  Minimizing cubic:
    
    0.0055418508 + 1.0197619z + 0.62521197z^2 - 0.60334332z^3
                                                        
  Method:

  1. Since 1/gamma(z) has a Taylor series around 0 with real coefficients, 
  so will the best approximating polynomial.

  2. By the maximum modulus theorem, the maximum approximation error
  in the unit disk will occur on the boundary.

  3. A rough check shows the max error occurs near t=0.2233, 0.36 and at 1/2,
  where z=exp(2*pi*i*t).

  4. Expand 1/gamma(z) around the first two of these points to get 
  accurate values to fit to.  (1/gamma(-1)=0 for t=1/2)

  5. Use a gradient optimizer to get the 4 coefficients of the cubic; 
  to find the max error near the above two approximate t values use a 
  Brent one-dimensional optimizer in the argument t for each gradient step.

Problem 6 NB: this solution is incorrect!
---------

  Particle walks on Z^2.  
  P(up)=1/4, P(dn)=1/4, P(rt)=1/4+eps, P(lt)=1/4-eps.
  P(eventual return to origin)=1/2.
  What is eps?
  
  Solution:

    0.108996971771294...

  Method:

    Compute probabilities on 4000*4000 lattice, use bisection to find eps.

Problem 7
---------

  A is a 20000*20000 integer matrix
  
  A(i,i) = prime(i)
  A(i,j) = 1 if |i-j|=1,2,4,8,...,16384
         = 0 else
  
  What is (A^-1)_{11}?

  Solution: 
  
    0.72507834626840116746868771925116096886918059447950895787816...

  Method:

    Solve A*x = [1,0,0,0,...,0]^T; solution is x_1.
    Conjugate gradient works well since A is sparse symmetric 
    positive definite.   Neither A nor A^-1 need be stored.

  
Problem 8
---------
  
  Square domain [0,2]*[0,2], 3 edges held at 0, one edge (x=2) at 5.
  u_t=u_{xx}+u_{yy}
  At what time does u at (1,1) reach 1?

  Solution: 

    0.42401138703368836379743366859326...
  
  Method:

   Solving the heat equation with Fourier series, we find that the
   temperature u at (1,1) at time t is

     (5/4) + 40/pi^2 sum_{m=1; m odd}^{infty} sum_{n=1; n odd}^{infty} 
       (-1)^{(m+n)/2} n/m 1/(m^2+n^2) exp[-(m^2+n^2)*pi^2*t/4]

   where the 5/4 term is the equilibrium temperature.
   Bisection easily finds the required value of t.

Problem 9
---------

  argmax_{a in [0,5]} (2+sin(10a)) * integrate_0^2 x^a * sin(a/(2-x)) dx

  Solution: 
  
    a = 0.78593367197930...

  Method: 
  
    The integral is easier to handle if we rewrite it as

      integrate_{1/2}^infinity (2-1/u)^a/u^2 * sin(a*u) du

    I used an adaptive integration algorithm optimized for
    oscillatory integrands, together with a golden-search search.


Problem 10
----------

  Rectangle R=[-5,5]*[-1/2,1/2]
  Particle starts Brownian motion at (0,0)
  What is probability of first exit through a short side of R?

  Solution: 
 
    2*1.918793989625613051703566593102...e-7
    = 0.00000038375879792512261034071331862048...

  Method: 
  
    Unless I have misunderstood something, this problem is easy!
    Solve Fokker-Planck equation on R with appropriate boundary conditions,
    p(x,y)=prob(first exit through right edge of R starting at (x,y))
    which just gives Laplace's equation p_{xx} + p_{yy} = 0 
    BC: 0 on left, bottom, top edges; 1 on right edge
    Solution by separation of variables is

    p(0,0) = sum_{m=1, m odd}^{5nfty} 4/(m*pi)*sinh(5*m*pi)/sinh(10*m*pi)

    Already 2 terms give almost 30 d.p.
    This is the probability of first exit through the right-hand edge;
    we double this for the probability of first exit through either short edge.