Re: bug in 2-adic sqrt

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

On Sat, 24 Nov 2001, John Cremona wrote:
> This would seem to be a bug in 2-adic sqrt, in 2.1.1.  I am well aware of the
> fact that one loses one bit of 2-adic precision doing this, but we seem to be
> losing two bits.

Indeed, that's a bug which has been with us from the start (was already there
in 1.39.15). Thanks for spotting this!

> Example:
>
> (12:25) gp > a=-71+O(2^5)
> %1 = 1 + 2^3 + 2^4 + O(2^5)
> (12:25) gp > b=sqrt(a)
> %2 = 1 + 2^2 + 2^3 + O(2^4)
> (12:25) gp > b^2-a
> %3 = 2^4 + O(2^5)
>
> Here a = -71 (mod 32) = 25 (mod 32), and b^2=a (mod 32)  has general solution 5
> (mod 16), so the "right" answer is surely
>
> (12:25) gp > c=5+O(2^4)
> %4 = 1 + 2^2 + O(2^4)

Actually, the one you now get [after a trivial fix: GP assumed that c = 1 mod 8
was a valid starting point in all cases; which should have been either
c = 1 mod 4, or c=1,3 mod 8] is

(17:36) gp > sqrt(25+O(2^5))
%1 = 1 + 2 + 2^3 + O(2^4)

You always get the square root which is smallest mod p (mod 8 if p = 2).

> [What I really would like is to be able to ask for sqrt(Mod(-71,32)) and get
> Mod(5,16); I am only using 2-adics to simulate this.]

The philosophy of INTMODs is to mimic the behaviour of a given prime field.
In particular, the domain is usually fixed ( Mod(0,2)^2 = Mod(0,2), not
Mod(0,4) ), unless you leave GP no other option ( Mod(1,2) * Mod(1,3) =
Mod(0,1) ). PADICs on the contrary have varying precisions.

Also PADICS have built-in valuation, whereas it would have to be computed for
INTMODs (e.g: ooops square root algorithm is failing, let's check... yes,
modulus is a prime power, OK, compute p and exponent, then call the padic
routine). And why support prime powers and not general integers ? So now we
may need to factor integers for each non-trivial modular functions etc, etc.
(which is of course very inefficient since it could be done once and for
all). So the current philosophy is: in case of possible trouble, assume the
modulus is prime.

Admittedly it's cumbersome (user's responsibility to work in the various
padic fields then apply chinese remaindering for instance), but INTMODs as
they currently stand are not the right way to do this. Also, it's never been
requested before :-)

Cheers,

    Karim.

P.S: I made a more extensive patch in the unstable branch, which makes the
thing about twice faster at low accuracies ( at high accuracies, everything
is dominated by the inverse in the last y <-- (y + x/y) / 2 standard
iteration ). The following patch is minimal and applies to version 2.1.1.

[ of course, and as always, the preferred way to apply patches is to use the
CVS server.  Does anybody still apply patches by hand from a saved mail
message ??? ]

*** src/basemath/trans1.c       2000/11/03 21:00:23     1.29
--- src/basemath/trans1.c       2001/11/25 16:48:30     1.29.2.1
***************
*** 675,684 ****
    pp = precp(x);
    if (egalii(gdeux, (GEN)x[2]))
    {
!     lp=3; y[4] = un; r = mod8((GEN)x[4]);
!     if ((r!=1 && pp>=2) && (r!=5 || pp!=2)) err(sqrter5);
!     if (pp <= lp) { setprecp(y,1); return y; }

      x = dummycopy(x); setvalp(x,0); setvalp(y,0);
      av2=avma; lim = stack_lim(av2,2);
      y[3] = lstoi(8);
--- 675,685 ----
    pp = precp(x);
    if (egalii(gdeux, (GEN)x[2]))
    {
!     lp=3; r = mod16((GEN)x[4]);
!     if (((r&7)!=1 && pp>=2) && ((r&3)!=1 || pp!=2)) err(sqrter5);
!     if (pp <= lp) { y[4] = un; setprecp(y,1); return y; }

+     y[4] = r==1? un: lstoi(3);
      x = dummycopy(x); setvalp(x,0); setvalp(y,0);
      av2=avma; lim = stack_lim(av2,2);
      y[3] = lstoi(8);

-- 
Karim Belabas                    Tel: (+33) (0)1 69 15 57 48
Dép. de Mathematiques, Bat. 425  Fax: (+33) (0)1 69 15 60 19
Université Paris-Sud             Email: Karim.Belabas@math.u-psud.fr
F-91405 Orsay (France)           http://www.math.u-psud.fr/~belabas
--
PARI/GP Home Page: http://www.parigp-home.de/