_S_class_group_and_units is mathematically incorrect

[rharron]

The output of _S_class_group_and_units is incorrect, and hence the output of selmer_group as well, in some cases where S contains non-principal ideals. Here's an example:

sage: K.<a> = NumberField(x^3 - 381 * x + 127)
sage: S = tuple(K.primes_above(13))
sage: K.selmer_group(S, 2)
[-7/13*a^2 - 140/13*a + 36/13,
 14/13*a^2 + 267/13*a - 85/13,
 7/13*a^2 + 127/13*a - 49/13,
 -1,
 1/13*a^2 + 20/13*a - 7/13,
 1/13*a^2 - 19/13*a + 6/13,
 121,
 10/13*a^2 + 44/13*a - 4555/13]

It's fairly easy, using Sage, to see that the S-2-Selmer group of K is an 8-dimensional F_2-vector space. The list of length 8 that is returned is supposed to be a basis of this (or rather a set of representatives in K^×). But the S-2-Selmer group is a subgroup of K^× mod squares, so 121 is the zero vector and hence the output is not linearly independent. The problem lies in the following:

sage: K._S_class_group_and_units(S)
([-7/13*a^2 - 140/13*a + 36/13,
  14/13*a^2 + 267/13*a - 85/13,
  7/13*a^2 + 127/13*a - 49/13,
  -1,
  1/13*a^2 + 20/13*a - 7/13,
  1/13*a^2 - 19/13*a + 6/13],
 [(Fractional ideal (11, a - 2), 2, 121),
  (Fractional ideal (19, 1/13*a^2 - 45/13*a - 332/13),
   2,
   10/13*a^2 + 44/13*a - 4555/13)])

The 121 in there is supposed to be such that (11, a-2)² = (121). But (11, a-2)² is not principal (in fact, (11, a-2) is a generator of the cyclic subgroup of order 6). It is just in the subgroup of the class group generated by the primes in S.

I'll shortly upload a patch that fixes this (partially suggested to me by Zev Klagsbrun).

Apply: attachment: trac_14489_fix_S_class_group_and_units.patch, attachment: trac_14489_doctests_32_64_bit.patch

Component: number fields

Keywords: S-class group

Author: Robert Harron

Reviewer: Peter Bruin

Merged: sage-5.11.rc0

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

[rharron]

Author: Robert Harron

[rharron]

comment:1

Note: technically this depends on #14391, but that ticket has already been merged into 5.9beta4, so I haven't listed it.

Also, I doctested the number field module with this patch applied, but haven't done the full long doctesting. I'll do that overnight.

[tscrim]

comment:2

Some things on the docstrings and code from looking at the patch file:

