Opened 5 years ago
Last modified 3 months ago
#11975 needs_work enhancement
ChowHeegner points
Reported by:  was  Owned by:  cremona 

Priority:  major  Milestone:  sage6.9 
Component:  elliptic curves  Keywords:  chow heegner 
Cc:  Merged in:  
Authors:  William Stein  Reviewers:  John Cremona 
Report Upstream:  N/A  Work issues:  add tolerance 
Branch:  public/ticket/11975 (Commits)  Commit:  eb042f0ba5c6848a8311d66cd8da5ee8b8c90d48 
Dependencies:  Stopgaps: 
Description (last modified by )
Include an implementation in Sage of the algorithm of [1] for computation of ChowHeegner points on elliptic curves. Assure that it can replicate the table of [1].
[1] http://wstein.org/papers/chow_heegner/chowheeg1.pdf
Apply all of the patches in order. I have not folded these patches because I think the revision history is valuable.
Ignore the attached patches below. I've rebased these for sage5.10 and put the rebased patches here:
Attachments (6)
Change History (48)
comment:1 Changed 5 years ago by
 Milestone sage4.7.3 deleted
comment:2 Changed 5 years ago by
 Description modified (diff)
 Milestone set to sage4.8
comment:3 followup: ↓ 4 Changed 5 years ago by
There is a lot of very useful stuff here, which could be moved elsewhere otherwise there's a danger that others might reinvent some of the wheels. I am thinking of the SL2Z and Gamma0Nequivalence functions, and also the approximatepoint class which could also be used for ordinary Heegner points, and also for points computed from a complex number modulo a period lattice (socalled elliptic exponential). This is not a criticism. And neither is it a review yet  sorry!
comment:4 in reply to: ↑ 3 Changed 5 years ago by
Replying to cremona:
There is a lot of very useful stuff here, which could be
moved elsewhere otherwise there's a danger that others
might reinvent some of the wheels. I am thinking of the SL2Z and Gamma0Nequivalence functions, and also the approximatepoint class which could also be used for ordinary Heegner points, and also for points computed from a complex number modulo a period lattice (socalled elliptic exponential).
I realize that. However, I didn't want a patch that would constantly bitrot, and at the same time, have to write way more code than necessary (e.g., if I do Gamma0N, then might as well do Gamma1N, etc.) I think it would be good for this code to safely go in as a step 1, then in the future have small parts of it (the ones you list above) refactored and extended to more general cases. I think it will be easier to referee now since it is fairly self contained and well documented (by a short paper). The alternative might be 10 patches all over Sage with confusing dependencies....
This is not a criticism. And neither is it a review yet  sorry!
The code isn't quite ready for review, but will be within a day or two.
comment:5 Changed 5 years ago by
Quick note: I'm going to post a new version of this patch soon that accounts for numerical noise issues. I tried to use #10952 (the tolerance doctest framework patch), but  much to my surprise  it sadly seems totally useless when complex numbers are involved?!
comment:6 Changed 5 years ago by
 Status changed from new to needs_review
OK, this is now ready for review. My only concern is that it's very numerical, and I haven't tested it on any 32bit Linux boxes.
comment:7 Changed 5 years ago by
I am starting to review this now...
comment:8 Changed 5 years ago by
Here is part 1 of a review which will take a while.
 Applies fine to 4.8.rc0; two test failures (numerical noise, 64bit):
sage t long devel/sagemain/sage/schemes/elliptic_curves/chow_heegner_fast.pyx ********************************************************************** File "/home/jec/sage4.8.rc0/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner_fast.pyx", line 397: sage: max([abs(v[i]w[i]) for i in range(len(v))]) Expected: 3...e13 Got: 2.96849663677e13 **********************************************************************
sage t long devel/sagemain/sage/schemes/elliptic_curves/chow_heegner.py ********************************************************************** File "/home/jec/sage4.8.rc0/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner.py", line 1587: sage: phi(v[0]), phi(v[1]) # long time Expected: (1.02140518266e14 + 1.0*I, 5.55111512313e16 + 1.0*I) Got: (5.44009282066e15 + 1.0*I, 5.55111512313e16 + 1.0*I)
 No problems building docs since the new files are not added to the
reference manual! If they are (now or one day) there may well be issues with formatting, but I will not check for these now.
 Testing sl2z_rep_in_fundom(z):
