Code from the Appendix

Supporting Interval Functions

mid[X_] := (Min[X] + Max[X])/2;
extremes[X_] := {Min[X], Max[X]};
diam[X_] := Max[X] - Min[X];

Subdivision of Intervals and Rectangles

subdivide1D[X_] := 
        Interval /@ {{Min[X], mid[X]}, {mid[X], Max[X]}};

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];

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)

Definition of ExpectedVisits (S. 169)

ExpectedVisits[ϵ_] := 
        1/ArithmeticGeometricMean[(Sqrt[1 - 4 ((((Subscript[p, EW]) + (Subscript[p, NS])))^2)]), (Sqrt[1 - 4 ((((Subscript[p, EW]) - (Subscript[p, NS])))^2)])] /. {(Subscript[p, EW]) ->  (Sqrt[((1/4) + ϵ) ((1/4) - ϵ)]), (Subscript[p, NS]) ->  (1/4)};

Sitzung S. 171

prec = 
      10010; ϵTry = ϵ /. 
        FindRoot[ExpectedVisits[ϵ] - 2, {ϵ, 0.06, 0.07}, 
          WorkingPrecision ->  prec, 
          AccuracyGoal -> prec]; ϵInterval = 
      IntervalBisection[
        ExpectedVisits[#] - 2 &, ϵTry + {-1, 1} 
            (10^(-prec + 5)), (10^(-prec))];
DigitsAgreeCount[ϵInterval]