Richard Guy asked about a(4), the prime for which sum(1/p) exceeds 4.
I responded:
Rich> It would be straightforward to set up a distributed computation
Rich> to sum 1/p for p < 1.8e18, and locate the specific P that bumps
Rich> the sum over 4. In a few more years the resources will be
Rich> available to many of us as individuals. I've thought about
Rich> doing similar computations myself, but always get bogged down in
Rich> "As long as you're touching all those primes, what else do you
Rich> want to compute?". I also want sum 1/(p-1), product P, etc.
...
Rich> Odlyzko & friends have pushed exact computations of pi(N) up to
Rich> around N = 10^21 using a zeta function based formula, so perhaps
Rich> sum(1/p) could be tackled similarly. This doesn't require
Rich> touching every prime.
which prompted victor(a)idaccr.org (Victor S. Miller) to elaborate:
Ah, this piqued my attention. Since I was one of the "friends"
involved in the combinatorial formula (which was the one used for the
above quoted values -- the analytic formula seems much more difficult
to get to work in a practical range), I can outline what needs to be
done for an algorithm which will probably be O(x^(2/3) + epsilson)
(but there's a serious issue of precision -- see below) to calculate
sum_{p <= x} 1/p to sufficient accuracy. Once you have an algorithm
for that, you can find the crossover point you want by binary search
(or even better, interpolation search):
for k a positive integer, let p_k denote the k'th prime, and set
S(x,k) = sum { 1/n | 1 <= n <= x, and p_j does not divide n for j<=k }.
Then, we have the recursion, for k>0:
S(x,k) = S(x,k-1) - (1/p_k) S(x/p_k,k-1) (*)
and
S(x,0) = sum_{n<=x} 1/n
Now S(N,sqrt(N)) = sum_{ sqrt(N) <= p <= N} 1/p
[nit: ^^^^^^^ should be pi(sqrtN) ]
(this is tricky this sum will be very close to log(2), so the "extra"
amount is important).
The key to the algorithm is to use the recursion (*) as long as k>=k_0
(for some convenient k_0, like 5) and x >= N^(2/3). We can then
evaluate the values of S(x,a) that we have left by means for a tree
data structure (see Lagarias, Miller and Odlyzko, Math. Comp. 1984 (I
think)).
Since the number of summands involved in all of the intermediate sums
is <=N, it means that we need to know the initial values of
sum_{1 <= i <= k} 1/k (for all k <= N^2/3) to at least an accuracy of
(1/(2N)). However, the fact that the recursion (*) involves
subtraction might well make this numerically unstable and force a lot
more precision to be needed. Maybe someone will be motivated to make it work.
I've included Victor's full message, since I had several comments.
Re Future Prospects, for a(5), somewhere around 10^50.
If we have an N^(2/3) algorithm, this would only be ~ 10^33.
Although the sum total of the worlds computing resources today is perhaps
only 10^25 cycles/year, we can't say 10^33 is forever out of range.
Returning to more immediate possibilities:
In the specific case of sum(1/p), the precision demands for Victor's
recursion
S(x,k) = S(x,k-1) - (1/p_k) S(x/p_k,k-1) (*)
and base case S(x,0) = sum_{n<=x} 1/n
seem manageable. Expanding the minuend recursively converts the
recursion to
S(x,k) = S(x,0) - sum (1/p_k) S(x/p_k,k-1).
In this form, the depth of the recursion is at most log N, where
N is the largest prime in our sum. The width is at most pi(sqrtN). Any
given leaf at the bottom of the tree has an implicit maximum coefficient
of (sqrtN)^logN in the final answer, even ignoring the 1/p_k reduction.
This would mean an extra log(sqrtN^logN) bits of precision in the leaf
calculations, around (logN)^2 bits. A more careful account might cut
this to O(logN).
Computing the base case harmonic numbers to moderate precision isn't
difficult. Gosper has a nice accelerated series, and even Euler-Maclaurin
summation should work.
I mentioned wanting some other results, like sum 1/(p-1) and product P.
The recursion only works for powers of P: we could use slight variations
to sum P^2 or sqrt P, while redoing the error analysis a little.
The required fixup for summing log P seems easy enough. We could also
do sum 1/P^K for K = 2...100 or so, which should be enough to handle
sum 1/(p-C) for small real and complex C. 1/(p-1) = 1/p + 1/p^2 + ...,
and we can get stuff like sum 1/(p^2 + 2) using partial fractions.
sum sqrt(P+3) can be done with another power series, 1/p^(K+1/2).
I don't see how to do the sum of (logP)^2, or loglogP or 1/logP.
Predicting future computing capability is hard:
In the late 1800s, the mathmatician Jevons published a 20 digit number,
the product of two 10 digit primes known only to himself. Anticipating
RSA by almost a century, he confidently asserted that noone would ever
know the factors. Cole exhibited the factors of M67 = 2^67-1, 21 digits,
at a math conference in 1903.
Rich rcs(a)cs.arizona.edu
