[sage-devel] Re: elliptic curve division polynomials
2008/6/10 William Stein <[EMAIL PROTECTED]>: > > On Tue, Jun 10, 2008 at 9:59 AM, John Cremona <[EMAIL PROTECTED]> wrote: >> >> 2008/6/10 William Stein <[EMAIL PROTECTED]>: >>> >>> On Tue, Jun 10, 2008 at 7:34 AM, John Cremona <[EMAIL PROTECTED]> wrote: The following should be of interest to David Harvey, David Kohel, Chris Wuthrich and WIlliam Stein (at least!). Following from my reviewing of Chris's #3377 (torsion subgroups of elliptic curves over number fields) where I found a bug (see #3383) in William's new full_division_polynomial() code (see #3109): The bug was caused by some horrendous Grobner basis reduction in a 2-variable polynomial ring over a number field. >>> >>> Is this a bug in Grobner basis or in my code that happens to use it? >> >> The latter. You coerce things into a ring of the form K[X,Y]/I where >> I is an ideal, in order to "reduce" them modulo the defining equation >> of the curve (which is the generator of I). That leads to a reduce(0 >> function in >> quotient_ring_element.py being called on the elements with respect to >> that ideal, which in turn leads to code in t-y_buchberger.py where the >> Error is raised. Since I is a principal ideal it should not be too >> hard to find its G. basis? > > It seems to me that you just wrote "the bug is in my code", but then > when you described the bug it seems to be that it occurs in > Martin Albrecht's Toy Buchburger implementation? You are right, I should have said "the former". > >> >> >>> I was able to eliminate it by reworking William's full_division_polynomial() to work in the genuine function field of the curve. Thus, if K is the base_ring of the curve I construct in order the polynomial ring K[x], then its field of fractions K(x), then -- via the polynomial ring K(x)[Y] -- the field K(x)(y) where y^2=polynomial in x is the equation of the curve. (I am currently assuming that the curve is in short Weierstrass form as WIlliam did, though that wilkl be fixed in due course). Elements g of that field have two components g[0] and g[1], both rational functions in x, representing g[0]+y*g[1]. It is possible though slightly awkward to evaluate such a function g at a point P. What works is lift(g).subs(Y=P[1]).subs(x=P[0]) Now full_division_polynomial(m) returns an element g of that field, which if m is odd happens to be a polynomial in x (no denominator, no y term, i.e. g[1]==0) and if m is even happens to be y times a poly on x (no denominator and g[0]==0). Using that I reworked multiplication_by_m() as necessary and hence division_points(). doctests all pass, and the example which previously failed in #3383 now passes. Problem 1: Over Q, multiplication_by_m(10) takes ages -- over 3 minutes! -- even with the x_only=True option. The time is almost all spent in dividing one element of QQ(x) by another. That is just a simple division of two rational functions (polynomials of degrees 99 and 100 with no common factor) -- possibly the time is taken computing their gcd? I know that ZZ[x] operations are faster than QQ[x] but it would be very ugly to try to detect this special case since the code is otherwise generic and works over any field (char. not 2 or 3 until we stop using short W. equations). >>> >>> It seems like you need to investigate this a bit further. If it is really >>> gcd, then we could just improve QQ[x]'s gcd code (e.g., clear denoms >>> and call something in ZZ[x] for now; I don't know.) >>> >> >> It is definitely the gcd. I was able to avoid the slowdown by >> constructing the element of the fraction field L not by x=n/d (where n >> and d are the numer- and denominator), nor by x=L(n,d) but by >> x = >> sage.rings.fraction_field_element.FractionFieldElement(L,n,d,reduce=False) >> which forces the gcd not to be taken. I am pretty sure that at this >> point in the code we know that the polynomials are coprime anyway, so >> I would vote for the constructor of the form L(n,d) to have a flag >> reduce whcih defaults to True but can be set to False for situations >> like that. > > +1 > >> Problem/question 2: looking at the code in ell_generic.py I find that there are now *3* different implementations of division polys of various kinds. (1) E.division_pol = E.torsion_pol is by David Kohel (2005-04-25) and returns a polynomial in x alone, with a factor of the 2-division poly in x when m is even. (2) E.full_division_pol is William's new code (which can be asked to call the preceding when m is odd) which returns a polynomial in x and y, i.e. in the 2-variable poly ring over K, which happens to be of the form p(x) or y*p(x) where p is a poly in x alone (but belonging to K[x,y]). (3) There is code which is all commented out due to David Har
[sage-devel] Re: elliptic curve division polynomials
On Tue, Jun 10, 2008 at 9:59 AM, John Cremona <[EMAIL PROTECTED]> wrote: > > 2008/6/10 William Stein <[EMAIL PROTECTED]>: >> >> On Tue, Jun 10, 2008 at 7:34 AM, John Cremona <[EMAIL PROTECTED]> wrote: >>> >>> The following should be of interest to David Harvey, David Kohel, >>> Chris Wuthrich and WIlliam Stein (at least!). >>> >>> Following from my reviewing of Chris's #3377 (torsion subgroups of >>> elliptic curves over number fields) where I found a bug (see #3383) >>> in William's new full_division_polynomial() code (see #3109): >>> >>> The bug was caused by some horrendous Grobner basis reduction in a >>> 2-variable polynomial ring over a number field. >> >> Is this a bug in Grobner basis or in my code that happens to use it? > > The latter. You coerce things into a ring of the form K[X,Y]/I where > I is an ideal, in order to "reduce" them modulo the defining equation > of the curve (which is the generator of I). That leads to a reduce(0 > function in > quotient_ring_element.py being called on the elements with respect to > that ideal, which in turn leads to code in t-y_buchberger.py where the > Error is raised. Since I is a principal ideal it should not be too > hard to find its G. basis? It seems to me that you just wrote "the bug is in my code", but then when you described the bug it seems to be that it occurs in Martin Albrecht's Toy Buchburger implementation? > > >> >>> I was able to >>> eliminate it by reworking William's full_division_polynomial() to work >>> in the genuine function field of the curve. Thus, if K is the >>> base_ring of the curve I construct in order the polynomial ring K[x], >>> then its field of fractions K(x), then -- via the polynomial ring >>> K(x)[Y] -- the field K(x)(y) where y^2=polynomial in x is the equation >>> of the curve. (I am currently assuming that the curve is in short >>> Weierstrass form as WIlliam did, though that wilkl be fixed in due >>> course). Elements g of that field have two components g[0] and g[1], >>> both rational functions in x, representing g[0]+y*g[1]. >>> >>> It is possible though slightly awkward to evaluate such a function g >>> at a point P. What works is lift(g).subs(Y=P[1]).subs(x=P[0]) >>> >>> Now full_division_polynomial(m) returns an element g of that field, >>> which if m is odd happens to be a polynomial in x (no denominator, no >>> y term, i.e. g[1]==0) and if m is even happens to be y times a poly on >>> x (no denominator and g[0]==0). >>> >>> Using that I reworked multiplication_by_m() as necessary and hence >>> division_points(). doctests all pass, and the example which >>> previously failed in #3383 now passes. >>> >>> Problem 1: Over Q, multiplication_by_m(10) takes ages -- over 3 >>> minutes! -- even with the x_only=True option. The time is almost all >>> spent in dividing one element of QQ(x) by another. That is just a >>> simple division of two rational functions (polynomials of degrees 99 >>> and 100 with no common factor) -- possibly the time is taken >>> computing their gcd? I know that ZZ[x] operations are faster than >>> QQ[x] but it would be very ugly to try to detect this special case >>> since the code is otherwise generic and works over any field (char. >>> not 2 or 3 until we stop using short W. equations). >> >> It seems like you need to investigate this a bit further. If it is really >> gcd, then we could just improve QQ[x]'s gcd code (e.g., clear denoms >> and call something in ZZ[x] for now; I don't know.) >> > > It is definitely the gcd. I was able to avoid the slowdown by > constructing the element of the fraction field L not by x=n/d (where n > and d are the numer- and denominator), nor by x=L(n,d) but by > x = sage.rings.fraction_field_element.FractionFieldElement(L,n,d,reduce=False) > which forces the gcd not to be taken. I am pretty sure that at this > point in the code we know that the polynomials are coprime anyway, so > I would vote for the constructor of the form L(n,d) to have a flag > reduce whcih defaults to True but can be set to False for situations > like that. +1 > >>> >>> Problem/question 2: looking at the code in ell_generic.py I find that >>> there are now *3* different implementations of division polys of >>> various kinds. >>> >>> (1) E.division_pol = E.torsion_pol is by David Kohel (2005-04-25) and >>> returns a polynomial in x alone, with a factor of the 2-division poly >>> in x when m is even. >>> (2) E.full_division_pol is William's new code (which can be asked to >>> call the preceding when m is odd) which returns a polynomial in x and >>> y, i.e. in the 2-variable poly ring over K, which happens to be of the >>> form p(x) or y*p(x) where p is a poly in x alone (but belonging to >>> K[x,y]). >>> (3) There is code which is all commented out due to David Harvey >>> (2006-09-24) consisting of 3 functions pseudo_torsion_polynomial(), >>> multiple_x_numerator(() and multiple_x_denominator(). The first is >>> the same as before for odd m but for even m just omi
[sage-devel] Re: elliptic curve division polynomials
2008/6/10 William Stein <[EMAIL PROTECTED]>: > > On Tue, Jun 10, 2008 at 7:34 AM, John Cremona <[EMAIL PROTECTED]> wrote: >> >> The following should be of interest to David Harvey, David Kohel, >> Chris Wuthrich and WIlliam Stein (at least!). >> >> Following from my reviewing of Chris's #3377 (torsion subgroups of >> elliptic curves over number fields) where I found a bug (see #3383) >> in William's new full_division_polynomial() code (see #3109): >> >> The bug was caused by some horrendous Grobner basis reduction in a >> 2-variable polynomial ring over a number field. > > Is this a bug in Grobner basis or in my code that happens to use it? The latter. You coerce things into a ring of the form K[X,Y]/I where I is an ideal, in order to "reduce" them modulo the defining equation of the curve (which is the generator of I). That leads to a reduce(0 function in quotient_ring_element.py being called on the elements with respect to that ideal, which in turn leads to code in t-y_buchberger.py where the Error is raised. Since I is a principal ideal it should not be too hard to find its G. basis? > >> I was able to >> eliminate it by reworking William's full_division_polynomial() to work >> in the genuine function field of the curve. Thus, if K is the >> base_ring of the curve I construct in order the polynomial ring K[x], >> then its field of fractions K(x), then -- via the polynomial ring >> K(x)[Y] -- the field K(x)(y) where y^2=polynomial in x is the equation >> of the curve. (I am currently assuming that the curve is in short >> Weierstrass form as WIlliam did, though that wilkl be fixed in due >> course). Elements g of that field have two components g[0] and g[1], >> both rational functions in x, representing g[0]+y*g[1]. >> >> It is possible though slightly awkward to evaluate such a function g >> at a point P. What works is lift(g).subs(Y=P[1]).subs(x=P[0]) >> >> Now full_division_polynomial(m) returns an element g of that field, >> which if m is odd happens to be a polynomial in x (no denominator, no >> y term, i.e. g[1]==0) and if m is even happens to be y times a poly on >> x (no denominator and g[0]==0). >> >> Using that I reworked multiplication_by_m() as necessary and hence >> division_points(). doctests all pass, and the example which >> previously failed in #3383 now passes. >> >> Problem 1: Over Q, multiplication_by_m(10) takes ages -- over 3 >> minutes! -- even with the x_only=True option. The time is almost all >> spent in dividing one element of QQ(x) by another. That is just a >> simple division of two rational functions (polynomials of degrees 99 >> and 100 with no common factor) -- possibly the time is taken >> computing their gcd? I know that ZZ[x] operations are faster than >> QQ[x] but it would be very ugly to try to detect this special case >> since the code is otherwise generic and works over any field (char. >> not 2 or 3 until we stop using short W. equations). > > It seems like you need to investigate this a bit further. If it is really > gcd, then we could just improve QQ[x]'s gcd code (e.g., clear denoms > and call something in ZZ[x] for now; I don't know.) > It is definitely the gcd. I was able to avoid the slowdown by constructing the element of the fraction field L not by x=n/d (where n and d are the numer- and denominator), nor by x=L(n,d) but by x = sage.rings.fraction_field_element.FractionFieldElement(L,n,d,reduce=False) which forces the gcd not to be taken. I am pretty sure that at this point in the code we know that the polynomials are coprime anyway, so I would vote for the constructor of the form L(n,d) to have a flag reduce whcih defaults to True but can be set to False for situations like that. >> >> Problem/question 2: looking at the code in ell_generic.py I find that >> there are now *3* different implementations of division polys of >> various kinds. >> >> (1) E.division_pol = E.torsion_pol is by David Kohel (2005-04-25) and >> returns a polynomial in x alone, with a factor of the 2-division poly >> in x when m is even. >> (2) E.full_division_pol is William's new code (which can be asked to >> call the preceding when m is odd) which returns a polynomial in x and >> y, i.e. in the 2-variable poly ring over K, which happens to be of the >> form p(x) or y*p(x) where p is a poly in x alone (but belonging to >> K[x,y]). >> (3) There is code which is all commented out due to David Harvey >> (2006-09-24) consisting of 3 functions pseudo_torsion_polynomial(), >> multiple_x_numerator(() and multiple_x_denominator(). The first is >> the same as before for odd m but for even m just omits the factor of >> the 2-division poly. I like this implementation: it uses a cache, >> and also the user supplies their own 'x' which need not be an >> indeterminate. This is *very* useful since it is expensive to compute >> (say) the 10th division poly as a polynomial and then substitute >> avalue for x, it is much faster to evaluate as you go along. The last >
[sage-devel] Re: elliptic curve division polynomials
On Tue, Jun 10, 2008 at 7:34 AM, John Cremona <[EMAIL PROTECTED]> wrote: > > The following should be of interest to David Harvey, David Kohel, > Chris Wuthrich and WIlliam Stein (at least!). > > Following from my reviewing of Chris's #3377 (torsion subgroups of > elliptic curves over number fields) where I found a bug (see #3383) > in William's new full_division_polynomial() code (see #3109): > > The bug was caused by some horrendous Grobner basis reduction in a > 2-variable polynomial ring over a number field. Is this a bug in Grobner basis or in my code that happens to use it? > I was able to > eliminate it by reworking William's full_division_polynomial() to work > in the genuine function field of the curve. Thus, if K is the > base_ring of the curve I construct in order the polynomial ring K[x], > then its field of fractions K(x), then -- via the polynomial ring > K(x)[Y] -- the field K(x)(y) where y^2=polynomial in x is the equation > of the curve. (I am currently assuming that the curve is in short > Weierstrass form as WIlliam did, though that wilkl be fixed in due > course). Elements g of that field have two components g[0] and g[1], > both rational functions in x, representing g[0]+y*g[1]. > > It is possible though slightly awkward to evaluate such a function g > at a point P. What works is lift(g).subs(Y=P[1]).subs(x=P[0]) > > Now full_division_polynomial(m) returns an element g of that field, > which if m is odd happens to be a polynomial in x (no denominator, no > y term, i.e. g[1]==0) and if m is even happens to be y times a poly on > x (no denominator and g[0]==0). > > Using that I reworked multiplication_by_m() as necessary and hence > division_points(). doctests all pass, and the example which > previously failed in #3383 now passes. > > Problem 1: Over Q, multiplication_by_m(10) takes ages -- over 3 > minutes! -- even with the x_only=True option. The time is almost all > spent in dividing one element of QQ(x) by another. That is just a > simple division of two rational functions (polynomials of degrees 99 > and 100 with no common factor) -- possibly the time is taken > computing their gcd? I know that ZZ[x] operations are faster than > QQ[x] but it would be very ugly to try to detect this special case > since the code is otherwise generic and works over any field (char. > not 2 or 3 until we stop using short W. equations). It seems like you need to investigate this a bit further. If it is really gcd, then we could just improve QQ[x]'s gcd code (e.g., clear denoms and call something in ZZ[x] for now; I don't know.) > > Problem/question 2: looking at the code in ell_generic.py I find that > there are now *3* different implementations of division polys of > various kinds. > > (1) E.division_pol = E.torsion_pol is by David Kohel (2005-04-25) and > returns a polynomial in x alone, with a factor of the 2-division poly > in x when m is even. > (2) E.full_division_pol is William's new code (which can be asked to > call the preceding when m is odd) which returns a polynomial in x and > y, i.e. in the 2-variable poly ring over K, which happens to be of the > form p(x) or y*p(x) where p is a poly in x alone (but belonging to > K[x,y]). > (3) There is code which is all commented out due to David Harvey > (2006-09-24) consisting of 3 functions pseudo_torsion_polynomial(), > multiple_x_numerator(() and multiple_x_denominator(). The first is > the same as before for odd m but for even m just omits the factor of > the 2-division poly. I like this implementation: it uses a cache, > and also the user supplies their own 'x' which need not be an > indeterminate. This is *very* useful since it is expensive to compute > (say) the 10th division poly as a polynomial and then substitute > avalue for x, it is much faster to evaluate as you go along. The last > two give the numerator and denominator of the x-coordinate of m*P, > i.e. of the first rational function returned by William's > multiplication_by_m(). [By the way, it is easy to get the > y-coordinate from the x-coordinate: if the x-coord is X(x) then the > y-coord is c*y*X'(x) where c is a constant. I think the constant is m > in this case (look at the invariant differential). > > Does anyone know why those three functions are commented out? No (since David Harvey doesn't know). > William, why did you write your own new functions instead of using > these? Since I didn't know how to use them to do what I wanted. > Can we try to clean all this up once and for all? Would you > trust me to do this? Yes. Yes. william --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-devel] Re: elliptic curve division polynomials
2008/6/10 David Harvey <[EMAIL PROTECTED]>: > > I can't remember why they are commented out. I'll try removing the comments and see what happens! > > However I should mention that there is yet *another* implementation > of something very similar, see sage.schemes.elliptic_curves.padics, > the function _multiply_point(). I'll leave that well alone, since although there is some trivial code duplication that function has a very specific purpose. John > > david > > > > > --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-devel] Re: elliptic curve division polynomials
On Jun 10, 2008, at 10:34 AM, John Cremona wrote: > Problem/question 2: looking at the code in ell_generic.py I find that > there are now *3* different implementations of division polys of > various kinds. > > (1) E.division_pol = E.torsion_pol is by David Kohel (2005-04-25) and > returns a polynomial in x alone, with a factor of the 2-division poly > in x when m is even. > (2) E.full_division_pol is William's new code (which can be asked to > call the preceding when m is odd) which returns a polynomial in x and > y, i.e. in the 2-variable poly ring over K, which happens to be of the > form p(x) or y*p(x) where p is a poly in x alone (but belonging to > K[x,y]). > (3) There is code which is all commented out due to David Harvey > (2006-09-24) consisting of 3 functions pseudo_torsion_polynomial(), > multiple_x_numerator(() and multiple_x_denominator(). The first is > the same as before for odd m but for even m just omits the factor of > the 2-division poly. I like this implementation: it uses a cache, > and also the user supplies their own 'x' which need not be an > indeterminate. This is *very* useful since it is expensive to compute > (say) the 10th division poly as a polynomial and then substitute > avalue for x, it is much faster to evaluate as you go along. The last > two give the numerator and denominator of the x-coordinate of m*P, > i.e. of the first rational function returned by William's > multiplication_by_m(). [By the way, it is easy to get the > y-coordinate from the x-coordinate: if the x-coord is X(x) then the > y-coord is c*y*X'(x) where c is a constant. I think the constant is m > in this case (look at the invariant differential). > > Does anyone know why those three functions are commented out? > William, why did you write your own new functions instead of using > these? Can we try to clean all this up once and for all? Would you > trust me to do this? I can't remember why they are commented out. However I should mention that there is yet *another* implementation of something very similar, see sage.schemes.elliptic_curves.padics, the function _multiply_point(). david --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
