This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 3

Consider the sums (Subscript[c, k])=(∑_(j=0)^k)(((-1))^j) (({{k}, {j}}))/ζ(2 j+2) [✶]. Visualize the dependence of (Subscript[c, k]) on k for 0≤k≤(10^4). Assuming |(Subscript[c, k])|∼(k^α), what is an approximate value for α?

c[k_] := Sum[(-1)^j Binomial[k, j]/Zeta[2j + 2], {j, 0, k}]

While for small k , this straightforward definition works, for larger k, it becomes exponentially slower. We observe that the results are of the order -(10^(-smallInteger)).

$MaxExtraPrecision = 10^5;
dataT = {#, Timing[N[c[#], 20]]}& /@ {10, 50, 100, 200, 500, 1000, 2000}

So we have to implement an alternative way to get a numerical value for (Subscript[c, k]) with k≈(10^4). Let us look at the size of the largest integers occurring in the exact calculations. The largest binomial coefficients has about 3000 digits for k=(10^4), but already for k=2000, the denominator of the zeta function has more than 5100 digits. This suggests using high-precision floating point numbers instead of exact arithmetic in the following.

{Binomial[10000, 5000], Denominator[Zeta[2000]]} // N

Next, we observe that ζ(2n) approaches 1 exponentially fast for large n. By observing that (∑_(j=0)^k)(((-1))^j) (({{k}, {j}}))=0, we subtract this sum to obtain summands of the form largeBinomial×smallZeta, which gives much smaller values than the original small. The following graphic shows sizes of the k+1 summands for k=500 and k=1000. The size of the summands vary over many hundred orders of magnitude to give a small final result. Approximately the first k/2 terms contribute significantly to the final result.

Show[Graphics[
{Hue[#1], Table[Point[{j, Log[10, Abs[N[(-1)^j Binomial[#2, j]*
                (1/Zeta[2j + 2] - 1), 20]]]}], {j, 0, #2}]}& @@@
 {{0, 500}, {0.8, 1000}}], PlotRange -> All, Axes -> True]

Now, let us implement two further computational optimizations. For a given k, we have to generate all binomial coefficients. While this could be, of course, done with k+1 calls to Binomial, a recursive calculation of the binomial coefficients is hundreds of times faster for larger k.

(* list of all k + 1 binomial coefficients Binomial[k, j], j, 0, 1, …, k *)
binomialList[k_] := FoldList[Times, 1, (k - Range[0, k - 1])/(1 + Range[0, k - 1])]

{Table[Binomial[10000, k], {k, 0, 10000}]; // Timing,
 binomialList[10000]; // Timing}

The second computational optimization refers to the zeta function values. Because the zeta function term in the definition of the (Subscript[c, k]) is independent of k, we can once and for all calculate a list of zeta function values and reuse parts of it for a given k. So, we calculate the following list of zeta function values zetaList. We also include the factor (((-1))^j) and calculate the function values to a precision depending on the size of the largest binomial coefficient with which they will be later multiplied. The following calculation is needed one time and will take a few minutes.

BL = binomialList[10001] // N;

(zetaList = Table[(-1)^j (1/Zeta[N[2j + 2, (* use adaptive precision *)
                         Max[Ceiling[Log[10, BL[[j + 1]]] + 300], 7/10 j]]] - 1), 
                  {j, 10000, 0, -1}] // Reverse;) // Timing

So, we end with the following optimized definition for the (Subscript[c, k]).

cF[k_] := binomialList[k].Take[zetaList, k + 1]

Here is a quick check showing that c[k] and cF[k] give identical results.

(cF[#]/c[#] - 1)& /@ {10, 100, 1000}

cF is much faster c and we can calculate (Subscript[c, 10000]) in a few seconds.

{cF[1000] // Timing // N, cF[10000] // Timing // N}

Next, we calculate a larger set of values of the (Subscript[c, k]) and visualize the logarithm of their absolute value.

tab = Table[C6 = Timing[{k, cF[k]}]; C6[[2]], {k, 10000, 0, -100}];

ListPlot[Log @ Drop[Abs[tab], -1], PlotRange -> All,
         AxesLabel -> {Log[k], Log[Subscript[c, k]]}]

Modulo some small oscillations visible for larger k, we see a dependence of the form (Subscript[c, k])∼(k^α). Fitting the data points to such a curve yields α≈-2.

Fit[Select[Log @ Abs[tab], #[[1]] > 5&], {1, x}, x] // N

Here are the resulting more pronounced oscillations of (Subscript[c, k])/(k^(-2)).

ListPlot[{#1, #2 #1^2}& @@@ Abs[tab], PlotRange -> All]

Interestingly, the condition |(Subscript[c, k])|≪(k^(-3/4+ε)) is equivalent to the Riemann hypothesis [✶].