Keith Briggs

This page was last modified 2013-05-12  


·maths notes « 
·ex libris 
·site map 

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:

Problem 1 

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



     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?


  1. Run a quick-and-dirty python program ( 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 (
     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 ( 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.


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

  More precisely:
     x(10)   = -0.736292698609618310775737171866391720451575491031024118015570
     y(10)   = -0.669642696363571375194376070967098466862227372928775421597411
     distance = 0.9952629194433541608903118094267216210294669227341543498032088

Problem 3

  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?


    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

Problem 4

  minimize  exp(sin(50x))


    x=-0.0244030796943785 y=2.10612427155358

     Random search followed by gradient optimization.

     Check with gp:
     41 100
     x=-0.0244030796962585; y=0.210612427157718

Problem 5

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

    min sup-norm error = 0.2143352346...

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

  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?



    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}?



    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.
  At what time does u at (1,1) reach 1?



   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

    a = 0.78593367197930...

    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?

    = 0.00000038375879792512261034071331862048...

    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.

This website uses no cookies. This page was last modified 2013-05-12 10:17 by Keith Briggs.