sage: z=CDF(110.1,0.0000000001) sage: zz,g = sl2z_rep_in_fundom(z) sage: g.acton(z) 5684.1 + 99999999.6275*I sage: zz 0.341867712408 + 99999999.6769*I sage: g [ 56841 6258194] [ 10 1101]
This suggests rounding problems  which are always a problem with this algorithm especially when the input point has very small imaginary part. (And that does happen in practice!) My own strategy here (when Im(z) is very small) is to start by taking a good rational approximation to Re(z) and first use the SL(2,Z) element taking that rational to infinity; only then apply this algorithm.
 Gamma_0(N)equivalence: very clever idea to check that (1) z1,z2
and (2) Nz1,Nz2 are SL(2,Z)equivalent! I am still trying to prove to myself that this works even when z1,z2 are elliptic, since in this case the g taking z1 to z2 is not unique (in PSL(2,Z)). So is it not possible for both pairs to be SL(2,Z)equivalent but the matrix in SL(2,Z) taking Nz1 to Nz2 is not the conjugate by [N,0;0,1] of the one taking z1 to z2?
 Typo:
l.228 "do" > "due"
 That's it as far as line 303.
comment:9 Changed 5 years ago by
John,
The attached patch addresses (the hell out of?) all of your points.
 William
comment:10 Changed 5 years ago by
This is John's second review, which he couldn't post because of a DB error. I'm posting it for him:
 With the second patch, all tests pass except for one (see below) and
the docs build OK. I have not yet actually looked at the refman pages. And the new files have 100% coverage.
sage t devel/sagemain/sage/schemes/elliptic_curves/chow_heegner.py ********************************************************************** File "/home/jec/sage4.8.rc0/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner.py", line 236: sage: sl2z_rep_in_fundom(complex(2.17,.19)) Exception raised: Traceback (most recent call last): File "/home/jec/sage4.8.rc0/local/bin/ncadoctest.py", line 1231, in run_one_test self.run_one_example(test, example, filename, compileflags) File "/home/jec/sage4.8.rc0/local/bin/sagedoctest.py", line 38, in run_one_example OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags) File "/home/jec/sage4.8.rc0/local/bin/ncadoctest.py", line 1172, in run_one_example compileflags, 1) in test.globs File "<doctest __main__.example_2[28]>", line 1, in <module> sl2z_rep_in_fundom(complex(RealNumber('2.17'),RealNumber('.19')))###line 236: sage: sl2z_rep_in_fundom(complex(2.17,.19)) File "/home/jec/sage4.8.rc0/local/lib/python/sitepackages/sage/schemes/elliptic_curves/chow_heegner.py", line 240, in sl2z_rep_in_fundom if isinstance(z, (float, int, long)) or z.imag() <= 0: TypeError: 'float' object is not callable
The problem is that when z has type complex, z.imag is not a function so z.imag() raises an error. I don't know what to do: imag(z) works for CDF and complex, but for complex it returns a complex number!!!
sage: imag(CDF(1,1)) 1.0 sage: imag(complex(1,1)) (1+0j)
 In the comp method for X0NPoint you get an assertion error when the
levels are different. Would it not be better to do something different?
 Your AtkinLehner method for X0NPoint is more general than the doc
says: it says that q should be a prime power, but it works fine for any qN coprime to N/q (as it should!).
sage: z = CDF(1,1) sage: P = X0NPoint(z,90,30) # N = 90 = 10*9 sage: P.atkin_lehner(10) [0.11049724 + 0.00055248619*I] sage: P.atkin_lehner(10).atkin_lehner(10) [0.10554020 + 0.000013888696*I] sage: P.atkin_lehner(10).atkin_lehner(10) == P True
Since you'll be editing the INPUT anyway, I suggest adding an example such as this.
 Function close_points: (1) missing explanation of the search_bound
input param. (2) I guess you mean to use this function only for rank1 curves? You only use the first generator. And it's inefficient to compute n*Q for each n in the range [k..k] separately, which is *precisely* why Iwrote this a few years ago:
sage: multiples? File: /home/jec/sage4.8.rc0/local/lib/python2.6/sitepackages/sage/groups/generic.py Docstring: Return an iterator which runs through "P0+i*P" for "i" in "range(n)".
For example:
sage: P = EllipticCurve('389a1').gens()[0] sage: [Q[0] for Q in multiples(P,11,5*P)] [22625407/11397376, 47503/16641, 26/361, 10/9, 1, 0, 1, 10/9, 26/361, 47503/16641, 22625407/11397376]
 The last example for the identify function calls point_exact() not
identify()! And point_exact seems to be a methd of a different class. Did you mean to say c.numerical_approx().identify()?
 I checked that EllipticCurve?([0,2011]) is still not in the database
