Chapter 7
Too Large to Be Easy, 
Too Small to Be Hard

Code from the Appendix

Count of Digit Agreement and Form for Pretty Interval Output

DigitsAgreeCount[a_, 
          b_] := (prec = 
            Ceiling@Min[Precision /@ {a, b}]; {{ad, ae}, {bd, 
                be}} = RealDigits[#, 10, prec] & /@ {a, b}; 
          If[ae ≠ be ∨ a b ≤ 0, Return[0]]; 
          If[ad == bd, Return@Length[ad]]; {{com}} = 
            Position[MapThread[Equal, {ad, bd}], False, 1, 1] - 1; 
          com);
DigitsAgreeCount[Interval[{a_, b_}]] := DigitsAgreeCount[a, b]; 
    IntervalForm[
        Interval[{a_, 
            b_}]] := (If[(com = DigitsAgreeCount[a, b]) == 0, 
          Return@Interval[{a, b}]]; 
        start = Sign[a] N[
              FromDigits@{ad[[
                    Range@com]], 1}, com]; {low, 
            up} = SequenceForm @@ Take[#, {com + 1, prec}] & /@ {ad, bd}; 
        If[ae == 0, start /= 10; ae++]; 
        SequenceForm[
          DisplayForm@SubsuperscriptBox[NumberForm@start, low, up], 
          If[ae ≠ 1, 
            Sequence @@ {" × ", 
                DisplayForm@
                  SuperscriptBox[10, ae - 1]}, ""]])

§7.1 A First Look: A Linear System

n = 20000; b = Table[0, {n}]; 
    b[[1]] = 1;
A = SparseArray[{{i_, i_} -> Prime[i]}, 
                n] + (# + Transpose[#]) &@
          SparseArray[
            Flatten@Table[{i, i + (2^j)} -> 1, {i, n - 1}, {j, 0, 
                  Log[2. , n - i]}], n];

§7.3.1 Preconditioning

diagonal = 
      Table[A[[i, i]], {i, n}]; 
    prec = 100; b = SetPrecision[b, prec + 5];
\(\(x = 
        LinearSolve[A, b, 
          Method -> {Krylov, Method -> ConjugateGradient, 
              Preconditioner -> \((#\/diagonal &)\), 
              Tolerance -> 10\^\(\(-prec\) - 1\)}];\)\)
N[x[[1]], 100]

§7.4.2 Validation: Interval Arithmetic

x[[1]] + 
        Interval[{-1, 1}] Norm[b - A . Interval /@ x] // 
      IntervalForm

Proof of Lemma 7.1

n = 1142;
A = SparseArray[{{i_, i_} -> Prime[i]}, 
                n] + (# + Transpose[#]) &@
          SparseArray[
            Flatten@Table[{i, i + (2^j)} -> 1, {i, n - 1}, {j, 0, 
                  Log[2. , n - i]}], n];

{λ, V} = 
        Eigensystem[Normal@N@A];
r = (λMin = λ[[
                  n]]) (x = 
                Interval /@ 
                  V[[n, All]]) - 
          x . A;
(λMin += Interval[{-1, 1}] Norm[r]/Norm[x]) // 
      IntervalForm

λ1 = Prime[n + 1] - 100; α2 = 
      11 × 
        100; (λFMin = ((2 ((λ0 = Interval[Min@λMin]) λ1 - α2))/(λ0 + λ1 + (Sqrt[4×α2 + (((λ1 - λ0))^2)])))) // IntervalForm

§7.5 Solving it Exactly

Here we present a straightforward implementation of the mathematics exposed in §7.5. Because Mathematica cannot be told to use decisive word-sized integer arithmetic, the run time is by no means comparable to LinBox. However, the implementation may serve educational purposes.

The Linear System at Hand

n = 40;
b = Table[0, {n}]; 
    b[[1]] = 1;
A = SparseArray[{{i_, i_} -> Prime[i]}, 
                n] + (# + Transpose[#]) &@
          SparseArray[
            Flatten@Table[{i, i + (2^j)} -> 1, {i, n - 1}, {j, 0, 
                  Log[2. , n - i]}], n];

§7.5.1 The Congruential Approach

Cramer's Rule

x1[A_] := 
        Determinant[Take[A, {2, -1}, {2, -1}]]/
          Determinant[A];

The Chinese Remaindering

Needs["NumberTheory`NumberTheoryFunctions`"];
WorkPrimes[A_] := 
        Module[{HadamardBound = Product[A[[i, i]], {i, Length@A}], 
            work = {}, prod = 1, k = PrimePi[(2^30)] + 1}, 
          While[prod < HadamardBound, work = {work, p = Prime[k++]}; 
            prod *= p]; Flatten@work];
Determinant[A_] := 
        Module[{primeList = WorkPrimes[A], detList}, 
          detList = Determinant[A, Modulus ->  #] & /@ primeList; 
          ChineseRemainder[detList, primeList]];

Wiedemann's Randomized Algorithm for Modular Determinants

Determinant[A_, 
          opts___Rule] := (m = Modulus /. opts; 
          Module[{n = Length@A}, 
            Mod[(((-1))^n) First@MinPol[A, n], m]]);
MinPol[A_, minDegree_] := 
        Module[{n = Length@A, pol = {}, u, v, b}, 
          While[Length@pol ≤ minDegree, {u, v} = 
              Table[Random[Integer, {0, m - 1}], {2}, {n}]; 
            a = Mod[NestList[u . (b = Mod[A . b, m]) &, u . (b = v), 
                  2 n - 1], m]; pol = LinearRecurrence[a, n]]; 
          pol];

Berlekamp-Massey algorithm as given in [KS91,p. 31]

LinearRecurrence[a_, d_] := 
        Module[{Λ = {1}, Σ = {}, 
            l = 0, δ = 1, r, δr, Λ1}, 
          For[r = 1, r ≤ 2 d, r++, 
            PrependTo[Σ, 0]; δr = 
              Take[a, {r, 
                    r - Length@Λ + 
                      1, -1}] . Λ; 
            If[(δr = Mod[δr, m]) > 0, Λ1 = 
                PadRight[Λ, 
                    Length@Σ] - δr × Σ/δ; 
              If[2 l < r, 
                l = r - 
                    l; δ = δr; Σ = Λ]; Λ = PolynomialMod[Λ1, m]]]; 
          Reverse@Λ];

The Solution

sol1 = x1[A]

§7.5.1 The p-adic Approach

The Iteration: Hensel Lifting

pAdicLinearSolve[A_, b_, p_, k_] := 
        Module[{lud = LUDecomposition[A, Modulus -> p], x, j, 
            r, Δx}, 
          x = LUBackSubstitution[lud, b, Modulus ->  p]; 
          For[j = 1, j < k, j++, 
            r = (b - A . x)/(p^j); Δx = 
              LUBackSubstitution[lud, r, Modulus ->  p]; 
            x = x + (p^j) Δx]; x];

The Rational Number Reconstruction as given in [vzGG99,§5.10]

RationalReconstruction[f_, m_] := 
      Module[{r0 = m, r = f, t0 = 0, t = 1}, 
        While[Abs[r] ≥  (Sqrt[m]), 
t = 
            t0 - Quotient[r0, r] (t0 = t); r = Mod[r0, (r0 = r)]]; 
        (r/t)]

The  Number of Lifting Steps

pAdicExponent[A_, p_] := 
        Ceiling[2 
             Log[p, 2 Product[
                  A[[i, i]], {i, 
                    Length@A}]]];

The Solution

p = 1125899906842597;
k = pAdicExponent[A, p];
x = pAdicLinearSolve[A, b, p, k];
sol2 = 
      RationalReconstruction[x[[1]], 
        (p^k)]

A Check

sol1 === sol2