Fastest algorithm to compute the sum of primes? Can anyone help me with references to the current fastest algorithms for counting the exact sum of primes less than some number n?  I'm specifically curious about the best case running times, of course.  I'm pretty familiar with the various fast algorithms for prime counting, but I'm having a harder time tracking down sum of prime algorithms...
 A: This has been great.  Some nice answers here - I'll accept one shortly.
I was looking for a lit review because I have implemented an approach that runs in $O(n^\frac{2}{3} \log n)$ time and $O(n^\frac{1}{3} \log n)$ space for summing primes to non-negative integer powers, and I wanted to see how it compared.  I wasn't exactly ready to write up what I'd done, but I guess I'll take a stab.  For the curious:
PART 1 - The General Identity
Start with the following identity
$\displaystyle\sum_{j=2}^n\frac{\Lambda(j)}{\log j}j^a = \sum_{j=2}^nj^a - \frac{1}{2}\sum_{j=2}^n\sum_{k=2}^{\frac{n}{j}}j^a \cdot k^a+ \frac{1}{3}\sum_{j=2}^n\sum_{k=2}^{\frac{n}{j}}\sum_{l=2}^{\frac{n}{jk}}j^a \cdot k^a \cdot l^a - \frac{1}{4}...$
where $\Lambda(n)$ is the Mangoldt function.  This is a generalization of Linnik's identity found here.  For ease of notation, using
$D_{a,k}(n) = \displaystyle\sum_{j=2}^n j^a \cdot D_{a,k-1}(\frac{n}{j})$
and
$D_{a,0}(n) = 1$
rewrite the nested sums on the right to give us
$\displaystyle\sum_{j=2}^n\frac{\Lambda(j)}{\log j}j^a = \sum_{j=1}^{\log_2 n} \frac{-1^{j+1}}{j}D_{a,j}(n)$
We evaluate $\log_2 n$ terms because $D_{a,k}(n) = 0$ when $n < 2^k$,
This counts a kind of sum of prime powers.  We can count prime sums by a kind of moebius inversion
$\displaystyle\sum_{p \leq n} p^a = \sum_{j=1}^{\log_2 n} \sum_{k=1}^{\log_2 (n^\frac{1}{j})} \frac{-1^{k+1}}{j \cdot k}\mu(j) D_{a \cdot j, k}(n^{\frac{1}{j}})$
Thus, if we can find fast methods of computing $D_{a,k}(n)$, we'll have fast ways of computing $\sum_{p \leq n} p^a$
In Mathematica, this is
DD[a_, k_, n_] := Sum[j^a DD[a, k - 1, n/j], {j, 2, n}]
DD[a_, 1, n_] := Sum[j^a, {j, 2, n}]
SumPrimes[a_, n_] :=
Sum[(-1)^(k + 1)/(j k) MoebiusMu[j] DD[j a, k, n^(1/j)], {j, 1,
Log[2, n]}, {k, 1, Log[2, (n^(1/j))]}]
SumPrimes[0,n] will return the count of primes, SumPrimes[1,n] will return the sum of primes less than n, and so on.
PART 2 - The Specific Algorithm
Here's my fastest approach for calculating $D_{a,k}(n)$  It is inspired by this method for calculating the Mertens function.  (If you find the following writeup interesting but rushed, I've written about my approach with somewhat more detail here and, for the special case of just prime counting, here.  A C++ implementation is here)
Bear with me; this is a bit involved.  It requires two observations.
First, we can use sieving to calculate values of $D_{a,k}(n)$ up to $n^\frac{2}{3}$ in roughly $O(n^\frac{2}{3} \log n)$ time and $O(n^\frac{1}{3} \log n)$ space with the following approach.
We want to sieve numbers up to $n^\frac{2}{3}$ in such a way that we have their full power signature, with $ n = {p_1} ^ {a_1} \cdot {p_2} ^ {a_2} \cdot {p_3} ^ {a_3} \cdot ... $.  With this, the normal count of divisor function is $d_k(n) = \binom {{a_1} + k - 1} {a_1} \cdot \binom {{a_2} + k - 1} {a_2} \cdot \binom {{a_3} + k - 1} {a_3} \cdot ...$  The strict count of divisors function, a count of divisor function excluding 1 as a factor, is, in turn, $d_k'(n) = \sum_{j=0}^k -1^{k-j} \binom {k}{j} d_j(n)$.  With the strict count of divisors function, we can say $ D_{a,k}(n) = D_{a,k}(n-1) + d_k'(n) \cdot n^a $ .  Thus, if we sieve all the numbers from 1 to $n^\frac{2}{3}$, for each number we can compute its values for $d_k(n)$, from which we can calculate $d_k'(n)$, which in turn lets us keep a running total of $D_{a,k}(n)$.  We will segment our sieve in blocks of size $n^\frac{1}{3}$ to stay in our memory bound.
The second observation is that, using the following identity, we can express $D_{a,k}(n)$ as sums that rely only on $D_{a,k}(j)$ where $j < n^\frac{2}{3}$ and of arbitrary values of $D_{a,1}(j) = \sum_{k=2}^j k^a$, which can generally be computed in constant or log time for any values of j for the values of a that we're working with.
If we let $d_{a,k}(n) = D_{a,k}(n) - D_{a,k}(n-1)$, we start with this combinatorial identity, whose derivation I have written up elsewhere but don't have the room to show:
$D_{a,k}(n) = $ $\displaystyle\sum_{j=t+1}^n d_{a,1}(j) \cdot D_{a, k-1}(\frac{n}{j}) $ $\displaystyle + \sum_{j=2}^t d_{a,k-1}(j) \cdot D_{a,1}( \frac{n}{j}) $
$+ \displaystyle \sum_{j=2}^t \sum_{s=\frac{t}{j} + 1}^\frac{n}{j} \sum_{m=1}^{k-2} d_{a,1}(s) \cdot d_{a,m}(j) \cdot D_{a, k-m-1}( \frac{n}{js} )$
For this identity, t is some value such that $1 < t < n$.  If you look through the sums closely, you should see that the right hand side doesn't rely on any values of $d_{a,k}(j)$ where $j > t$ and $D_{a,k}(j)$ where $ j > \frac{n}{t}$ except when $k = 1$, as desired.
Even if we set $t=n^\frac{2}{3}$ and sieve values of $D_{a,k}(j)$ up to $n^\frac{2}{3}$, these sums are too large to be computed in our time bound.  Instead, we have to take advantage of certain symmetries to reduce the calculations in these sums.  I'm going to hand wave over the process leading to this - you can find a description of it and much else here, but, taking one last leap, the following Mathematica code takes the identity we were just talking about and evolves it further as the function labeled DDFast - remember, when tracing through the following code, that in a full implementation any instances of DD and d in DDFast would be looked up thanks to sieving and caching within our time and space bounds, so what follows isn't an accurate reflection of its execution time, just the accurate mechanics of DDFast.
DD[A_, k_, n_] := Sum[j^A DD[A, k - 1, n/j], {j, 2, n}]
DD[A_, 1, n_] := Sum[j^A, {j, 2, n}]
d[ A_, k_, n_ ] := DD[ A, k, n ] - DD[ A, k, n - 1]
rng[ 0, start_, end_ ] := Floor[end] - (start  - 1)
rng[ 1, start_, end_ ] :=
Floor[end] (Floor[end] + 1)/2 - (start - 1) start/2
rng[ 2, start_, end_ ] :=
Floor[end] (Floor[end] + 1) (2 Floor[end] + 1)/
6 - (start - 1) start (2 start - 1)/6
rng[ A_, start_, end_ ] := Sum[m^A, {m, start, end}]
DDFast[A_, 1, n_] := rng[A, 2, n]
DDFast[A_, k_, n_] :=
Sum[j^A DD[A, k - 1, n/j], {j, Floor[n^(1/3)] + 1, n^(1/2)}] +
Sum[rng[A, Floor[n/(j + 1)] + 1, n/j] DD[A, k - 1, j], {j, 1,
n/Floor[n^(1/2)] - 1}] +
Sum[d[A, k - 1, j] rng[A, 2, n/j], {j, 2, n^(1/3)}] +
Sum[s^A d[A, m, j] DD[A, k - m - 1, n/(j s)], {j, 2, n^(1/3)}, {s,
Floor[Floor[n^(1/3)]/j] + 1, Floor[n/j]^(1/2)}, {m, 1, k - 2}] +
Sum[(rng[A, Floor[n/(j (s + 1))] + 1, n/(j s)]) (Sum[
d[A, m, j] DD[A, k - m - 1, s], {m, 1, k - 2}]), {j, 2,
n^(1/3)}, {s, 1, Floor[n/j]/Floor[Floor[n/j]^(1/2)] - 1}]
SumPrimesFast[A_, n_] :=
Sum[(-1)^(k + 1)/(j k) MoebiusMu[j] DDFast[j A, k, n^(1/j)], {j, 1,
Log[2, n]}, {k, 1, Log[2, (n^(1/j))]}]
The function labeled DDFast, which is really the key here, can be computed in $O(n^\frac{2}{3})$ time as long as all of the functions it references have been cached.
So there we have it.  If you compute SumPrimesFast and, in turn, DDFast, but interleave the computation of those sums with the sieving process described above, you end up with something like an $O(n^\frac{2}{3} \log n)$ time and $O(n^\frac{1}{3} \log n)$ space algorithm that can count prime sums to non-negative, integer powers, including as a prime counting function for $a = 0$.
It's worth taking a look, too, at rng[ A_, start_, end_ ] := Sum[m^A, {m, start, end}] - the reason this algorithm works for non-negative integers is that rng can be computed in constant time for non-negative, integer values with Faulhaber's Formula.  It's possible the algorithm could work for other powers if rng can be usefully extended to other powers with sufficient speed and precision.
I have an implementation of this algorithm in C++ here.  It has a few C math precision issues I haven't tracked down, so it's often off by a 2 or 3, but it does essentially implement what I've described here, and might be useful for stepping through with breakpoints.
A somewhat more fleshed out description of this can be found here.  Several details that I've hand-waved over here are described there.  Even better written up is a description of this as a prime counting algorithm, here.  It does a better job showing a few key derivations.
A: This is a difficult problem. I asked about it here on math.se and here on cstheory. In both cases my question was somewhat broader: allowing sums over different exponents rather than just 1.  In the second link I also allowed interactive proofs.
I noted that for exponent 0 there are efficient algorithms (superior to enumerating primes): the various analytical or combinatorial $\pi(x)$ algorithms. But no one was able to suggest an efficient means to solve the problem for any exponent.  An interactive proof was proposed that allows a person to check a proposed proof in linear time, but the confidence in the proof is asymptotically 0% against a cunning adversary.
I am not aware of any hardness results, though, so this seems to be wide-open.
A: Deléglise-Dusart-Roblot [1] give an algorithm which determines $\pi(x,k,l)$, the number of primes up to $x$ that are congruent to $l$ modulo $k,$ in time $O(x^{2/3}/\log^2x).$ Using this algorithm to count the number of primes in all residue classes $k<2\log x$ takes
$$1+\sum_{p<2\log x}(p-2)\sim\frac{2\log^2x}{\log\log x}$$
invocations of Deléglise-Dusart-Roblot for a total of $O(x^{2/3}/\log\log x)$ time.
This allows one to determine the value of $\sum_{p\le x}p$ mod all primes up to $2\log x$ and hence, by the Prime Number Theorem and Chinese Remainder Theorem, the value of the sum mod $\exp(\vartheta(2\log x))=x^2(1+o(1)).$ Together with bounds on the value of $\sum_{p\le x}p$ [2], this allows the computation of the sum.
Note that the primes slightly beyond $2\log x$ may be required depending on the value of $\vartheta(2\log x).$ Practically speaking, except for $x$ tiny, $2\log x+\log x/\log\log x$ suffices. This does not change the asymptotics.
I do not know if it is possible to modify the Lagarias-Odlyzko analytic formula [3] to count in residue classes. If so, this would allow an $O(x^{1/2+o(1)})$ algorithm.
References
[1] Marc Deléglise, Pierre Dusart, and Xavier-François Roblot, Counting primes in residue classes, Mathematics of Computation 73:247 (2004), pp. 1565-1575. doi 10.1.1.100.779
[2] Nilotpal Kanti Sinha, On the asymptotic expansion of the sum of the first n primes (2010).
[3] J. C. Lagarias and A. M. Odlyzko, Computing $\pi(x)$: An analytic method, Journal of Algorithms 8 (1987), pp. 173-191.
A: Edit Nov 22: Changed the condition of the test function $\Phi$ somewhat to simplify my argument (removed one sum in the identity below).
Although I like Charles' answer and its application of the Chinese remainder theorem, let me give you a different take on the problem. It is possible to use an analytic method more directly (similar to the Lagarias-Odlyzko method), just using the Riemann zeta-function and not any theory of the Dirichlet L-functions (or primes in arithmetical sequences). This method can treat rather arbitrary sums
$$ \sum_{ p < x } f(p) $$ for smooth functions $f$.
The idea is to consider the identity 
$$\sum_{  p < x } f(p) =\sum_{p} f(p)\Phi \left( \frac {p-x} {\sqrt x} \right )
- \sum_{x< p < x + \sqrt x } f(p) \Phi \left(\frac {p-x} {\sqrt x} \right), $$
where $\Phi(x)$ is smooth test function such that $\Phi(x)=1$ for $x<0$ and $\Phi(x)=0$ for $x>1$. It is clear that the last sum can be calculated in $O(  x^{1/2} )$ time (say sieving out the primes in the interval, and explicit calculation). The idea is now to use the Mellin-transform identity
$$ \sum_{p}f(p)\Phi \left( \frac {p-x} {\sqrt x} \right )= \frac 1 {2 \pi i} \int_{c-\infty i}^{c+\infty i} \sum_{p} p^{-s} \int_0^\infty \Phi \left(\frac{y-x}{\sqrt x} \right)f(y) y^{s-1} dyds.  (c>1)$$
The point here is that the Mellin transform $\int_0^\infty \Phi(\frac{y-x}{\sqrt x})f(y) y^{s-1}dy$ is small when $|\Im(s)|>x^{1/2+\epsilon}$ and $\Re(s)=c$ for any $\epsilon>0$, so that part of the integral can be discarded. The Dirichlet series $\sum_{p} p^{-s}$ can be expressed in terms of the logarithm of the zeta function for the arguments $ks$ (summing over $k$). Now the Odlyzko-Schönhage algorithm (Quite nice algorithm, uses fast fourier transform) allows us to calculate the values of the Riemann zeta-function for say $s=c+it$ for $ |t| < x^{1/2+\epsilon}$ in $O(x^{1/2+\epsilon+o(1)})$ time. We also remark that the Weil explicit formula can be used instead of this complex integral (then the zeros of the Riemann zeta-function needs to be calculated fast, but this can also  be done by the Odlyzko-Schönhage algorithm). This means that the total time to calculate $\sum_{  p < x } f(p)$ will be $O(x^{1/2+\epsilon})$ for any $\epsilon>0$.
Note that the same argument applies when primes in arithmetical progression are concerned (see my comment on Charles answer). Since the Odlyzko-Schönhage algorithm also holds for the Dirichlet L-functions, this case can be treated in the same way.
A: I'll put in a plug for my original paper with Lagarias and Odlyzko, as well as a recent paper by Bach, Klyve and Sorenson: http://www.ams.org/journals/mcom/2009-78-268/S0025-5718-09-02249-2/home.html Computing prime harmonic sums Math. Comp. 78 (2009), 2283-2305.  Although the algorithms of Lagarias and Odlyzko (with the extra algorithm of Odlyzko-Schonhage for computing good approximations to a bunch of values of $\zeta(s)$) are asymptotically the best, having complexity $O(x^{1/2 + \epsilon})$ the combinatorial algorithms will probably work best for any reasonable ranges.  Though one should look at J. Buethe, J. Franke, A. Jost, and T. Kleinjung, "Conditional Calculation of pi(10^24)", Posting to the Number Theory Mailing List, Jul 29 2010.  A little more detail (that's all I have right now) is contained in the talk of David Platt: www.maths.bris.ac.uk/~madjp/junior%20talk.pdf .   The idea in the combinatorial method is the following:
Suppose that $f(n)$ is a completely multiplicative function of the positive integers.  In the case of calculating $\pi(x)$, we take $f(n) = 1$.  In the case of calculating the sum of the primes $\le x$ we take $f(n) = n$.  In the case of Calculating $\pi(x,a,q)$ we take $f(n)$ to be a Dirichlet character mod $q$.
Define $\phi(x,a)$ as the sum of $f(n)$ for all integers $n \le x$ which are not divisible by the first $a$ primes.  Clearly $\phi(x,0) = \sum_{n \le x} f(n)$.  We have the recursion:
$\phi(x,a+1) = \phi(x,a) - f(p_{a+1})\phi(x/p_{a+1},a).$
Imagine that you have a labeled tree, where the labels are $\pm f(k) \phi(x/k,b)$ for some $k$ and $b$.  You can expand any node in the tree you want by applying the above formula.  This has the property that it leaves the sum of the values at the leaves the same.  So the idea is to start with a tree consisting of one node labeled $\phi(x,\pi(\sqrt x))$ and keep expanding until either $b=0$ or $x/k < $ some cutoff.  If you accumulate all such nodes you can evaluate all of them by sieving (and a special data structure -- see the Lagarias-Miller-Odlyzko paper for details), and the optimal value for the cutoff is $x^{2/3}$ perhaps multiplied by some logarithmic factors.
Here's the reference to my paper with Lagarias and Odlyzko: Computing $\pi(x)$: The Meissel-Lehmer Method
Author(s): J. C. Lagarias, V. S. Miller, A. M. Odlyzko
Source: Mathematics of Computation, Vol. 44, No. 170 (Apr., 1985), pp. 537-560
