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