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

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

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.

Cheers

Jerome

-----Message d'origine-----
De : Karim Belabas [mailto:Karim.Belabas@math.u-bordeaux.fr] 
Envoyé : mercredi 10 janvier 2018 23:12
À : 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-10 21:01]:
[...]
> (20:23) gp > a=10^1000000; \\ 1 million decimal digits !
> time = 27 ms.
> (20:24) gp > sqrtint(a);
> time = 39 ms.
> (20:24) gp > sqrtnint(a,12);
> time = 208 ms.
> (20:24) gp > sqrtnint(a,123);
> time = 127 ms.
> (20:24) gp > sqrtnint(a,1234);
> time = 172 ms.
> (20:24) gp > sqrtnint(a,12345);
[...]
>   *** sqrtnint: Warning: increasing stack size to 640000000.
> time = 10,236 ms.
> (20:24) gp > sqrtnint(a,123456);
[...]
>   ***   at top-level: sqrtnint(a,123456)
>   ***                 ^------------------
>   *** sqrtnint: the PARI stack overflows !
>   current stack size: 1000001536 (953.676 Mbytes)
>   [hint] you can increase 'parisizemax' using default()

Indeed, the initial guess was far off in this case. I committed a fix to
master by making it ceil( exp(log(a)/n) ) [ computed with 64 bits of
accuracy ]:

(23:07) gp > a=10^1000000;
time = 12 ms.
(23:07) gp > sqrtint(a);
time = 24 ms.
(23:07) gp > sqrtnint(a,12);
time = 156 ms.
(23:07) gp > sqrtnint(a,123);
time = 88 ms.
(23:07) gp > sqrtnint(a,1234);
time = 72 ms.
(23:07) gp > sqrtnint(a,12345);
time = 56 ms.
(23:07) gp > sqrtnint(a,123456);
time = 24 ms.

Thanks for the hint and the detailed analysis !

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