Zhendong Wan's Approach to the Exact Solution of #7

F. Bornemann, October 22, 2004

The linear system at hand

n = 20000;
b = Table[0, {n}]; 
    b[[1]] = 1;
A = (SparseArray[{{ii_, ii_} -> Prime[ii]}, 
                  n] + (# + Transpose[#]) &)[
          SparseArray[
            Flatten[Table[{i, i + (2^j)} -> 1, {i, n - 1}, {j, 0, 
                  Log[2. , n - i]}]], n]];
\(diagonal = 
      Table[A[[i, i]], {i, n}]; 
    precond = #\/diagonal &;\)

Hadamard's bound of the determinant

Bnd = (∏_(i = 1)^n) A[[i, 
            i]];

Wan's iterative refinement algorithm (Algorithm 1 adopted to Problem 7)

Info[err_, step_, t_, flag_] := 
        If[Mod[step, 100] == 0 ∨ flag, 
          Print["#", step, ": Error ≤ ", 
            NumberForm[err, {4, 3}], "; Time: ", 
            NumberForm[TimeUsed[] - t, {10, 2}]]];
LinSolve[A_, b_, precond_, tol_] := 
        LinearSolve[A, b, 
          Method -> {Krylov, Method -> ConjugateGradient, 
              Preconditioner -> precond, 
              Tolerance -> tol}];
\(Wan[A_, b_, tol_] := 
      Module[{k = 0, r = b, den = 1, αα = {}, xx = {}, err, rn, 
          xn, x, α, t = TimeUsed[]}, 
        While[\((err = N[Norm[r, 2]\/den])\) > tol, \(k++\); 
          Info[err, k, t, False]; 
          xn = LinSolve[A, rn = N[r], precond, 1\/10\^15]; αα = 
            Append[αα, α = 
                Floor[Norm[rn, ∞]\/Norm[rn - A . xn, ∞]]]; 
          If[α < 2, Print["Insufficient accuracy"]; Break[]]; 
          xx = Append[
              xx, \((x = 
                    Floor[α\ xn])\)[[1]]]; den *= α; r = α\ r - A . x]; 
        Info[err, k, t, True]; 
        num = ∑\+\(i = 1\)\%k\((∏\+\(j = i + 1\)\%k αα[[j]])\)\ xx[[
                i]]; num\/den]\)

Check it out on 10.000 digits first

r = Wan[A, b, (10^(-10005))]; // Timing

N[r, 10002]

Now, let it run for 194.779 digits accuracy; accurately enough to reconstruct the rational solution

r = Wan[A, b, (1/(2.  (Bnd^2)))]; // Timing

The rational reconstruction step

RationalReconstruction[r_, Bnd_] := 
        Module[{a = Numerator[r], b = Denominator[r], u0 = 0, u1 = 1, v0 = 1, 
            v1 = 0}, 
          While[Abs[v1] ≤ Bnd, q = Quotient[a, b]; 
            b = a - q (a = b); u1 = u0 + q (u0 = u1); 
            v1 = v0 + q (v0 = v1)]; (u0/v0)];

sol = RationalReconstruction[r, Bnd]; // Timing

Numerator@sol

Denominator@sol

MaxMemoryUsed[]/(1024. ^2)