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

[Jérôme Raulin <jerome@raulinfoissac.fr>]

Thanks a lot Karim, it seems to work correctly now and the performance are
excellent (in the 10^10000 range I got a 40x improvement on my script).
Cheers.
Jerome

-----Message d'origine-----
De : Karim Belabas [mailto:Karim.Belabas@math.u-bordeaux.fr] 
Envoyé : jeudi 11 janvier 2018 10:17
À : pari-dev@pari.math.u-bordeaux.fr
Objet : Re: sqrtnint algorithm improvement for huge arguments (plus memory
leak issue)

* 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]
`