C.5 Mathematica Code

C.5.1 Utilities

Kronecker Tensor Product of Matrices

SetAttributes[KroneckerTensor, OneIdentity]; 
    KroneckerTensor[u_?MatrixQ, v_?MatrixQ] := 
      Module[{w = Outer[Times, u, v]}, 
        Partition[
          Flatten[Transpose[
              w, {1, 3, 2, 4}]], Dimensions[
                w][[2]] Dimensions[
                w][[4]]]];
KroneckerTensor[u_, v_, w__] := 
        Fold[KroneckerTensor, u, {v, w}];
CircleTimes = KroneckerTensor;

Supporting Interval Functions

mid[X_] := (Min[X] + Max[X])/2;
diam[X_] := Max[X] - Min[X];
diam[{X__Interval}]  :=  
        Max[diam  /@  {X}];
extremes[X_] := {Min[X], Max[X]};
hull[X_] := Interval[extremes[X]];
MidRad[x_, r_] := 
        x + Interval[{-1, 1} ] r;
IntervalMin[{Interval[{a_,  b_}],  
            Interval[{c_,  d_}]}]   :=   Interval[{Min[a, c],  
            Min[b, d]}];
IntervalMin[{}]   :=  ∞;

Subdivision of Intervals and Rectangles

subdivide1D[X_] := 
        Interval /@ {{Min[X], mid[X]}, {mid[X], 
              Max[X]}};
subdivide2D[{X_, Y_}] := 
        Distribute[{subdivide1D[X], subdivide1D[Y]}, List];

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]}, ""]])

C.5.2 Problem-Dependent Functions and Routines

Reliable Photon Trajectory (Chapter 2)

Options[
          ReliableTrajectory]  :=  {StartIntervalPrecision  -> 
            Automatic, 
          AccuracyGoal  ->  12};
