eclib
eclib copied to clipboard
JohnCremona
32 bit test failure on h1bsdcurisog and h1bsd
Testing h1bsdcurisog...h1bsdcurisog completed
4c4
< 37 a 1 [0,0,1,-1,0] 1 1 1 5.9869172924639192596640199589 0.30599977383405230182048368335 0.051111408239968840235886099757 1.0000000000000000000000000001
---
> 37 a 1 [0,0,1,-1,0] 1 1 1 5.9869172924639192596640199589 0.30599977383405230182048368337 0.051111408239968840235886099757 1.0000000000000000000000000001
31c31
< 389 a 1 [0,1,1,-2,0] 2 1 1 4.9804251217101101506427155839 0.75931650028842677023019260781 0.15246017794314375162432475705 0.99999999999999999999999999989
---
> 389 a 1 [0,1,1,-2,0] 2 1 1 4.9804251217101101506427155839 0.75931650028842677023019260774 0.15246017794314375162432475705 0.9999999999999999999999999998
and
Testing h1bsd...h1bsd completed
6c6
< 53 a Rank = 1 L^(r)(f,1)/r! = 0.43586382417785716 L(f,1)/period = 0
---
> 53 a Rank = 1 L^(r)(f,1)/r! = 0.43586382417785719 L(f,1)/period = 0
15c15
< 58 a Rank = 1 L^(r)(f,1)/r! = 0.46370416478825057 L(f,1)/period = 0
---
> 58 a Rank = 1 L^(r)(f,1)/r! = 0.46370416478825061 L(f,1)/period = 0
The behaviour is not new to 20240408. The test failure shows now because the precision was changed from 30 to 100 bits.
I traced this to a difference in myg1(..) coming from Euler(), as shown by the following program
#include <eclib/periods.h>
int main(void)
{
set_precision(100);
bigfloat g1 = myg1(TWOPI/sqrt(to_bigfloat(37)));
cout << "myg1(2π/√37) = " << g1 << endl;
cout << " = 0.2076512696155670653383817926 + "
<< g1 - to_RR("0.2076512696155670653383817926") << endl;
bigfloat gamma = Euler();
cout << "Euler() = " << gamma << endl;
cout << " = 0.57721566490153286060651209 + "
<< gamma - to_RR("0.57721566490153286060651209") << endl;
}
On 64 bit, the output is
$ g++ myg1.cc -o myg1 -lec -lntl && ./myg1
myg1(2π/√37) = 0.20765126961556706533838179268
= 0.2076512696155670653383817926 + 0.83027610274511492519584431496e-28
Euler() = 0.57721566490153286060651209006
= 0.57721566490153286060651209 + 0.63108872417680944432938285223e-28
but on 32 bit the output is
$ g++ myg1.cc -o myg1 -lec -lntl && ./myg1
myg1(2π/√37) = 0.20765126961556706533838179268
= 0.2076512696155670653383817926 + 0.75139001222301374465467145843e-28
Euler() = 0.57721566490153286060651209007
= 0.57721566490153286060651209 + 0.70997481469891062487055570875e-28
A workaround to get the same output on 32 bits is:
--- a/libsrc/interface.cc
+++ b/libsrc/interface.cc
@@ -103,7 +103,7 @@ void Compute_Euler(RR& y)
l = RR::precision();
- x = 1 + static_cast<long>((0.25 * (l - 3)) * (NTL_BITS_PER_LONG * LOG2));
+ x = 1 + static_cast<long>((0.25 * (l - 3)) * (64 * LOG2));
n = 1 + static_cast<long>(3.591 * x);
a=x;
However, I wonder if there's something amiss here. Computing to 120 bits instead gives
Euler() = 0.57721566490153286060651209008240243
= 0.57721566490153286060651209 + 0.82402431597388244500026313156100052e-28
The difference between the two values is on the order of 1e-29 which is 95-96 bits of precision instead of 100 (I think one would expect ...900082 for 100 bits).
I cannot test on a 32-bit machine, and would not want to if I could, so why would I not just say somewhere in the documentation that 32-bits is not supported? Otherwise I would have to work out how to do the tests you are doing and double up all the make check code. It is already bad enough that the option --disable-mpfp results are differenc using clang than gcc. I do not want to spend the rest of my life thinking about such things.
I cannot test on a 32-bit machine, and would not want to if I could, so why would I not just say somewhere in the documentation that 32-bits is not supported? Otherwise I would have to work out how to do the tests you are doing and double up all the make check code. It is already bad enough that the option --disable-mpfp results are differenc using clang than gcc. I do not want to spend the rest of my life thinking about such things.
That is fine and wise. I don't use 32-bit either. We do support 32 bit on void linux and we build and test our packages on i686 (which is really a x86_64 box running a 32-bit userspace chroot).
However, differences might be a symptom of something else. In this case, Euler() doesn't give the correct 100 bit result, but I don't know if this is important.
If this is not important I'll make a PR with the workaround I found so the tests pass in i686.
Thanks Gonzalo @tornaria . To be honest I am not sure where the code for Euler() came from, possibly LiDiA from years ago. I have never used static_cast anywhere else, which is telling. The only place in the code where the gamma constant is used is in computing higher L-values, and one reason I am not so concerned to fix it is that when I compute data for the LMFDB now I used a Sage script (calling either Magma or pari) for these values anyway.
One day I will move to making much more use of FLINT in eclib, including for all its floating point stuff, and then all these issues will (I expect) go away). But it will be rather a lot of work.
@tornaria in case you are still intending to make a PR, please checkout the latest master branch first. I have made a lot of mostly trivial changes after running cppcheck systematically on all source files. But have not (yet) added a new tag or made a new release.
@JohnCremona if Euler() is not important, I could make a PR to switch NTL_BITS_PER_LONG to 64 so the result is the same on 32 bits, if that's ok with you.
Otherwise, maybe this can be switched to use pari mpeuler() if available? Or from flint/arb?
An alternative, is to run the code in Euler() with a bit of extended precision (seems 8-10 extra bits is ok, maybe add 20 to be safe?). This is not very hard and maybe preferable to the first workaround I mentioned, as it would be an improvement on 64 bits as well.
Something like this is what I have in mind:
--- a/libsrc/interface.cc
+++ b/libsrc/interface.cc
@@ -102,6 +102,7 @@ void Compute_Euler(RR& y)
bigfloat u, v, a, b, c;
l = RR::precision();
+ RR::SetPrecision(l+20);
x = 1 + static_cast<long>((0.25 * (l - 3)) * (NTL_BITS_PER_LONG * LOG2));
n = 1 + static_cast<long>(3.591 * x);
@@ -124,6 +125,7 @@ void Compute_Euler(RR& y)
add(u, u, a);
add(v, v, b);
}
+ RR::SetPrecision(l);
div(y, u, v);
}
@tornaria Thanks. I am cleaning up the code in periods.cc some more as well. myg0(x) is just exp(-x) and for r>0 (not just r>1) mygr(x) is redundant, begin handled better by G(r,x), implemented a la Cohen. I will add in your trick of locally increasing precision. If you look at the function Q() you'll see that I hard-wire a fixed largish precision for zeta(2), zeta(3), zeta(4) anyway, perhaps I should just do that for Euler().
The plan is to switch to flint/arb for all floating points stuff (and more), but not right away. We do have libpari available, but I would need a conversion function from a pari real to an NTL RR for that to work.
@tornaria if you can test the current commit d209a2d please would you do so and report here? I think it is better, pending more work to use a library which does this much better than I even will. Thanks!
@tornaria if you can test the current commit d209a2d please would you do so and report here? I think it is better, pending more work to use a library which does this much better than I even will. Thanks!
Sure, I will do later today.
If you apply this:
--- a/libsrc/interface.cc
+++ b/libsrc/interface.cc
@@ -125,8 +125,8 @@ void Compute_Euler(RR& y)
add(u, u, a);
add(v, v, b);
}
- div(y, u, v);
RR::SetPrecision(l);
+ div(y, u, v);
}
long prec() {return RR::precision();}
the last division is done with the original precision, and the result is then more stable.
BTW, something like this:
--- a/libsrc/interface.cc
+++ b/libsrc/interface.cc
@@ -104,7 +104,7 @@ void Compute_Euler(RR& y)
l = RR::precision();
RR::SetPrecision(l+20);
- x = 1 + static_cast<long>((0.25 * (l - 3)) * (NTL_BITS_PER_LONG * LOG2));
+ x = 1 + static_cast<long>(0.25 * (l - 3));
n = 1 + static_cast<long>(3.591 * x);
a=x;
also works, and it's 30-40 times faster. Simple program to test:
#include <eclib/interface.h>
#include <eclib/timer.h>
#define MAX_PREC 1000
int main(void)
{
set_precision(2*MAX_PREC);
bigfloat gamma = Euler();
bigfloat log2 = log(to_bigfloat(2));
int count = 0;
start_time();
for(int i=0; i<=MAX_PREC; i++) {
set_precision(i);
bigfloat a = Euler() - gamma*1;
if (a != 0) {
bigfloat b = floor(log(abs(a)) / log2);
cout << "prec = " << i << ", bit error = " << b << endl;
if (++count >= 20) break;
}
}
stop_time();
cout << "Tested Euler() with precision 1 to " << MAX_PREC;
show_time();
cout << endl;
}
I've built and tested the void-linux package for d209a2d, everything works fine (including i686).
Thanks -- I put in both your suggestions. I may tag and make a new release soon.