• I think it would be better to put the import of prod and Matrix at the module level since they are imported at startup anyways. This will have a small, but noticable effect on the speed.
• Could you add a doctest to check the output when len(S) == 0 and/or J.is_principle() is true (if these tests aren't already present)?
• Could you indent the block starting at line 3486 by 2 spaces (the first char should be in-line with the first non-whitespace character after the bullet/dash)?

I can't do the mathematical review.

Thanks,

Travis

[rharron]

comment:3

Thanks, Travis, I'll try to get around to those changes tonight. One comment and one question:

(1) Re the doctests: the first doctest (which was already there) is an example of the case len(S) == 0. The second doctest (which was already there) is when S contains only principal ideals and so the J you come across is necessarily already principal. Also, in the doctest I added, for the second ideal (namely (19, 1/13a² - 45/13a - 332/13)) the J you get is principal. So, your second point is already covered. Should I add comments to this effect in the docs?

(2) Re import: So, I have no real idea when to import in the module and when to import in the function. Is there somewhere to read about that. I guess you're saying that in the case of the two functions at hand, these are imported by sage when you boot it and so really are already loaded and so I might as well import them at the module level?

Thanks!

[tscrim]

comment:4

Hey Robert,

Replying to @rharron:

(1) Re the doctests: the first doctest (which was already there) is an example of the case len(S) == 0. The second doctest (which was already there) is when S contains only principal ideals and so the J you come across is necessarily already principal. Also, in the doctest I added, for the second ideal (namely (19, 1/13a² - 45/13a - 332/13)) the J you get is principal. So, your second point is already covered. Should I add comments to this effect in the docs?

No, this is likely clear to those in the field.

(2) Re import: So, I have no real idea when to import in the module and when to import in the function. Is there somewhere to read about that.

It's a matter of when the module with the class/function is loaded into sage. If it's at the module level, it is imported when the module is imported (which is anytime any part of a module is imported [at startup, this is in an all.py file]). Thus there are two reasons not to import something at the module level:

1. It is a very large import that is not happening otherwise and is typically not needed until called. These large imports causes sage to take longer to startup (ex. numpy is not imported at startup).

2. It would cause a circular dependency, in which case sage will fail to start.

I don't know of a specific place offhand to read on this...

I guess you're saying that in the case of the two functions at hand, these are imported by sage when you boot it and so really are already loaded and so I might as well import them at the module level?

Correct, this will make the function (marginally) faster each time it's called.

Best,

Travis

[rharron]

comment:5

2. It would cause a circular dependency, in which case sage will fail to start.

That would explain why I was having some trouble. I was able to load prod at the module level, but whenever I put Matrix in sage wouldn't start.

[rharron]

comment:6

Running the long tests, I got a doctest failure in polynomial_quotient_ring.py (since it has a function _S_class_group_and_units that uses the function I changed to compute such things for quotients of polynomial rings over number fields by not necessarily irreducible polynomials):

**********************************************************************
File "/Users/rharron/sages/sage-5.8/devel/sage/sage/rings/polynomial/polynomial_quotient_ring.py", line 984:
    sage: S.S_class_group([])
Expected:
    [((1/4*xbar^2 + 31/4, (-1/8*a + 1/8)*xbar^2 - 31/8*a + 31/8, 1/16*xbar^3 + 1/16*xbar^2 + 31/16*xbar + 31/16, -1/16*a*xbar^3 + (1/16*a + 1/8)*xbar^2 - 31/16*a*xbar + 31/16*a + 31/8), 6, -1/16*xbar^3 + 1/16*xbar^2 - 31/16*xbar + 47/16), ((-1/4*xbar^2 - 23/4, (1/8*a - 1/8)*xbar^2 + 23/8*a - 23/8, -1/16*xbar^3 - 1/16*xbar^2 - 23/16*xbar - 23/16, 1/16*a*xbar^3 + (-1/16*a - 1/8)*xbar^2 + 23/16*a*xbar - 23/16*a - 23/8), 6, 1/16*xbar^3 + 3/16*xbar^2 + 23/16*xbar + 85/16), ((-5/4*xbar^2 - 115/4, 5/4*a*xbar^2 + 115/4*a, -5/16*xbar^3 + 5/16*xbar^2 - 115/16*xbar + 115/16, 1/16*a*xbar^3 + 7/16*a*xbar^2 + 23/16*a*xbar + 161/16*a), 2, 5/16*xbar^3 + 37/16*xbar^2 + 115/16*xbar + 867/16)]
Got:
    [((1/4*xbar^2 + 31/4, (-1/8*a + 1/8)*xbar^2 - 31/8*a + 31/8, 1/16*xbar^3 + 1/16*xbar^2 + 31/16*xbar + 31/16, -1/16*a*xbar^3 + (1/16*a + 1/8)*xbar^2 - 31/16*a*xbar + 31/16*a + 31/8), 6, 1/16*xbar^3 - 5/16*xbar^2 + 31/16*xbar - 139/16), ((-1/4*xbar^2 - 23/4, (1/8*a - 1/8)*xbar^2 + 23/8*a - 23/8, -1/16*xbar^3 - 1/16*xbar^2 - 23/16*xbar - 23/16, 1/16*a*xbar^3 + (-1/16*a - 1/8)*xbar^2 + 23/16*a*xbar - 23/16*a - 23/8), 6, 1/16*xbar^3 + 3/16*xbar^2 + 23/16*xbar + 85/16), ((-5/4*xbar^2 - 115/4, 5/4*a*xbar^2 + 115/4*a, -5/16*xbar^3 + 5/16*xbar^2 - 115/16*xbar + 115/16, 1/16*a*xbar^3 + 7/16*a*xbar^2 + 23/16*a*xbar + 161/16*a), 2, 5/16*xbar^3 + 37/16*xbar^2 + 115/16*xbar + 867/16)]
**********************************************************************
File "/Users/rharron/sages/sage-5.8/devel/sage/sage/rings/polynomial/polynomial_quotient_ring.py", line 990:
    sage: S.S_class_group([K.ideal(a)])
Expected:
    [((1/4*xbar^2 + 31/4, (-1/8*a + 1/8)*xbar^2 - 31/8*a + 31/8, 1/16*xbar^3 + 1/16*xbar^2 + 31/16*xbar + 31/16, -1/16*a*xbar^3 + (1/16*a + 1/8)*xbar^2 - 31/16*a*xbar + 31/16*a + 31/8), 6, -1/16*xbar^3 + 1/16*xbar^2 - 31/16*xbar + 47/16), ((-1/4*xbar^2 - 23/4, (1/8*a - 1/8)*xbar^2 + 23/8*a - 23/8, -1/16*xbar^3 - 1/16*xbar^2 - 23/16*xbar - 23/16, 1/16*a*xbar^3 + (-1/16*a - 1/8)*xbar^2 + 23/16*a*xbar - 23/16*a - 23/8), 2, -1/8*xbar^2 - 15/8)]
Got:
    [((1/4*xbar^2 + 31/4, (-1/8*a + 1/8)*xbar^2 - 31/8*a + 31/8, 1/16*xbar^3 + 1/16*xbar^2 + 31/16*xbar + 31/16, -1/16*a*xbar^3 + (1/16*a + 1/8)*xbar^2 - 31/16*a*xbar + 31/16*a + 31/8), 6, -1/16*xbar^3 + 1/16*xbar^2 - 31/16*xbar + 47/16), ((-1/4*xbar^2 - 23/4, (1/8*a - 1/8)*xbar^2 + 23/8*a - 23/8, -1/16*xbar^3 - 1/16*xbar^2 - 23/16*xbar - 23/16, 1/16*a*xbar^3 + (-1/16*a - 1/8)*xbar^2 + 23/16*a*xbar - 23/16*a - 23/8), 2, -1/16*a*xbar^3 + (-13/16*a + 1/8)*xbar^2 - 23/16*a*xbar - 299/16*a + 31/8)]
**********************************************************************

It is understandable that the output has changed given that it uses a my new version of the function for number fields. However(!) when I open sage and run the code in the doctest, I get the expected output of the doctest, not the new one. Why is the doctest giving a different result than running the same commands in sage?

[rharron]

comment:7

Ugh, my comment was just swallowed by trac.

Try again, but quicker: I figured out the problems.

Line 984: Both answers are correct. At some point, my code returns -pr whereas the old code returns pr. The discrepancy between the doctest illustrates the subtlety in this: if I open sage and enter all the commands in the doctest for S_class_group in order, then I get the same answer the doctest got; whereas just entering the commands from line 981 to 984 gives the answer the doctest expected. Anyway, I'm just going to update the doctest. (Is there an issue with the difference between 32 and 64 bit, I'm on a 64 bit machine).

