Gerhard Niklasch on Wed, 1 Sep 1999 01:42:39 +0200 (MET DST) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
ellisoncurve bug & suggested patch |
Xavier Roblot recently pointed out the following weirdness to me: > BTW, someone working on elliptic curves sent me the following > bug report (which is still present in the new version): > > ? E40 = [0,0,0,-7,-6]; > ? pol = x^3-7*x-7; > // clearly, [a,1] is on the curve for any root a of pol; > // but PARI doesn't agree > ? r = polroots(pol); > ? ellisoncurve(E40,[r[1],1]) > %1 = 1 > ? ellisoncurve(E40,[r[2],1]) > %2 = 0 > ? ellisoncurve(E40,[r[3],1]) > %3 = 0 There were in fact two bugs. One was of the one-character-typo kind and had been introduced waaaaaaaaay back in 2.0.2.alpha. The other was the logic in the internal function oncurve(). It computes the lefthand side y^2 + a_1*x*y + a_3*y and the righthand side (cubic in x) of the equation, and then their difference. When the difference is zero, it returns true (1). When the difference is nonzero and both sides are exact, it returns false (0). So far, so good. When one of the sides is exact and the other isn't, as was the case above (a_1 is zero and x is the only inexact quantity in sight, so the lefthand side is exact), it used to return false, which accounts for the strange output above. With the patch below, results look more intuitive: (01:28) gp > ellisoncurve(E40,[r[1],1]) %4 = 1 (01:28) gp > ellisoncurve(E40,[r[2],1]) %5 = 1 (01:28) gp > ellisoncurve(E40,[r[3],1]) %6 = 1 (01:29) gp > ellisoncurve(E40,[r[1],1.000000000000000000001]) %9 = 0 (01:29) gp > ellisoncurve(E40,[r[1],1.0000000000000000000000001]) %10 = 0 (01:29) gp > ellisoncurve(E40,[r[1],1.00000000000000000000000000001]) %11 = 1 (01:29) gp > ellisoncurve(E40,[r[3]+0.000000000000000001,1]) %12 = 0 (01:29) gp > ellisoncurve(E40,[r[3]+0.00000000000000000000000001,1]) %13 = 0 (01:29) gp > ellisoncurve(E40,[r[3]+0.00000000000000000000000000000001,1]) %14 = 1 Diffs are against the latest pre-2.0.17 snapshot, but it should go in with mild protest if you have 2.0.16.beta (line numbers are off by 2 in that case). Apply by hand against older versions; this particular corner of code hasn't changed during the last several beta versions. Please do let me know if you find this breaking anything else. Enjoy, Gerhard $ diff -c src/modules/elliptic.c.orig src/modules/elliptic.c *** src/modules/elliptic.c.orig Wed Sep 1 01:03:15 1999 --- src/modules/elliptic.c Wed Sep 1 01:24:14 1999 *************** *** 449,455 **** static long ellexpo(GEN e) { ! long i,k2, k = -gexpo((GEN)e[1]); for (i=2; i<6; i++) { k2 = gexpo((GEN)e[i]); --- 449,455 ---- static long ellexpo(GEN e) { ! long i,k2, k = gexpo((GEN)e[1]); for (i=2; i<6; i++) { k2 = gexpo((GEN)e[i]); *************** *** 458,463 **** --- 458,468 ---- return k; } + /* Exactness of lhs and rhs in the following depends in non-obvious ways + on the coeffs of the curve as well as on the components of the point z. + Thus if e is exact, with a1==0, and z has exact y coordinate only, the + lhs will be exact but the rhs won't.-- Changed to get rid of the resulting + weirdness; hope this doesn't break anything p-adically. --GN1999Sep01 */ int oncurve(GEN e, GEN z) { *************** *** 469,476 **** p2 = ellRHS(e,(GEN)z[1]); x = gsub(p1,p2); if (gcmp0(x)) { avma=av; return 1; } p = precision(p1); ! q = precision(p2); if (p < q) q = p; ! if (!q) { avma=av; return 0; } /* one of p1, p2 is exact */ /* the constant 0.93 is arbitrary */ q = (gexpo(x)-ellexpo(e) < - bit_accuracy(q) * 0.93); avma = av; return q; --- 474,482 ---- p2 = ellRHS(e,(GEN)z[1]); x = gsub(p1,p2); if (gcmp0(x)) { avma=av; return 1; } p = precision(p1); ! q = precision(p2); ! if (!p && !q) { avma=av; return 0; } /* both of p1, p2 are exact */ ! if (!q || (p && p < q)) q = p; /* min among nonzero elts of {p,q} */ /* the constant 0.93 is arbitrary */ q = (gexpo(x)-ellexpo(e) < - bit_accuracy(q) * 0.93); avma = av; return q;