\(\(ReliableTrajectory[p_, v_, tMax_, \ opts___Rule] := 
        Module[{error = ∞, ϵ, s0, \ s, P, V, t, Trem, 
            path, S, 
            T}, 
{ϵ, s, 
              s0}\  = \({10. ^\((\(-AccuracyGoal\))\), AccuracyGoal, \ 
                  StartIntervalPrecision} /. {opts}\) /. 
              Options[ReliableTrajectory]; 
If[
            s0 === Automatic, \ s0\  = \ s]; 
s = s0; 
          wp = Max[17, s + 2]; 
While[
            error > ϵ, 
P = 
              N[\((MidRad[#, 10\^\(-s\)] &)\) /@ p, 
                wp]; 
V = 
              N[\((MidRad[#, 10\^\(-s\)] &)\) /@ v, 
                wp]; 
path = {P}; 
Trem = 
              Interval[tMax]; 
While[Trem > 0, 
              M = Round[P + 2  V/3\ ]; 
S = 
                t /. \((Solve[\((P + t\ V - M)\) . \((P + 
                              t\ V - M)\) == 1/9, 
                      t])\); 
If[
                FreeQ[S, 
                  Power[Interval[{_?Negative, \ _?
                          Positive}], _]], 
\ \ \ \ \ \ \ \ \ \ \ T = IntervalMin[Cases[S, _?Positive]], 
                Break[]]; 
Which[
T ≤ Trem, P += T\ V; V = H[P - M] . V; \ 
                Trem -= T, 
T > 
                    Trem\  && \ \ Trem ≥ 2/3, P += 2  V/3; \ 
                Trem -= 2/3, 
T > Trem\  && \ Trem < 2/3, 
                P += Trem\ V; Trem = 0, 
True, \ 
                Break[]]; 
AppendTo[path, 
                P]; 
If[
                Precision\ [\ {Trem, \ P, \ V, \ T}] < \ ag, \ 
                Break[]]]; 
wp\  = 
              Max[17, \((\(++s\))\) + 2]; 
error = 
              diam[P + 
                  Table[MidRad[\(-Max[Abs[Trem]]\), 
                      Max[Abs[Trem]]], {2}]]]; 
Print[
            StringForm[\*""Initial condition interval radius is (10^(-``))."", s]]; 
path];\)\)

Return Probability (Chapter 6)

Matrices[
          n_] := Matrices[n] = (m = 2 n + 1; (Subscript[p, E]) = (1/4) + ϵ0; 
            (Subscript[p, W]) = (1/4) - ϵ0; (Subscript[p, N]) = (1/4); (Subscript[p, S]) = (1/4); 
            Id = SparseArray[{i_, i_} -> 1, {m, 
                  m}]; 
PEW = 
              SparseArray[{{i_, j_} /; j == i + 1 -> 
                    (Subscript[p, E]), {i_, j_} /; j == i - 1 -> (Subscript[p, W])}, {m, 
                  m}]; 
PNS = 
              SparseArray[{{i_, j_} /; j == i + 1 -> 
                    (Subscript[p, N]), {i_, j_} /; j == i - 1 -> (Subscript[p, S])}, {m, 
                  m}]; {PEW⊗Id + Id⊗PNS, 
              Id⊗Id});
ReturnProbability[ϵ_Real, 
          n_Integer] := ({A, Id} = Matrices[n]; 
          Block[{ϵ0 = ϵ}, m = 2 n + 1; ctr = n m + n + 1; 
            r = A[[All, ctr]]; 
            A[[All, ctr]] = 0; 
            q = LinearSolve[Id - A, r]; 
            q[[ctr]]]);

Multiplication with the Matrix (Subscript[A, n]) (Chapter 7)

n = 20000;
BitLength := Developer`BitLength;
indices[i_] := indices[i] = 
        i + Join[
            2^(Range@BitLength[n - i] - 
                  1), -2^(Range@BitLength[i - 1] - 1)]
diagonal = Prime[Range@n];
A[x_] := 
        diagonal* x + 
          MapIndexed[
            Plus @@ x[[
                  indices[#2[[1]]]]] &, x];

Proof of the Rigorous Bound for (Subscript[λ, min])((Subscript[A, 1142])) (Chapter 7)

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

Interval Theta Function (Chapter 8)

SetAttributes[{θ0, θ1}, 
        Listable];
θ0[q_, K_] := 
        2  (q^(1/
                4)) (∑_(k = 0)^K)((((-1))^k)/(2  k + 1)) 
              (q^(k (k + 1)));
θ1[
        q_, K_] := ((q^(-3/4))/2) (∑_(k = 0)^K)(((-1))^k) (2  k + 1) (q^(k (k + 1))); θ[
        q_Interval, K_Integer] := hull[θ0[q, {K - 1, K}]];
Derivative[1, 0][θ][q_Interval, K_Integer]:=hull[θ1[q, {K - 1, K}]];

C.5.3 General Functions and Routines

Two-Dimensional Minimization

LowestCriticalPoint[f_,  {x_,  a_,  b_},  {y_,  c_,  d_}, 
          upperbound_,  
          tol_] := (rects = 
            N[{ {Interval[{a, b}], 
                  Interval[{c, d}]}}]; 
fcn[{xx_, 
                yy_}]  :=  
            f /. { x ->  xx, 
                y ->  yy}; 
gradf[{xx_, yy_}] := 
            Evaluate[{D[f, x], D[f, y]} /. {x -> xx,  
                  y ->  yy}]; 
{low, 
              upp} = {-∞, 
              upperbound}; 
While[(upp - low) > 
              tol, 
rects = 
              Join  @@  subdivide2D  /@  
                  rects; 
fvals = 
              fcn  /@  rects; 
upp = 
              Min[upp, Min[Max /@ fvals]]; 
pos = 
              Flatten[Position[
                  Min /@ fvals, _?(# ≤ 
                          upp  &)]]; 
rects = 
              rects[[pos]];  
            fvals = fvals[[
                pos]]; 
pos =  
              Flatten[Position[
                  Apply[And, IntervalMemberQ[gradf /@  rects, 0],  {1}],  
                  True]]; 
rects = 
              rects[[pos]]; 
            low = Min[
                fvals[[
                  pos]]]]; 
{{low, upp}, 
            Map[mid,  rects,  {2}]});

Interval Arithmetic-Geometric Mean

AGMStep[{a_, b_}] := {(Sqrt[a b]), 
          (1/2) (a + b)};
EnclosingAGM[{a_Real, b_Real}] := 
        Interval@
          FixedPoint[
            extremes@AGMStep[Interval /@ #] &, {a, 
              b}];
Unprotect[ArithmeticGeometricMean];
ArithmeticGeometricMean[A_Interval, B_Interval] := 
        Block[{Experimental`$EqualTolerance = 0, 
            Experimental`$SameQTolerance = 0}, 
          Interval@
            extremes[EnclosingAGM /@ {Min /@ #, Max /@ #} &@{A, 
                  B}]];
Protect[ArithmeticGeometricMean];

Interval Bisection Method

IntervalBisection[f_, {a_, b_}, tol_] := (X = Interval[{a, b}]; 
        pos = {1}; 
        While[Max[diam /@ X] > tol, 
          pos = Flatten@
              Position[
                f /@ (X = 
                      Join @@ subdivide1D /@ 
                          X), _?(IntervalMemberQ[#, 0] &)]; 
          X = X[[pos]]]; 
        IntervalUnion @@ X)

Rational Reconstruction

RationalReconstruction[p_, q_, Bnd_] := 
        Module[{a = p, b = q, u0 = 0, u1 = 1, v0 = 1, v1 = 0}, 
          While[Abs[v1] ≤ Bnd, qt = Quotient[a, b]; 
            b = a - qt (a = b); u1 = u0 + qt (u0 = u1); 
            v1 = v0 + qt (v0 = v1)]; (u0/v0)];

Interval Newton Iteration

IntervalNewton[f_, {a_, b_}]:=Block[{Experimental`$EqualTolerance = 0},X1 = Interval[{a, b}];X = Interval[a];While[X =!= X1,X = X1;x = Interval[mid[X]];X1=IntervalIntersection[x-((f[x])/(f'[X])),X]];X]