Line 990: The expected answer is incorrect. The ideal in question when squared is not principal (in fact, the line 984 example says that the order is 6 not 2). So, I'll update the doctest with the new answer.

[rharron]

comment:8

Attachment: trac_14489_fix_S_class_group_and_units.patch.gz

Alright, new patch uploaded.

[rharron]

comment:9

The patch applies cleanly against 5.9 and passes long tests.

[tscrim]

comment:10

Hey Robert, one more comment about the doc, I think it might be better to use a LaTeX (single-letter) variable for the gen/order since they are used so much in an exponential form. Otherwise I'd suggest \mathrm{} (but \text{} is fine too, as long as it's used everywhere). Thanks.

[pjbruin]

Reviewer: Peter Bruin

[pjbruin]

comment:12

I would tend to give this a positive review. The patch fixes the problem, I believe the code is correct, all tests pass, and the documentation builds correctly. I have a few comments, though:

1. The changes in the docstrings are not very clearly written. The sentences are very long and contain awkward constructions ("a fractional ideal representative of the S-class group generator whose order (in the S-class group) is order"; "principal generator").

2. The new docstring of _S_class_group_and_units in number_field.py suggests that to obtain a principal ideal, gen^order can be multiplied by any fractional ideal J whose class is in the subgroup of the class group generated by ideals in S. However, the condition is more strict: J must be in the subgroup of the ideal group generated by ideals in S.

