I want to calculate the product over 1-1/p , where p runs over the primes upto 10^10
I know the approximation exp(-gamma)/ln(10^10) , where gamma is the Euler-Mascheroni-constant and ln the natural logarithm, but I want to calculate the exact product to see how close the approximation is.
The problem is that PARI/GP takes very long to calculate the prime numbers from about 4.2 * 10^9 to 10^10. The prodeuler-command also takes very long.
Is there any method to speed up the calculation with PARI/GP ?
I'm inclined to think the performance issue has mostly to do with the rational numbers rather than the generation of primes up to 10^10.
As a quick test I ran
a(n)=my(t=0);forprime(p=1,n,t+=p);t
with a(10^10) and it computed in a couple of minutes which seems reasonable.
The corresponding program for your request is:
a(n)=my(t=1);forprime(p=1,n,t*=(1-1/p));t
and this runs much slower than the first program, so my question would be to ask if there is a way to reformulate the computation to avoid rationals until the end? Is my formulation above even as you intended? - the numbers are extremely large even for 10^6, so it is no wonder it takes a long time to compute and perhaps the issue has less to do with the numbers being rational but just their size.
One trick I have used to compute large products is to split the problem so that at each stage the numbers on the left and right of the multiplication are roughly the same size. For example to compute a large factorial, say 8! it is much more efficient to compute ((1*8)*(2*7))*((3*6)*(4*5)) rather than the obvious left to right approach.
The following is a quick attempt to do what you want using exact arithmetic. It takes approximately 8mins up to 10^8, but the size of the numerator is already 1.9 million digits so it is unlikely this could ever get to 10^10 before running out of memory. [even for this computation i needed to increase the stack size].
xvecprod(v)={if(#v<=1, if(#v,v[1],1), xvecprod(v[1..#v\2]) * xvecprod(v[#v\2+1..#v]))}
faster(n)={my(b=10^6);xvecprod(apply(i->xvecprod(
apply(p->1-1/p, select(isprime, [i*b+1..min((i+1)*b,n)]))), [0..n\b]))}
Using decimals will definitely speed things up. The following runs reasonably quickly for up to 10^8 with 1000 digits of precision.
xvecprod(v)={if(#v<=1, if(#v,v[1],1), xvecprod(v[1..#v\2]) * xvecprod(v[#v\2+1..#v]))}
fasterdec(n)={my(b=10^6);xvecprod(apply(i->xvecprod(
apply(p->1-1.0/p,select(isprime,[i*b+1..min((i+1)*b,n)]))),[0..n\b]))}
The fastest method using decimals is the simplest:
a(n)=my(t=1);forprime(p=1,n,t*=(1-1.0/p));t
With precision set to 100 decimal digits, this produces a(10^9) in 2 minutes and a(10^10) in 22 minutes.
10^9: 0.02709315486987096878842689330617424348105764850
10^10: 0.02438386113804076644782979967638833694491163817
When working with decimals, the trick of splitting the multiplications does not improve performance because the numbers always have the same number of digits. However, I have left the code, since there is a potential for better accuracy. (at least in theory.)
I am not sure I can give any good advice on the number of digits of precision required. (I'm more of a programmer type and tend to work with whole numbers). However, my understanding is that there is a possibility of losing 1 binary digit of precision with every multiplication, although since rounding can go either way on average it won't be quite so bad. Given that this is a product of over 450 million terms, that would imply all precision is lost.
However, using the algorithm that splits the computation, each value only goes through around 30 multiplications, so that should only result in a loss of at most 30 binary digits (10 decimal digits) of precision so working with 100 digits of precision should be sufficient. Surprisingly, I get the same answers either way, so the simple naive method seems to work.
During my tests, I have noticed that using forprime is much faster than using isprime. (For example, the fasterdec version took almost 2 hours compared with the simple version which took 22 minutes to get to the same result.). Similary, sum(p=1,10^9,isprime(p)) takes approximately 8 minutes, compared with my(t=1);forprime(p=1,10^9,t++);t which takes just 11 seconds.