Wow, I would be very interested in the code if you can paste it here. What is the largest prime number sage can handle on whatever machine you are doing this on? Are you saying that sage can work with Prime[Prime[10^14]]? And what about accuracy and how many digits of computation are you carrying?

I only work with primes up to a certain cutoff. Beyond that cutoff, I simply use the upper bound for $p_n$.

Okay, that's what I figured. That's the problem with the upper bound, it isn't that tight so your result of 1.004 could be too big of an overestimate and the actual sum might still be less than 1. I tried doing this in mathematica carrying a 100 digits with a hundred million terms using primes themselves. It took its sweet time (another reason to be surprised when you said it ran in about a minute) and the result was less than one. Don't quote me but I think it was something like 0.98 when I was just playing around with some of this stuff without any serious effort.

+1 nice work though, even computationally it is quite a challenge.

No, you're misunderstanding. 1.004 is a lower bound for the sum. Because we have an upper bound for $p_n$, we have an upper bound for $p_{p_n}$, so we have a lower bound for $1/p_{p_n}$, so we have a lower bound for the whole sum. It is completely rigorous.

By the way, I had previously computed an upper bound of 1.06, and with a longer computation I found a better lower bound of 1.02. This sum isn't that hard to work with, and I'm sure others have better results as well.

Okay so now as a chat, I'll have a chance to explain what I was doing and why it is not as naive as your approach...

Basically computing the sum term by term is very slow, because this series has not been "optimized" to converge quickly. Essentially, the term $p_{p_n}$ has very slow growth. So I exploit this by "chunking" up the sum.

I use chunks of the form (n,n+log(n)) because on intervals of length O(log n), the terms remain "constant" enough for our purposes.

So I replace sum_n^{n+log n}1/p_{p_n} with a chunked lower bound: log n/p_{p_{n+log n}}

Here, we have used the fact that p_{p_n} is monotonic increasing.

Okay, I see, nice work. I just played a bit with it last night trying to do it term by term and I could only go up a hundred million terms or so and it was still a bit less than one so left it at that.

Yes, so I am going way beyond a hundred million terms, but I am not sampling every term.

Basically, I had to go out to about 10^14 (a factor of 10^6 more terms than you).

The mathematica implementation on my machine only goes up to some certain number but my simple program either would break down before that or would take way too long and it was pretty late so I just left it at that.

So I am using arbitrary precision computation (enabled by default in Sage). I am sure there is a way to do this in Mathematica as well.

Yeah, I was carrying a hundred digits.

So after chunking, I think I had about 500 chunks total. So I only summed 500 things, in order to get a lower bound for summing 10^14 things. Does that make sense? That's why it ran so quickly.

In fact, with a little more math you could turn these "computational" tricks into a new series that converges rapidly to a lower bound of the sum in question. These are the kinds of tricks people use for approximating pi and stuff like that.

Sure, the idea was to see if goes above one or not. So if using the lower bounds takes the sum above one, then I guess the question is answered. Can you repeat this with more refined chunks, with more chunks but the width of each chunk being smaller?

Yes, but it takes longer :)

In fact, if you refine the chunks all the way down to 1, then you're just summing the series as usual.

How does the complexity/computational time increase as a function of number of chunks? Like how much longer would it be to double the number of chunks

Roughly linear in the number of chunks.

Cool, I would just try it refine until it hits a reasonable time like 10-15 mins and see what the lower bound comes to. I am curious if it will drop below one. Because your estimate is just so close to one. That would be funny if quadrupling the number of chunks drops the estimate to below one.

No, since it is a lower bound it will not drop, it will only increase.

Think of it this way: if S is the sum in question, I have proven that S > 1.004. So if T is a new lower bound for S that refines the old lower bound, T will be closer to S than 1.004, so S > T > 1.004. Does that make sense?

I'm not sure if you saw my comment, but I let the program run with more chunks and I got a lower bound of 1.02

Which was a refinement of the lower bound of 1.004

I saw that comment. I don't know why I keep confusing it with an upper bound. Well looks like you have nailed at least that part of the question then.

In that case, refining the chunks then will take you closer to the actual value of the sum

Any idea what that is? Like how large is it?

Anyway it's really late for me. I wanted to post the code but couldn't figure out how to publish a sage notebook on my website. So good night for now.

Take care, thanks!

Not sure about the exactly value. Probably near 1.03 or 1.04, I don't know