3. In general (not because of this patch), the code for computing Selmer groups is a bit convoluted. Conceptually, the computation of the generators for principal ideals of the form gen^order belongs in selmer_group, not _S_class_group_and_units. In my opinion it would be more correct to return, as S-class group generators, pairs (gen, order) instead of triples (gen, order, pr), and leave the computation of the principal ideal generators to selmer_group.

I made a patch on top of this one to address these issues. An additional reviewer would probably be needed, but I think it should not be too hard to review. The two patches together have a net effect of simplifying the Selmer group code somewhat.

Alternatively, I could file a new ticket for the additional changes, if the second patch would delay closing this ticket too long.

Any opinions about which option would be better?

[pjbruin]

Attachment: trac_14489_simplify_S_class_group_and_units.patch.gz

suggested enhancements; apply after trac_14489_fix_S_class_group_and_units.patch

[rharron]

comment:13

Thanks for reviewing this, Peter. In instances where sage is silently returning a mathematically incorrect answer, I try to create a minimal patch that fixes the problem to speed up the review process and get the fix in quickly. So, I'd suggest that since you think the patch I've written works properly, then you should accept this ticket and open a new ticket (which would be an enhancement ticket rather than a defect ticket) to clean up the code. My rationale is that the patch I wrote only changed a couple dozen lines of code and is fairly simple and yet has taken almost 2 months to be reviewed, so who knows how long this new patch would take to review? Makes sense?

(Also, note that the long sentences in the docstrings and the term "principal generator" were already present before my patch. My philosophy of minimal changes made me leave them in.)

[pjbruin]

comment:14

OK, let's do it like that.

[jdemeyer]

comment:16

On 32-bit systems:

sage -t devel/sage/sage/rings/polynomial/polynomial_quotient_ring.py
**********************************************************************
File "devel/sage/sage/rings/polynomial/polynomial_quotient_ring.py", line 987, in sage.rings.polynomial.polynomial_quotient_ring.PolynomialQuotientRing_generic.S_class_group
Failed example:
    S.S_class_group([])
Expected:
    [((1/4*xbar^2 + 31/4, (-1/8*a + 1/8)*xbar^2 - 31/8*a + 31/8, 1/16*xbar^3 + 1/16*xbar^2 + 31/16*xbar + 31/16, -1/16*a*xbar^3 + (1/16*a + 1/8)*xbar^2 - 31/16*a*xbar + 31/16*a + 31/8), 6, 1/16*xbar^3 - 5/16*xbar^2 + 31/16*xbar - 139/16), ((-1/4*xbar^2 - 23/4, (1/8*a - 1/8)*xbar^2 + 23/8*a - 23/8, -1/16*xbar^3 - 1/16*xbar^2 - 23/16*xbar - 23/16, 1/16*a*xbar^3 + (-1/16*a - 1/8)*xbar^2 + 23/16*a*xbar - 23/16*a - 23/8), 6, 1/16*xbar^3 + 3/16*xbar^2 + 23/16*xbar + 85/16), ((-5/4*xbar^2 - 115/4, 5/4*a*xbar^2 + 115/4*a, -5/16*xbar^3 + 5/16*xbar^2 - 115/16*xbar + 115/16, 1/16*a*xbar^3 + 7/16*a*xbar^2 + 23/16*a*xbar + 161/16*a), 2, 5/16*xbar^3 + 37/16*xbar^2 + 115/16*xbar + 867/16)]
