Re: sqrtnint algorithm improvement for huge arguments (plus memory leak issue)

[Karim Belabas <Karim.Belabas@math.u-bordeaux.fr>]

* Jérôme Raulin [2018-01-11 08:29]:
> Hi Karim and thanks for your prompt answer.
> 
> Indeed performance improvement is impressive, my solution to latest project
> Euler problem got a 10x performance increase.
> Still there must be a side effect in your fix since sqrtnint givers wrong
> answer in some specific cases :
>   Latest development version including your fix :
>   gp > sqrtnint(10^500,27)
>   %1 = 3300034791125284633
>   Previous version :
>   gp > sqrtnint(10^500,27)
>   %1 = 3300034791125284639
> The second being the correct answer.
> Other cases are : (10^200, 11), (10^300, 8), (10^400, 11), (10^400, 22)
> In some cases the error is very small but for example for sqrtnint(10^300,
> 8) the error is huge :
>   gp > sqrtnint(10^300,8)
>   %1 = 31622776601683793296491743452850028544
>   Previous version :
>   gp > sqrtnint(10^300,8)
>   %1 = 31622776601683793319988935444327185337
> 
> I add a look at the new code and it may be due to the fact that the bailout
> criteria makes the hypothesis that the original approximation is always
> greater or equal than the result and that now ceil(exp(log(a)/n)) can
> sometimes give approximation lower than result. Still you seems to take care
> of this case while copying x to xs.

Or so I thought. Indeed, there is no guarantee that exp(log(a)/n) is
correct to 1 ulp: atomic operations are, not their composition, and the
difference is bounded by log(2^BIT_INLONG) ~ 44.36 in our case. So my + 1
was definitely not enough.

I just increased the accuracy for this initial floating point computation
to 128 bits if the result is > 2^32, before truncating to the requested
64 bits rounding up, and the above problems disappear.

Thanks for checking the patch !

Cheers,

    K.B.
--
Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
Universite de Bordeaux         Fax: (+33) (0)5 40 00 21 23
351, cours de la Liberation    http://www.math.u-bordeaux.fr/~kbelabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux.fr/  [PARI/GP]
`