Problem 7:

Let be the 20,000 x 20,000 matrix whose entries are zero everywhere except for the primes 2, 3, 5, 7, ..., 224737 along the main diagonal and the number 1 in all the positions with ,2,4,8,...,16384. What is the (1,1) entry of A^(-1) ?

Solution:

This is equivalent to computing after solving , where . The system is diagonally dominant, and we use an iterative method: , where D is the diagonal of (but without actually constructing the matrices). Gonnet used a Darwin program, equivalent to the following in Maple.

> b:= Array(1..20000,datatype=float[8]): b[1]:= 1:
x:= Array(1..20000,datatype=float[8]):
primes:= Array(1..20000,datatype=float[8]):
for i from 1 to 20000 do primes[i]:= ithprime(i) od:

>    solvproc:= proc(N)
local iter, toterr, i, rhs, ij, newxi;
global x;
for i from 1 to N do x[i]:= 0 od:
for iter to 40 do
    toterr := 0;
    for i to N do
    rhs := b[i];
      ij := 1;
    while ij < N do
      if i-ij > 0 then rhs := rhs - x[i-ij] fi;
      if i+ij <= N then rhs := rhs - x[i+ij] fi;
      ij := 2*ij
    od;
    newxi := rhs/primes[i];
    toterr := toterr + abs(newxi-x[i]);
    x[i] := newxi
  od;
    if toterr=0 then break fi;
    lprint(toterr);
od:
x[1];
end;