kimwalisch
Riemann R approximation is a little off
The output of --Ri gives different numerical results than Wolfram Alpha and python's mpmath. For example,
For powers of 10 starting at 10^0 primecount --Ri gives:
0 4 25 168 1226 9587 78527 664667 5761551 50847455 455050683 4118052494 37607910542 346065531065 3204941731601 29844570495886 279238341360977 2623557157055978 24739954284239494 234057667300228940 2220819602556027016 21127269485932299722 201467286689188773984 1925320391607837270272 18435599767347541916672 176846309399141934432256 1699246750872419992338432 16352460426841662907482112 157589269275973235320553472 1520698109714271829697232896
Using mpmath riemannr I get:
from mpmath import riemannr, mp mp.dps = 40 # set high deliberately
for i in range(0, 30): print(riemannr(10**i))
1.0 4.564583141005090239865774695584049809675 25.66163326692418259322679794035569814997 168.3594462811673480649133109867329108466 1226.93121834343310855421625817211837034 9587.431738841973414351612923908229431098 78527.39942912770485887029214095925103488 664667.4475647477679853466998874388326848 5761551.86732016956230886495973466773615 50847455.42772142751394887577256049489582 455050683.3068469244631532415819991388608 4118052494.63140044176104610770875738881 37607910542.22591023474569601742946144013 346065531065.8260271978929257301899631105 3204941731601.689034750500754116280826964 29844570495886.92737822228672779202881061 279238341360977.1872302539272988649936199 2623557157055978.003875460015662127309679 24739954284239494.40252165144480328003146 234057667300228940.2346566885611062528803 2220819602556027015.401217592243516061431 21127269485932299723.73386404400837512661 201467286689188773625.1590118748480817499 1925320391607837268776.080252873251190086 18435599767347541878146.80335901928241809 176846309399141934626965.8309690195779633 1699246750872419991992147.2218584438395 16352460426841662910939464.57821529157634 157589269275973235652219770.5690766206166 1520698109714271830281953370.160723380532
Ignoring the few first values, these differ in the following ways,
a) The values are much the same until 2220819602556027015.4 which ends with a 6 in the primecount version. b) From there the values start to drift apart. For example the final value 1520698109714271830281953370.16 is 584720474.16 far from the primecount version.
I checked wolfram alpha and it gives results very close to mpmath and hence quite far from primecount.
Yes, I know about this issue. If you are on Linux and you use the GCC compiler, then you can build primecount using 128-bit floating point numbers which increases the Riemann R precision:
git clone https://github.com/kimwalisch/primecount.git
cd primecount
cmake . -DWITH_FLOAT128=ON
make -j8
The Riemann R precision is accurate enough for primecount's use case, i.e. we use to compute the nth prime. Because of this issue the Riemann R function is also not part of primecount's public C/C++ API. It is more of an internal helper function...
That’s good to know. The option is currently listed in the readme at https://github.com/kimwalisch/primecount and also when you use —help. Maybe adding a comment there would make sense.
What stops cmake from choosing WITH_FLOAT128=ON where possible?
Maybe adding a comment there would make sense.
I have done that: https://github.com/kimwalisch/primecount/commit/fdecedce6dfc58ecc353df8f015460c65ec8892b
What stops cmake from choosing WITH_FLOAT128=ON where possible?
I did some research and found a GCC bug that I reported some time ago: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=92826. The problem is that there is a GCC warning when using libquadmath with the -pedantic flag and there is no clean way to suppress this warning. I did implement a runtime CMake check to automatically enable float128 on the cmake_check_float128_support branch and it seems to work, but because of the GCC warning issue I am unsure whether it is worth merging this. If we merge this then there will be compiler warnings when compiling primecount with GCC & -pedantic.