PeriodicRegion.__div__: use integer arithmetic

[jdemeyer]

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.

CC: @JohnCremona

Component: elliptic curves

Author: Jeroen Demeyer

Branch/Commit: 15b7588

Reviewer: John Cremona

Issue created by migration from https://trac.sagemath.org/ticket/17088

[jdemeyer]

Author: Jeroen Demeyer

[jdemeyer]

Branch: u/jdemeyer/ticket/17088

[vbraun]

Commit: 15b7588

[vbraun]

comment:5

Looks good to me. John, do you want to review this?

---

New commits:

15b7588

Use integer arithmetic in PeriodicRegion.__div__

[jdemeyer]

comment:6

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.

[JohnCremona]

comment:7

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.

[jdemeyer]

comment:8

Replying to @JohnCremona:

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.

[jdemeyer]

comment:9

This patch does fix the problem at #16938 (but it can be reviewed independently of that).

[JohnCremona]

comment:10

Replying to @jdemeyer:

Replying to @JohnCremona:

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 rc 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 nn times to fill the whole parallelogram with n*n copies of the original configuration.

Back to the code...

[JohnCremona]

comment:11

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 (nr)x(nc) smaller tiles, then shrinks these by a factor of n. Before shrinking there are 2 integer indices, ar+i and bc+j running through range(nr) and range(nc) 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 ar+i is exactly divisible by n then using (ar+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

[jdemeyer]

comment:12

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.

[jdemeyer]

comment:13

It seems like the old code was confused about whether or not n should be an integer...

[JohnCremona]

comment:14

Thanks for taking the time to think about this!

The trouble is that before and after the transformation there are exactly rc tiles. If the new configuration had (nr)(nc) tiles then everything would be exact. As it is, there is a loss of resulution: there has to be because a pattern drawn on rc tiles is redrawn nn times, still on the same number rc of tiles. Think of this in two steps: first make (nr)(nc) minitiles which relpicate the original pattern nn times on rc blocks of minitiles, so minitile with coordinates (ar+i,bc+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.

[jdemeyer]

comment:15

Replying to @JohnCremona:

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.

[JohnCremona]

Reviewer: John Cremona

[JohnCremona]

comment:16

OK, I am convinced!

[vbraun]

Changed branch from u/jdemeyer/ticket/17088 to 15b7588