Riemann R approximation is a little off

[lesshaste]

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.

[kimwalisch]

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...

[lesshaste]

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?

[kimwalisch]

Maybe adding a comment there would make sense.

I have done that: fdecedc

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.