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 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.

Change History (17)

comment:1 Changed 6 years ago by jdemeyer

  • Authors set to Jeroen Demeyer
  • Description modified (diff)

comment:2 Changed 6 years ago by jdemeyer

  • Description modified (diff)

comment:3 Changed 6 years ago by jdemeyer

  • Description modified (diff)

comment:4 Changed 6 years ago by jdemeyer

  • 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 vbraun

  • Commit set to 15b7588e82c7d67bbbab14fec4c14c862bcc4c1f

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


New commits:

15b7588Use integer arithmetic in PeriodicRegion.__div__

comment:6 Changed 6 years ago by jdemeyer

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: Changed 6 years ago by cremona

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: Changed 6 years ago by 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.

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 jdemeyer

  • 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 cremona

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 cremona

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 jdemeyer

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 jdemeyer

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

comment:14 follow-up: Changed 6 years ago by cremona

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 jdemeyer

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 cremona

  • Reviewers set to John Cremona
  • Status changed from needs_review to positive_review

OK, I am convinced!

comment:17 Changed 6 years ago by vbraun

  • Branch changed from u/jdemeyer/ticket/17088 to 15b7588e82c7d67bbbab14fec4c14c862bcc4c1f
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.