Got:
    [((1/4*xbar^2 + 31/4, (-1/8*a + 1/8)*xbar^2 - 31/8*a + 31/8, 1/16*xbar^3 + 1/16*xbar^2 + 31/16*xbar + 31/16, -1/16*a*xbar^3 + (1/16*a + 1/8)*xbar^2 - 31/16*a*xbar + 31/16*a + 31/8), 6, -1/16*xbar^3 + 1/16*xbar^2 - 31/16*xbar + 47/16), ((-1/4*xbar^2 - 23/4, (1/8*a - 1/8)*xbar^2 + 23/8*a - 23/8, -1/16*xbar^3 - 1/16*xbar^2 - 23/16*xbar - 23/16, 1/16*a*xbar^3 + (-1/16*a - 1/8)*xbar^2 + 23/16*a*xbar - 23/16*a - 23/8), 6, 1/16*xbar^3 + 3/16*xbar^2 + 23/16*xbar + 85/16), ((-5/4*xbar^2 - 115/4, 5/4*a*xbar^2 + 115/4*a, -5/16*xbar^3 + 5/16*xbar^2 - 115/16*xbar + 115/16, 1/16*a*xbar^3 + 7/16*a*xbar^2 + 23/16*a*xbar + 161/16*a), 2, 5/16*xbar^3 + 37/16*xbar^2 + 115/16*xbar + 867/16)]
**********************************************************************
File "devel/sage/sage/rings/polynomial/polynomial_quotient_ring.py", line 1063, in sage.rings.polynomial.polynomial_quotient_ring.PolynomialQuotientRing_generic.class_group
Failed example:
    S.class_group()
Expected:
    [((2, -a + 1, 1/2*xbar + 1/2, -1/2*a*xbar + 1/2*a + 1), 6, -1/2*xbar + 3/2)]
Got:
    [((2, -a + 1, 1/2*xbar + 1/2, -1/2*a*xbar + 1/2*a + 1), 6, 1/2*xbar - 3/2)]
**********************************************************************

[pjbruin]

comment:17

I don't have access to a 32-bit machine on which to test this. The expected output and the actual output differ only in the third element of each triple, and they differ by a unit (which is -1 in the second case; I haven't checked in the first case). This is consistent with the fact that this element is a somewhat arbitrary generator of a principal ideal. This doctest failure should go away when applying both this patch and #14746, which eliminates this third element, among other things.

Perhaps the quickest fix would be to mark the output of these doctests as random in this patch and revert to non-random in the patch of #14746?

[jdemeyer]

comment:18

Or better: instead of # random, use # 32-bit and # 64-bit as explained in http://sagemath.org/doc/developer/conventions.html#further-conventions-for-automated-testing-of-examples

[pjbruin]

correct results for doctests on 32-bit machines

[pjbruin]

comment:19

Attachment: trac_14489_doctests_32_64_bit.patch.gz

Hopefully the fact that I submitted a simple patch to fix the 32-bit/64-bit doctest issue does not preclude me from giving a positive review again.

[pjbruin]

comment:20

For patchbot:

Apply trac_14489_fix_S_class_group_and_units.patch, trac_14489_doctests_32_64_bit.patch

[vbraun]

comment:21

Test failure in my 32-bit VM:

[vbraun@localhost hg]$ sage -t sage/rings/polynomial/polynomial_quotient_ring.py  # 2 doctests failed
Running doctests with ID 2013-07-18-21-23-12-2dbc373e.
Doctesting 1 file.
sage -t sage/rings/polynomial/polynomial_quotient_ring.py
**********************************************************************
File "sage/rings/polynomial/polynomial_quotient_ring.py", line 987, in sage.rings.polynomial.polynomial_quotient_ring.PolynomialQuotientRing_generic.S_class_group
Failed example:
    S.S_class_group([])