(no: it has conductor 1.7 billion). I won't ask for this example to be changed to [0,2012] (which has a smaller conductor!).
 check_optimal: special case needed for 990h3. sorry.....
 "Return the elliptic curve that this is the modular parametrization of."
should be "Return the elliptic curve of which this is the modular parametrization." (There are some lapses of grammar up with which I will not put!)
 reached line 1681 of the first patch, just before ChowHeegnerPoint? class 
comment:11 Changed 5 years ago by
The story so far: patches part2 and part3 address all the points Imade in the frist two parts of my review, with one exception. As I emailed William independently the test for Gamma0(N)equivalence of points in the upper halfplane is incorrect and can return True when it should be False, for example:
sage: N=5 sage: z1=CDF(2,1)/5 sage: z2=CDF(2,1)/5 sage: is_gamma0N_equivalent(z1,z2,5,1e10) True
This is wrong since (up to sign) the only matrices in SL(2,Z) taking z1 to z2 are 1,0],[4,1? and 2,1],[3,2?, neither of which is in Gamma0(5):
sage: g1=SL2Z([1,0,4,1]) sage: g2=SL2Z([2,1,3,2]) sage: abs(g1.acton(z2) z1) 6.20633538312e17 sage: abs(g2.acton(z2) z1) 2.00148302124e16
I'll now continue the review...
comment:12 Changed 5 years ago by
 Reviewers set to John Cremona
 Status changed from needs_review to needs_work
 Work issues set to Gamma0(N) equivalence
No more comments! As soon as the Gamma0(N) equivalence test is fixed I'll give the triple of patches a combined positive review.
Changed 5 years ago by
part 4, addressing the major issue John Cremona found in equivalence checking.
comment:13 followup: ↓ 14 Changed 5 years ago by
 Status changed from needs_work to needs_review
I fixed the Gamma0(N) equivalence test mainly by taking old code I had already written for a previous version and putting it back in. This ended up revealing some other issues, interesting examples, etc., which is why the patch is bigger than expected.
comment:14 in reply to: ↑ 13 Changed 5 years ago by
Replying to was:
I fixed the Gamma0(N) equivalence test mainly by taking old code I had already written for a previous version and putting it back in. This ended up revealing some other issues, interesting examples, etc., which is why the patch is bigger than expected.
I'll look at it as soon as I can.
comment:15 followup: ↓ 16 Changed 5 years ago by
 Status changed from needs_review to needs_info
The fix looks good to me. I am testing now (done & fine). I am slightly worried that the equality tests for i and rho using == would be better done using equality up to some precision?
One could cache the "constants" i, rho and their stabilisers, but it is hardly worth while since it will only be rare that these last lines of code are used in practice.
Lastly, do you need to test for both rho and rhobar given that you have normalised the fund.dom.reps.? If so, the third of your special cases will not happen.
sage: z = (2*CC(pi)*CC(I)/6).exp() sage: z 0.500000000000000 + 0.866025403784439*I sage: canonicalize_sl2z(z) (0.500000000000000 + 0.866025403784439*I, None)
comment:16 in reply to: ↑ 15 ; followup: ↓ 18 Changed 5 years ago by
Replying to cremona:
The fix looks good to me. I am testing now (done & fine). I am slightly worried that the equality tests for i and rho using == would be better done using equality up to some precision?
It is, in a sense, because it happens in C=ComplexField(prec)
, where prec is the parameter that gets passed in.
One could cache the "constants" i, rho and their stabilisers, but it is hardly worth while since it will only be rare that these last lines of code are used in practice.
Since the precision changes a lot, this might not be worth it, given how much other work happens and how rare this case is.
Lastly, do you need to test for both rho and rhobar given that you have normalised the fund.dom.reps.?
Nope  very good point. I've removed that.
Is there anything else you want me to change?
 William
Changed 5 years ago by
comment:17 Changed 5 years ago by
 Status changed from needs_info to needs_review
comment:18 in reply to: ↑ 16 Changed 5 years ago by
 Status changed from needs_review to positive_review
Replying to was:
Is there anything else you want me to change?
No! Here comes a positive review. You'll need to list the patches to be applied in the description (easy: all of them). I did fold them into one patch, but that made me the owner so I have not uploaded it.
 William
comment:19 Changed 5 years ago by
 Description modified (diff)
comment:20 Changed 5 years ago by
 Status changed from positive_review to needs_work
 Work issues Gamma0(N) equivalence deleted
Administrative note: trac_11975part2.patch needs a commit message.
comment:21 Changed 5 years ago by
There is numerical noise and a 32/64bit issue on hawk (OpenSolaris? 06.200932):
sage t long force_lib devel/sage/sage/schemes/elliptic_curves/chow_heegner_fast.pyx ********************************************************************** File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner_fast.pyx", line 469: sage: s.newton(.6*I) Expected: [(0.0720852069059  0.638326735148*I, 6, 1.3877787807814457e17)] Got: [(0.0720852069059  0.638326735148*I, 6, 0.0)] ********************************************************************** File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner_fast.pyx", line 473: sage: s.newton([0, .6*I]) Expected: [(0.605829586188, 7, 0.0), (0.0720852069059  0.638326735148*I, 6, 1.3877787807814457e17)] Got: [(0.605829586188, 7, 0.0), (0.0720852069059  0.638326735148*I, 6, 0.0)] ********************************************************************** sage t long force_lib devel/sage/sage/schemes/elliptic_curves/chow_heegner.py ********************************************************************** File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner.py", line 86: sage: zz Expected: 0.341867712408 + 99999999.6769*I Got: 0.34186771241 + 99999999.6769*I ********************************************************************** File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner.py", line 207: sage: sl2z_rep_in_fundom(CC(1+I), eps=1) # must be small Expected: Traceback (most recent call last): ... ValueError: prec (=0) must be >= 2 and <= 9223372036854775807. Got: Traceback (most recent call last): File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/local/bin/ncadoctest.py", line 1231, in run_one_test self.run_one_example(test, example, filename, compileflags) File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/local/bin/sagedoctest.py", line 38, in run_one_example OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags) File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/local/bin/ncadoctest.py", line 1172, in run_one_example compileflags, 1) in test.globs File "<doctest __main__.example_2[22]>", line 1, in <module> sl2z_rep_in_fundom(CC(Integer(1)+I), eps=Integer(1)) # must be small###line 207: sage: sl2z_rep_in_fundom(CC(1+I), eps=1) # must be small File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/local/lib/python/sitepackages/sage/schemes/elliptic_curves/chow_heegner.py", line 265, in sl2z_rep_in_fundom w = ComplexIntervalField(prec)(z) File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/local/lib/python/sitepackages/sage/rings/complex_interval_field.py", line 70, in ComplexIntervalField C = ComplexIntervalField_class(prec) File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/local/lib/python/sitepackages/sage/rings/complex_interval_field.py", line 144, in __init__ ParentWithGens.__init__(self, self._real_field(), ('I',), False, category = Fields()) File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/local/lib/python/sitepackages/sage/rings/complex_interval_field.py", line 246, in _real_field self.__real_field = real_mpfi.RealIntervalField(self._prec) File "real_mpfi.pyx", line 277, in sage.rings.real_mpfi.RealIntervalField (sage/rings/real_mpfi.c:3097) File "real_mpfi.pyx", line 443, in sage.rings.real_mpfi.RealIntervalField_class.__init__ (sage/rings/real_mpfi.c:3302) ValueError: prec (=0) must be >= 2 and <= 2147483647. ********************************************************************** File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner.py", line 949: sage: f(t[0][0]) Expected: 4.4408920985e16  1.11022302463e16*I Got: 4.4408920985e16  8.66548586359e17*I ********************************************************************** File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner.py", line 1524: sage: phi(v[1]) Expected: 1.0 Got: 1.0  8.49472402142e17*I ********************************************************************** File "/export/home/buildbot/build/sage/hawk1/hawk_full/build/sage5.0.beta3/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner.py", line 2011: sage: P.numerical_approx() Expected: (1.44444444444440...: 1.03703703703... : 1.00000000000000) Got: (1.44444444444447 + 2.74248939530208e15*I : 1.03703703703708  3.49918222885543e15*I : 1.00000000000000) **********************************************************************
comment:22 Changed 5 years ago by
Could you also change the headings of the newly added files to be compliant with http://sagemath.org/doc/developer/conventions.html#headingsofsagelibrarycodefiles
comment:23 Changed 5 years ago by
In module_list.py
, you should not use cblas
but the variables BLAS, BLAS2.
comment:24 Changed 5 years ago by
On cleo (RHEL 5.3 ia64):
sage t long force_lib devel/sage/sage/schemes/elliptic_curves/chow_heegner.py ********************************************************************** File "/home/buildbot/build/sage/cleo1/cleo_full/build/sage5.0.beta3/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner.py", line 949: sage: f(t[0][0]) Expected: 4.4408920985e16  1.11022302463e16*I Got: 4.4408920985e16  8.97307024383e17*I ********************************************************************** File "/home/buildbot/build/sage/cleo1/cleo_full/build/sage5.0.beta3/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner.py", line 395: sage: is_sl2z_equivalent(z1, z2, 40) Expected: True Got: False ********************************************************************** File "/home/buildbot/build/sage/cleo1/cleo_full/build/sage5.0.beta3/devel/sagemain/sage/schemes/elliptic_curves/chow_heegner.py", line 1524: sage: phi(v[1]) Expected: 1.0 Got: 1.0  1.00505377226e16*I **********************************************************************
comment:25 Changed 4 years ago by
 Description modified (diff)
comment:26 Changed 4 years ago by
 Description modified (diff)
comment:27 followup: ↓ 28 Changed 4 years ago by
Sorry but I am now confused about what to apply. Presumably the 5 patches on githhub in the order part1,...,part5, but the last of those still does not address Jeroen's points?
I am really hoping that I can give this another positive review without having to read all except one last patch! I had completely forgotten that this never made it. SO as soon as you tag it "needs review" I'll look at it.
comment:28 in reply to: ↑ 27 Changed 4 years ago by
Replying to cremona:
Sorry but I am now confused about what to apply. Presumably the 5 patches on githhub in the order part1,...,part5,
Yes.
but the last of those still does not address Jeroen's points?
Yes.
I am really hoping that I can give this another positive review without having to read all except one last patch! I had completely forgotten that this never made it. SO as soon as you tag it "needs review" I'll look at it.
It's not "needs review" yet. The issues are just some numerical noise, and I'm just updating the patches so people can use them (not so that they can get into sage, just yet).
comment:29 Changed 4 years ago by
 Milestone changed from sage5.11 to sage5.12
comment:30 Changed 3 years ago by
 Milestone changed from sage6.1 to sage6.2
comment:31 Changed 3 years ago by
 Branch set to u/haochen_uw/ticket/11975
comment:32 Changed 3 years ago by
 Commit set to 5fc9462412f1d24f194de1308c572205fec8b8e2
comment:33 Changed 3 years ago by
 Milestone changed from sage6.2 to sage6.3
comment:34 Changed 3 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:35 Changed 20 months ago by
 Branch changed from u/haochen_uw/ticket/11975 to public/ticket/11975
 Commit 5fc9462412f1d24f194de1308c572205fec8b8e2 deleted
 Keywords chow heegner added
 Milestone changed from sage6.4 to sage6.9
 Work issues set to add tolerance
comment:36 Changed 20 months ago by
 Commit set to 23328cd74700d0d916ce59dc00884e55dd88457d
Branch pushed to git repo; I updated commit sha1. New commits:
9d3b47c  rebased trac 11975: chowheegner points

3a11242  fixed the doctests of chow_heegner.py

5fc9462  fixed all doctests; fixed infinite loop issue in points_in_h

f67069b  Merge branch 'u/haochen_uw/ticket/11975' into 6.9.b3

23328cd  trac #11975 make it compile again

comment:37 Changed 20 months ago by
 Commit changed from 23328cd74700d0d916ce59dc00884e55dd88457d to 3b0956990a127e4f47418de715f0d72eb95a4ad9
Branch pushed to git repo; I updated commit sha1. New commits:
3b09569  trac #11975 full pep8 and numerical tolerance

comment:38 Changed 19 months ago by
 Commit changed from 3b0956990a127e4f47418de715f0d72eb95a4ad9 to 2fd9d8f020934d8008fdad45c19a38f2d89cbec3
Branch pushed to git repo; I updated commit sha1. New commits:
2fd9d8f  trac #11975 more abs tol added

comment:39 Changed 8 months ago by
 Commit changed from 2fd9d8f020934d8008fdad45c19a38f2d89cbec3 to decb3bdc51b9d49c8a3a1a9f8dc73314812987bb
comment:40 Changed 8 months ago by
 Commit changed from decb3bdc51b9d49c8a3a1a9f8dc73314812987bb to 322f01a651f3b953e608a89e377b185598be3f3f
Branch pushed to git repo; I updated commit sha1. New commits:
322f01a  trac 11975 one detail

comment:41 Changed 8 months ago by
 Commit changed from 322f01a651f3b953e608a89e377b185598be3f3f to 7a069fedfbc4880f6a1916fb37e56400d104cdd3
Branch pushed to git repo; I updated commit sha1. New commits:
7a069fe  trac 11975 using cysignals

comment:42 Changed 3 months ago by
 Commit changed from 7a069fedfbc4880f6a1916fb37e56400d104cdd3 to eb042f0ba5c6848a8311d66cd8da5ee8b8c90d48
Milestone sage4.7.3 deleted