 mersenneforum.org (https://www.mersenneforum.org/index.php)
-   y-cruncher (https://www.mersenneforum.org/forumdisplay.php?f=159)
-   -   Prime representing constant (https://www.mersenneforum.org/showthread.php?t=26298)

 Viliam Furik 2020-12-11 00:38

Prime representing constant

Inspired by [URL="https://youtu.be/_gCKX6VMvmU"]recent Numberphile video[/URL] (very short, 12 minutes long), I have started programming a Python program that could compute the constant to as many places as possible, in at most 24 hours, which is what I consider a reasonable time frame for the development phase. So far, the longest completed run computed 3 million decimal digits and took roughly 6 hours to complete. The collected time data suggests that if runtime for n digits is x, then runtime for an is xa[SUP]2[/SUP]. It also seems to be lower than that, by a factor of about 1.1, on average.

My program uses a slight variation of the algorithm used in the video (the infinite sum) -> As we know, what is the highest prime we will use, we can end the sum there, and sum up the fractions. Then we can factor out all the way to the beginning, and we get something like this:

( ( ( (2-1)*2 + (3-1) )*3 + (5-1) )*5 + (7-1) )*7 + (11-1)... / 2*3*5*7...

If you look at the bracket part, you will notice, that for each prime, you add it to the previous result, subtract 1, and then you multiply the whole thing by the prime, and then do the same for the next prime. At each step, you also multiply the denominator by the prime, continuing the primorial. When you are done, you simply divide the fraction into a number with a given number of digits.

That should be a lot faster than the original algorithm because you divide only in the end, thus keeping the numbers you work with as integers, until the end. Plus, it opens the possibility of backup files so that when you end, you can continue the computation later.

I have also found, that for n decimal digits, the highest prime I need is at least (1+x)*n/log[SUB]10[/SUB](e), where x tends to zero as n grows – for n = 25000, x = 0.02 should be sufficient. That results from Chebyshev’s function and some quick mathematics using logarithms.

A friend of mine is going to help me convert it into a C program, using hopefully a simple FFT multiplication (current Python program uses Karatsuba multiplication for simplicity). Until then, I will try if using NumPy speeds it up (it should, hopefully), and maybe I will get it to run on CUDA, through the Numba library.

My question to Alex Yee is, whether it is worth including in the y-cruncher. Another question to anyone, apart from FFT and C, are there any other improvements that can be done?

BTW, my program reads the primes from a sieved file containing primes under 1,000,000,000, enough to compute 400,000,000 digits. I found it is faster to read them from a file, rather than sieving them on the start of each run, i.e. the total time it takes it to read them is lower than the time it takes it to sieve them.

 retina 2020-12-11 00:52

[QUOTE=Viliam Furik;565898]Another question to anyone, apart from FFT and C, are there any other improvements that can be done?[/QUOTE]An obvious way to improve the speed is to use assembly.

But what is your goal here? Just to make it really fast and nothing else? If all the effort is just for a single run then it seems not worthwhile to go to a lot of trouble when you can let the computer do the hard work, albeit a little less efficiently and for longer, instead of you using your time while the computer sits idle.

 Dr Sardonicus 2020-12-11 02:01

See also the arxiv submission [url=https://arxiv.org/pdf/2010.15882.pdf]A Prime-Representing Constant[/url]. This includes a proof that the constant is irrational.

 Viliam Furik 2020-12-11 02:05

[QUOTE=Dr Sardonicus;565907]See also the arxiv submission [url=https://arxiv.org/pdf/2010.15882.pdf]A Prime-Representing Constant[/url]. This includes a proof that the constant is irrational.[/QUOTE]

But the same paper includes an algorithm to calculate the constant to an arbitrary precision (bottom of page 1).

 casevh 2020-12-11 04:24

Have you tried using the gmpy2 library? It provides access to the GMP library for multiple-precision arithmetic. It should be very easy to convert your existing program to use gmpy2.

I maintain gmpy2. If you have any questions, let's start a new thread under Programming.

casevh

 Viliam Furik 2020-12-11 11:17

[QUOTE=casevh;565915]Have you tried using the gmpy2 library? It provides access to the GMP library for multiple-precision arithmetic. It should be very easy to convert your existing program to use gmpy2.

I maintain gmpy2. If you have any questions, let's start a new thread under Programming.

casevh[/QUOTE]

I tried, but I couldn't install it. It errored out when installing. I will send you details in PM.

 Stargate38 2020-12-11 17:24

I hope you release this to the public when it's finished. I would love to try it out.

 Viliam Furik 2020-12-11 18:29

[QUOTE=Stargate38;565946]I hope you release this to the public when it's finished. I would love to try it out.[/QUOTE]

Sure, why not. It will be most probably open-source. I am now tweaking it with the help of casevh and his ultra-fast library.

 Dr Sardonicus 2020-12-11 20:15

[QUOTE=Viliam Furik;565909]But the same paper includes an algorithm to calculate the constant to an arbitrary precision (bottom of page 1).[/QUOTE]From the definitions of g[sub]n[/sub] and f[sub]1[/sub] as the limiting value of the g[sub]n[/sub], and using Bertrand's Postulate os in the paper to estimate the "tail" of the infinite series for the limiting value f[sub]1[/sub] (the constant 2.92005...), we obtain for n > 1

g[sub]n[/sub] < f[sub]1[/sub] < g[sub]n[/sub] + 2/p[sub]n-1[/sub]#.

Using the Prime Number Theorem, this estimate indicates that, for a given positive integer K,

0 < f[sub]1[/sub] - g[sub]n[/sub] < 10[sup]-K[/sup] if p[sub]n[/sub] > K*ln(10), approximately.

If this estimate is right, then for 3 million decimal digits' accuracy you need to take primes p up to about 6.9 million.

Of course, being able to maintain accuracy to K decimal digits has its own price if K is large.

 Viliam Furik 2020-12-11 21:09

[QUOTE=Dr Sardonicus;565961]Using the Prime Number Theorem, this estimate indicates that, for a given positive integer K,

0 < f[sub]1[/sub] - g[sub]n[/sub] < 10[sup]-K[/sup] if p[sub]n[/sub] > K*ln(10), approximately.

If this estimate is right, then for 3 million decimal digits' accuracy you need to take primes p up to about 6.9 million.[/QUOTE]
[QUOTE=Viliam Furik;565898]
I have also found, that for n decimal digits, the highest prime I need is at least (1+x)*n/log[SUB]10[/SUB](e), where x tends to zero as n grows – for n = 25000, x = 0.02 should be sufficient. That results from Chebyshev’s function and some quick mathematics using logarithms.[/QUOTE]
Yes, I know. It is right, because the ratio of x to theta(x) - first Chebyshev's function, tends to 1, as x goes to infinity. As the theta(x) is asymptotically equal to ln(p#), and that ln(x#) = log10(x#) / log10(e), the log10(x#) (about the number of digits we want to calculate) is equal to ln(x#) * log10(e), which is equal to x * log10(e). And because the value of the function x/theta(x) is slightly higher than 1, even when x = 100,000, then we can multiply the expression by 1+x, where x is some small number.

[QUOTE=Dr Sardonicus;565961]Of course, being able to maintain accuracy to K decimal digits has its own price if K is large.[/QUOTE]
[QUOTE=Viliam Furik;565898]If you look at the bracket part, you will notice, that for each prime, you add it to the previous result, subtract 1, and then you multiply the whole thing by the prime, and then do the same for the next prime. At each step, you also multiply the denominator by the prime, continuing the primorial. When you are done, you simply divide the fraction into a number with a given number of digits.[/QUOTE]
That's why I am calculating the approximation by computing numerator and denominator separately. I divide only after all the computation has been done.

 Viliam Furik 2020-12-15 01:02

[QUOTE=Stargate38;565946]I hope you release this to the public when it's finished. I would love to try it out.[/QUOTE]

I have some great development news:

Program is fast enough to compute 100,000,000 digits in less than 2 hours, and 1,000,000 in one second. Unfortunately, the runtime scaling of the generating algorithm is quadratic (imperfectly), because it iterates through all primes less than a maximum bound, which I set to be 2.5 times the number of digits.

All times are measured for a single-threaded run because the multi-core parallelization hasn't been successfully implemented yet.

I am hugely thankful to casevh not only for his gmpy2 library but also for his help with programming. Without his contribution, I wouldn't have got this far.

The first publicly available version will be made when the program reaches 1 billion (10^9) decimal places under a day, which, I hope, will be soon.

All times are UTC. The time now is 04:39.