Expected:
    [((1/4*xbar^2 + 31/4, (-1/8*a + 1/8)*xbar^2 - 31/8*a + 31/8, 1/16*xbar^3 + 1/16*xbar^2 + 31/16*xbar + 31/16, -1/16*a*xbar^3 + (1/16*a + 1/8)*xbar^2 - 31/16*a*xbar + 31/16*a + 31/8), 6, 1/16*xbar^3 - 5/16*xbar^2 + 31/16*xbar - 139/16), ((-1/4*xbar^2 - 23/4, (1/8*a - 1/8)*xbar^2 + 23/8*a - 23/8, -1/16*xbar^3 - 1/16*xbar^2 - 23/16*xbar - 23/16, 1/16*a*xbar^3 + (-1/16*a - 1/8)*xbar^2 + 23/16*a*xbar - 23/16*a - 23/8), 6, 1/16*xbar^3 + 3/16*xbar^2 + 23/16*xbar + 85/16), ((-5/4*xbar^2 - 115/4, 5/4*a*xbar^2 + 115/4*a, -5/16*xbar^3 + 5/16*xbar^2 - 115/16*xbar + 115/16, 1/16*a*xbar^3 + 7/16*a*xbar^2 + 23/16*a*xbar + 161/16*a), 2, 5/16*xbar^3 + 37/16*xbar^2 + 115/16*xbar + 867/16)]
Got:
    [((1/4*xbar^2 + 31/4, (-1/8*a + 1/8)*xbar^2 - 31/8*a + 31/8, 1/16*xbar^3 + 1/16*xbar^2 + 31/16*xbar + 31/16, -1/16*a*xbar^3 + (1/16*a + 1/8)*xbar^2 - 31/16*a*xbar + 31/16*a + 31/8), 6, -1/16*xbar^3 + 1/16*xbar^2 - 31/16*xbar + 47/16), ((-1/4*xbar^2 - 23/4, (1/8*a - 1/8)*xbar^2 + 23/8*a - 23/8, -1/16*xbar^3 - 1/16*xbar^2 - 23/16*xbar - 23/16, 1/16*a*xbar^3 + (-1/16*a - 1/8)*xbar^2 + 23/16*a*xbar - 23/16*a - 23/8), 6, 1/16*xbar^3 + 3/16*xbar^2 + 23/16*xbar + 85/16), ((-5/4*xbar^2 - 115/4, 5/4*a*xbar^2 + 115/4*a, -5/16*xbar^3 + 5/16*xbar^2 - 115/16*xbar + 115/16, 1/16*a*xbar^3 + 7/16*a*xbar^2 + 23/16*a*xbar + 161/16*a), 2, 5/16*xbar^3 + 37/16*xbar^2 + 115/16*xbar + 867/16)]
**********************************************************************
File "sage/rings/polynomial/polynomial_quotient_ring.py", line 1063, in sage.rings.polynomial.polynomial_quotient_ring.PolynomialQuotientRing_generic.class_group
Failed example:
    S.class_group()
Expected:
    [((2, -a + 1, 1/2*xbar + 1/2, -1/2*a*xbar + 1/2*a + 1), 6, -1/2*xbar + 3/2)]
Got:
    [((2, -a + 1, 1/2*xbar + 1/2, -1/2*a*xbar + 1/2*a + 1), 6, 1/2*xbar - 3/2)]
**********************************************************************
2 items had failures:
   1 of  23 in sage.rings.polynomial.polynomial_quotient_ring.PolynomialQuotientRing_generic.S_class_group
   1 of  25 in sage.rings.polynomial.polynomial_quotient_ring.PolynomialQuotientRing_generic.class_group
    [357 tests, 2 failures, 6.25 s]
----------------------------------------------------------------------
sage -t sage/rings/polynomial/polynomial_quotient_ring.py  # 2 doctests failed
----------------------------------------------------------------------
Total time for all tests: 6.4 seconds
    cpu time: 5.8 seconds
    cumulative wall time: 6.2 seconds

[vbraun]

comment:22

The doctest failure seems to be fixed in the follow-up #14746, however that introduces another 32-bit failure.

[vbraun]

Dependencies: #14746

[vbraun]

comment:24

Sorry, just noticed that I was one off in my patch queue. Never mind, all tests pass.

[vbraun]

Changed dependencies from #14746 to none

[jdemeyer]

Merged: sage-5.11.rc0