test failure in 32 bit due to bad rounding in set_precision

[tornaria]

The first test failure is this:

Testing comptest...
7c7
< AGM((1,1),(2,1)) = (1.47134936286465357,1.02005412633766002)
---
> AGM((1,1),(2,1)) = (1.471349362864653566,1.02005412633766002)
11c11
< Enter a real or complex: Main cube root = (1.62893714592217588,0.52017450230454584)
---
> Enter a real or complex: Main cube root = (1.628937145922175875,0.5201745023045458395)
15c15
< Next cube root = (-0.363984239564424216,-1.67078820068899634)
---
> Next cube root = (-0.3639842395644242155,-1.670788200688996339)
18,21c18,21
< (-1.28781547955764799,-0.857896758328490286)
< (-1.28781547955764799,0.857896758328490286)
< (0.287815479557647989,1.41609308017190794)
< (0.287815479557647989,-1.41609308017190794)
---
> (-1.287815479557647989,-0.8578967583284902864)
> (-1.287815479557647989,0.8578967583284902864)
> (0.2878154795576479889,1.416093080171907939)
> (0.2878154795576479889,-1.416093080171907939)

There are others but it seems they all arise from set_precision(70) which evaluates log(0.3*70) as 20 instead of 21. I suppose the issue is that 0.3 rounds up in double precision (which happens in 64bit) but it rounds down in single precision (which happens 32bit).

Workaround: replace 0.3 by 0.3L so it's always a double. That's enough to make all tests pass for me in 32bit.

Might as well change it to 0.30103L which is more accurate, 0.3010299956639812L, or M_LN2/M_LN10.

[JohnCremona]

Thanks for the report, Gonzalo. I'll make the suggested changes, though I do not have a 32-bit machine to test it on. I'll do it via a pull request to the master branch and would appreciate your testing it when I do.

[tornaria]

Sure. If it is of any help, the patch I'm using right now is below, with this all tests pass. I can also do a PR if you want.

This changes 0.3 by 0.3L, otherwise tests fail in i686 due to different rounding.

See: https://github.com/JohnCremona/eclib/issues/68

--- libsrc/eclib/interface.h    2021-03-18 12:22:52.000000000 -0300
+++ libsrc/eclib/interface.h    2021-04-05 15:15:52.355310137 -0300
@@ -240,7 +240,7 @@
 
 // Set internal precision to n bits and output precision to (0.3*n)-1 decimal places
 inline void set_precision(long n)
-{RR::SetPrecision(n); RR::SetOutputPrecision(long(0.3*n)-1);}
+{RR::SetPrecision(n); RR::SetOutputPrecision(long(0.3L*n)-1);}
 
 // Mostly for backward compatibility (used in saturate.cc) or for
 // temporarily changing internal precision when no output is relevant:
@@ -252,7 +252,7 @@
   {long n; cerr<<prompt<<": "; cin>>n; set_precision(n);}
 
 // read current precision converted to decimal (approximately)
-inline long decimal_precision() {return long(RR::precision()*0.3);}
+inline long decimal_precision() {return long(RR::precision()*0.3L);}
 
 // read current bit precision
 inline long bit_precision() {return RR::precision();}
@@ -289,7 +289,7 @@
 inline int is_approx_zero(double x) {return fabs(x)<1e-10;}
 
 // We cannot set internal bit precision in this mode, so we just set the output decimal precision
-inline void set_precision(long n) {cout.precision(min(15,long(0.3*n)));}
+inline void set_precision(long n) {cout.precision(min(15,long(0.3L*n)));}
 inline void set_precision(const string prompt)  {cout.precision(15);}
 #define Pi()    3.1415926535897932384626433832795028841
 #define Euler() (0.57721566490153286060651209008240243104)

[JohnCremona]

Please test the constants branch. Instead of putting 0.3L into the code in a few places I used a const double set to a higher precision value.

[tornaria]

Please test the constants branch. Instead of putting 0.3L into the code in a few places I used a const double set to a higher precision value.

All tests pass now, thanks!