Karim BELABAS on Sun, 25 Nov 2001 18:10:01 +0100 (MET) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: bug in 2-adic sqrt |
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/