Re: usqrtsafe isn't

[Bill Allombert <Bill.Allombert@math.u-bordeaux1.fr>]

On Mon, Nov 16, 2009 at 12:28:19AM -0500, Charles Greathouse wrote:
> In writing a program using the Pari library today I was able to trace
> a problem to usqrtsafe.  In particular,
> usqrtsafe(18446744073709551615) = 4294967296
> but 4294967296^2 = 18446744073709551616 > 18446744073709551615.

Thanks a lot for reporting this problem!

The code suppose the error is always by default by actually it is always
by excess, because the conversion to 'double' round to the nearest and
usqrtsafe should round to the lowest.

The patch below should fix that.

Cheers,
Bill.


diff --git a/src/kernel/none/mp_indep.c b/src/kernel/none/mp_indep.c
index ef5f98b..e9691e6 100644
--- a/src/kernel/none/mp_indep.c
+++ b/src/kernel/none/mp_indep.c
@@ -700,7 +700,7 @@ usqrtsafe(ulong a)
 {
   ulong x = (ulong)sqrt((double)a);
 #ifdef LONG_IS_64BIT
-  ulong y = x+1; if (y <= LOWMASK && y*y <= a) x = y;
+  if (x > LOWMASK || x*x > a) x--;
 #endif
   return x;
 }