# Using Fold to show stages of Euler product formula

I would like to recursively replicate each stage of the Euler product proof. This does it:

ap[r1_, r2_, e_] :=  With[{r = Range[r1, r2]}, Pick[r, PrimeOmega[r], e]]
ap2 = ap[1, 10^3, 1];
β = 10^3;
a = Table[1/n, {n, β}]


$\quad\quad\quad\quad\quad\zeta(s)=1+\dfrac{1}{2^{s}}+\dfrac{1}{3^{s}}+\dfrac{1}{4^{s}}+\dfrac{1}{5^{s}}+\ldots$

b = Reverse@Complement[a, 1/ap2[[1]] a]


$\quad\quad\quad\quad\quad\dfrac{1}{2^{s}}\zeta(s)=\dfrac{1}{2^{s}}+\dfrac{1}{4^{s}}+\dfrac{1}{6^{s}}+\dfrac{1}{8^{s}}+\dfrac{1}{10^{s}}+\ldots$

$\quad\quad\quad\quad\quad\left(1-\dfrac{1}{2^{s}}\right)\zeta(s)=1+\dfrac{1}{3^{s}}+\dfrac{1}{5^{s}}+\dfrac{1}{7^{s}}+\dfrac{1}{9^{s}}+\dfrac{1}{11^{s}}+\dfrac{1}{13^{s}}+\ldots$

c1 = Reverse@Sort[Join[Complement[b, 1/ap2[[2]] b], -Complement[1/ap2[[2]] b, b]]];
c2 = DeleteCases[Table[If[Denominator[c1[[α]]] > β, 0, c1[[α]]], {α, Length@c1}], 0]


$\quad\quad\quad\quad\quad\dfrac{1}{3^{s}}\left(1-\dfrac{1}{2^{s}}\right)\zeta(s)=\dfrac{1}{3^{s}}+\dfrac{1}{9^{s}}+\dfrac{1}{15^{s}}+\dfrac{1}{21^{s}}+\dfrac{1}{27^{s}}+\dfrac{1}{33^{s}}+\ldots$

$\quad\quad\quad\quad\quad\left(1-\dfrac{1}{3^{s}}\right)\left(1-\dfrac{1}{2^{s}}\right)\zeta(s)=1+\dfrac{1}{5^{s}}+\dfrac{1}{7^{s}}+\dfrac{1}{11^{s}}+\dfrac{1}{13^{s}}+\dfrac{1}{17^{s}}+\ldots$

d1 = Reverse@Sort[Join[Complement[c2, 1/ap2[[3]] c2], -Complement[1/ap2[[3]] c2, c2]]];
d2 = DeleteCases[Table[If[Denominator[d1[[α]]] > β, 0, d1[[α]]], {α, Length@d1}], 0];


$\quad\quad\quad\quad\quad\text{etc}\ldots$

I can do it by naming each stage, but am finding difficulty coding it using Fold.

Note: The "minus Complement" is not necessary for the primes, but for other Euler products, it is necessary.

Hi Martin: Here is a partial solution using FoldList, along with a helper function fz defined as

fz[z_, p_Integer] :=
Join[{1},
Cases[Apply[List,Expand[Total[z(1-p^-s)]]/.Power[m_,Times[e_,s]]:>(m^-e)^-s],
_Integer^-s]
]


I use fz to multiply the current expansion z by 1-p^-s, for prime p, then mess with the expression to get a list of single integers to the power -s. The initial value for z is Table[1/k^s,{k,1,n}], where n is a maximum. In the question, you used n=1000.

Then use EulerProductRecursion[n] to get a table of successive iterations.

EulerProductRecursion[n_Integer] :=
Module[{z = Table[1/k^s, {k, 1, n}],
p = Prime[Range[PrimePi[Sqrt[n]]]]},
Print[p];
Total[FoldList[fz[#1, #2] &, z, p], {2}]]


The prime cutoff using PrimePi[Sqrt[n]] is not fundamental to the solution, and can be set to simply PrimePi[n], but this latter setting is boring for p>Sqrt[n].

• thanks for your answer - will try to have a play around with it :) – martin Dec 30 '14 at 4:25

I actually went for this in the end:

recursionNumber = 10; range = 100;
f[list_, a_] :=  list - Take[Riffle[Table[0, {x, plist[[a]] range}], list/plist[[a]],
plist[[a]]], range]
plist = Prime[Range[range]];
list = DeleteCases[Fold[f, Table[1/n, {n, range}], Range[recursionNumber]], 0];
Plus @@ Join[{1}, Rest@Table[DisplayForm[FractionBox[Numerator[list[[n]]],
Denominator[list[[n]]^s]]], {n, Length@list}]]


as it was more what I was after (should really have been a bit clearer in question), but hats off to KennyColnago for his excellent answer.