Opened 6 years ago
Closed 6 years ago
#17088 closed defect (fixed)
PeriodicRegion.__div__: use integer arithmetic
Reported by: | jdemeyer | Owned by: | |
---|---|---|---|
Priority: | critical | Milestone: | sage-6.4 |
Component: | elliptic curves | Keywords: | |
Cc: | cremona | Merged in: | |
Authors: | Jeroen Demeyer | Reviewers: | John Cremona |
Report Upstream: | N/A | Work issues: | |
Branch: | 15b7588 (Commits) | Commit: | 15b7588e82c7d67bbbab14fec4c14c862bcc4c1f |
Dependencies: | Stopgaps: |
Description (last modified by )
Replace the unreliable floating-point arithmetic of PeriodicRegion.__div__
in src/sage/schemes/elliptic_curves/period_lattice_region.pyx
by pure integer arithmetic to avoid floating-point rounding issues.
Note that I don't understand the purpose of the code, I just want to fix the fact that the results are not reproducible.
Background
this is causing doctest failures in #16938 due to slightly different floating-point rounding. The failure is the following:
sage -t src/sage/schemes/elliptic_curves/period_lattice_region.pyx ********************************************************************** File "src/sage/schemes/elliptic_curves/period_lattice_region.pyx", line 395, in sage.schemes.elliptic_curves.period_lattice_region.PeriodicRegion.__div__ Failed example: (S / 3).plot() Expected: Graphics object consisting of 109 graphics primitives Got: Graphics object consisting of 121 graphics primitives **********************************************************************
(this output of plot()
commands was added by #16746). As you can see, there are more pixels in the debug version.
Change History (17)
comment:1 Changed 6 years ago by
- Description modified (diff)
comment:2 Changed 6 years ago by
- Description modified (diff)
comment:3 Changed 6 years ago by
- Description modified (diff)
comment:4 Changed 6 years ago by
- Branch set to u/jdemeyer/ticket/17088
- Created changed from 10/02/14 09:31:42 to 10/02/14 09:31:42
- Modified changed from 10/02/14 09:34:57 to 10/02/14 09:34:57
comment:5 Changed 6 years ago by
- Commit set to 15b7588e82c7d67bbbab14fec4c14c862bcc4c1f
comment:6 Changed 6 years ago by
I haven't set the ticket to needs_review because I haven't tested that it actually fixes the issue at #16938.
But you can certainly review the code changes itself.
comment:7 follow-up: ↓ 8 Changed 6 years ago by
OK, I am looking at the code. Note that this was implemented by Robert Bradshaw based on (a small) part of a thesis by a student of mine (published at http://www.ams.org/journals/mcom/2010-79-272/S0025-5718-10-02352-5/). I worked hard to get RB's code into acceptable shape as regards doctests etc. The application for this is in finding a lower bound for the canonical height of non-torsion points on elliptic curves over number fields, and this particular bit of code is dealing with the worst case, a complex embedding, with period parallelograms being recursively divided into sub-parallelogams to solve a minimization problem. It is by far the nastiest part of the code.
At the moment I cannot see how the new code manages to do in one line what the old code took 4 lines to do, but Jeroen is very clever so I will look again.
comment:8 in reply to: ↑ 7 ; follow-up: ↓ 10 Changed 6 years ago by
Replying to cremona:
At the moment I cannot see how the new code manages to do in one line what the old code took 4 lines to do, but Jeroen is very clever so I will look again.
Part of the problem with the old code is that it's not clear what it is supposed to do.
Anyway, I claim the following: if you assume that all the floating point operations are exact, the 4 lines all do exactly the same so they can easily be replaced by 1 line. I guess the addition of the 3 extra lines was simply to hide floating-point rounding errors.
comment:9 Changed 6 years ago by
- Status changed from new to needs_review
This patch does fix the problem at #16938 (but it can be reviewed independently of that).
comment:10 in reply to: ↑ 8 Changed 6 years ago by
Replying to jdemeyer:
Replying to cremona:
At the moment I cannot see how the new code manages to do in one line what the old code took 4 lines to do, but Jeroen is very clever so I will look again.
Part of the problem with the old code is that it's not clear what it is supposed to do.
I agree. It took me doxens of hours (literally) to work out what RB's code did, and I asure you that this whole file is now very much more understandable than it was before. That was a while ago though.
Anyway, I claim the following: if you assume that all the floating point operations are exact, the 4 lines all do exactly the same so they can easily be replaced by 1 line. I guess the addition of the 3 extra lines was simply to hide floating-point rounding errors.
Something like that, but I am now unconvinced in the correctness of either the old or the new.
There's a parallelogram, whose vertices in C are 0, w1, w2, w1+w2. The state of the object at any time is held by a 2-dimensional array "data", say of dimensions r x c, indicating, out of all the r*c sub-parallelograms which the big one is divided into, those for which data[i,j]=True have some property. The plots show this configuration. When the structure is divided by n the entire thing is supposed to be shrunk by 1/n and then replicated n*n times to fill the whole parallelogram with n*n copies of the original configuration.
Back to the code...
comment:11 Changed 6 years ago by
I just looked at the original Magma code but it has been changed too much by RB to be of any help here.
For rigour it is important the overestimate, in the sense that if any part of a tile contains a point for which the condition is True then the whole tile should be included. If one overestimates and includes a tile unnecessarily, then this may slow down the procedure a little but rigour is preserved. The original subdivision is into r x c tiles. Division first creates (n*r)x(n*c) smaller tiles, then shrinks these by a factor of n. Before shrinking there are 2 integer indices, a*r+i and b*c+j running through range(n*r) and range(n*c) respectively. [r and c are called rows and cols in the code.] These then undergo rounded division to give a pair of indices in range(r), range(c).
Here is the important part: when a*r+i is exactly divisible by n then using (a*r+i)/n is fine, and similarly with column indices, BUT if a*r+i is between two multiples of n then BOTH the quotients (rounded down AND up) should be used. That is why up to 4 tiles may get tagged.
The old code was supposed to do this, but I don't think it did so correctly. The following is, I think, what it should be like:
u = (a*rows+i)//n v = (b*cols+j)//n u_exact = (n*u==a*rows+i) v_exact = (n*v==b*cols+j) dat = data[i,j] new_data[u,v] = dat if not u_exact: new_data[u+1,v] = dat if not v_exact: new_data[u+1,v+1] = dat else: if not v_exact: new_data[u,v+1] = dat
comment:12 Changed 6 years ago by
I think I understand your explanation, but shrinking a tile by a factor n (where n is integer, this is crucial!) should ensure that a shrinked tile lies in exactly one big tile: every big tile contains exactly n x n shrinked tiles.
So I don't think there is a problem.
comment:13 Changed 6 years ago by
It seems like the old code was confused about whether or not n should be an integer...
comment:14 follow-up: ↓ 15 Changed 6 years ago by
Thanks for taking the time to think about this!
The trouble is that before and after the transformation there are exactly r*c tiles. If the new configuration had (n*r)*(n*c) tiles then everything would be exact. As it is, there is a loss of resulution: there has to be because a pattern drawn on r*c tiles is redrawn n*n times, still on the same number r*c of tiles. Think of this in two steps: first make (n*r)*(n*c) minitiles which relpicate the original pattern n*n times on r*c blocks of minitiles, so minitile with coordinates (a*r+i,b*c+j) gets the value of data[i,j] with a,b in range(n), i in range(r), j in range(c). Next we rearrange the coordinate blocks: instead of n identical blocks of r we have r blocks of n (in the x-coordinate and similarly for the y). These new blocks make up the new tiles, each a union of n*n minitiles with (in general) different values, and we must give the new tile the value True if any of those minitiles was True.
Enough for one day.
comment:15 in reply to: ↑ 14 Changed 6 years ago by
Replying to cremona:
These new blocks make up the new tiles, each a union of n*n minitiles
That's the crucial point here. Each minitile lies in exactly one big block, which is why it works.
comment:16 Changed 6 years ago by
- Reviewers set to John Cremona
- Status changed from needs_review to positive_review
OK, I am convinced!
comment:17 Changed 6 years ago by
- Branch changed from u/jdemeyer/ticket/17088 to 15b7588e82c7d67bbbab14fec4c14c862bcc4c1f
- Resolution set to fixed
- Status changed from positive_review to closed
Looks good to me. John, do you want to review this?
New commits:
Use integer arithmetic in PeriodicRegion.__div__