Opened 11 years ago

Closed 4 years ago

# reduced_basis for number field multiples wrong

Reported by: Owned by: schilly davidloeffler major sage-6.4 number fields John Cremona Jeroen Demeyer N/A 8f1f51d 8f1f51d83fc198a31683740001bcffe9a3d30bbc todo

This was reported by the "report a problem" link:

The reduced_basis for number fields applies LLL to basis of the maximal order in order to get a reduced basis. The problem is that after applying LLL, it multiples the transformation matrix by the vector [1,a,a2,...], where 'a' is generator of the field - instead of using the vector it actually used for the maximal order.

How to reproduce:

```x = QQ['x'].0
k = NumberField(x^6 + 2218926655879913714112*x^4 - \
32507675650290949030789018433536*x^3 + \
4923635504174417014460581055002374467948544*x^2 - \
36066074010564497464129951249279114076897746988630016*x + \
264187244046129768986806800244258952598300346857154900812365824,'a')
new_basis = k.reduced_basis()
print new_basis[0].minpoly()
```

This prints a crazy big polynomial. Looking a bit at what is returned it is clear that reduced_basis is multiplying the transformation matrix by the wrong vector.

To solve you just need to multiply by the vector of generators of the maximal order:

```mp = pari( k.Minkowski_embedding(k.maximal_order().gens(), prec=120) )
ml = sage.libs.pari.gen.gen.qflll(mp)
T = matrix([[ZZ(y) for y in list(x)] for x in list(ml)]).transpose()
r = Matrix(O.gens())* T
for rr in r[0]:
print rr.minpoly()
```

This works fine.

### comment:1 Changed 11 years ago by fwclarke

Looks like the definition of `O` was omitted. The following slightly simplified code works ok for me:

```sage: O = k.ring_of_integers()
sage: mp = pari( k.Minkowski_embedding(O.gens(), prec=120) )
sage: ml = sage.libs.pari.gen.gen.qflll(mp)
sage: T = ml.sage()
sage: r = vector(O.gens())*T
sage: for rr in r:
....:     rr.minpoly()
```

### comment:2 Changed 8 years ago by jdemeyer

• Milestone changed from sage-5.11 to sage-5.12

### comment:3 Changed 8 years ago by vbraun_spam

• Milestone changed from sage-6.1 to sage-6.2

### comment:4 Changed 8 years ago by vbraun_spam

• Milestone changed from sage-6.2 to sage-6.3

### comment:5 Changed 7 years ago by vbraun_spam

• Milestone changed from sage-6.3 to sage-6.4

### comment:6 Changed 6 years ago by jakobkroeker

• Stopgaps set to todo

### comment:7 Changed 4 years ago by cremona

While the original diagnosis is not wrong, there is a more basic problem as this smaller example shows:

```sage: K.<a> = NumberField(x^3-10)
sage: ZK = K.ring_of_integers()
sage: ZK.basis()
[1/3*a^2 + 1/3*a + 1/3, a, a^2]
sage: ZK([1,2,3])
3*a^2 + 2*a + 1
```

Although the ring of integers is an order which knows its own Z-basis, when you construct an element of the order from an integer vector it ignores that and uses the power_basis of the field and not its integral basis.

I think the problem is in these lines of _element_constructor() (order.py lines 1142-1143):

```       if not is_Element(x) or x.parent() is not self._K:
x = self._K(x)
```

There appears to be a related problem in the constructor of OrderElementAbsolute? in number_field_element:

```       K = order.number_field()
NumberFieldElement_absolute.__init__(self, K, f)
self._number_field = K
```

where the parameter f is not documented and there are no relevant tests, but this seems to be saying "take the data in f and construct a field element from it".

Perhaps the simplest solution is to go to those lines in order.py and insert code to handle the case where the input data is an integer vector of the right length.

### comment:8 Changed 4 years ago by cremona

After doing the above "simplest solution", which did fix the `x^3-10` example above, I saw that it's not the fix needed for the original which is to take the vectors output by the LLL reduction and push them into the ring of integers rather than the field. Now I get

```sage: [c.minpoly() for c in new_basis]

[x - 1,
x^6 + 3*x^5 + 258*x^4 + 511*x^3 + 3564*x^2 + 3309*x + 2347,
x^6 - 24*x^5 + 6126*x^4 - 312664*x^3 + 5407566*x^2 - 33643572*x + 95921443,
x^6 + 18*x^5 + 3366*x^4 + 82008*x^3 + 886962*x^2 - 840726*x + 5521647,
x^6 + 27*x^5 + 9579*x^4 + 623358*x^3 + 5060091*x^2 - 139224285*x + 880944177,
x^6 - 72*x^5 + 65286*x^4 - 10762768*x^3 + 473072922*x^2 - 2502686322*x + 54227921641]
```

which is a lot better.

Patch up when I have added doctests...

### comment:9 Changed 4 years ago by cremona

• Authors set to John Cremona
• Branch set to u/cremona/10017
• Commit set to 6bee3fca54bf50cfbd7a82e71d899b7ddafacb89
• Status changed from new to needs_review

The doctest based on the original post works fine for me when run manually, but when run as a doctest something weird happens -- it takes ages and then gives an enormous python crash. I have no idea what this is but hope that a human can help before this breaks all the patchbots out there...

New commits:

 ​6bee3fc `#10017: fix construction of order elements from list, and number field reduced basis`

### comment:10 Changed 4 years ago by git

• Commit changed from 6bee3fca54bf50cfbd7a82e71d899b7ddafacb89 to 5d51b3260b6bec7eacc24aca4f469ee1ba8b4572

Branch pushed to git repo; I updated commit sha1. New commits:

 ​0000be4 `Merge branch 'develop' into 10017` ​5d51b32 `#10017 added 2 lines to new doctest`

### comment:11 Changed 4 years ago by jdemeyer

Quick comments (more might come):

1. Replace `See (:trac:`10017`)::` by `Check that :trac:`10017` is fixed::` or something along those lines.
1. Replace
```sage: assert R.discriminant()==k.discriminant()
```

by

```sage: R.discriminant() == k.discriminant()
True
```

### comment:12 Changed 4 years ago by git

• Commit changed from 5d51b3260b6bec7eacc24aca4f469ee1ba8b4572 to 8b3eb02ede88e3de6bf03cc54f8565b2a0ef3bcf

Branch pushed to git repo; I updated commit sha1. New commits:

 ​8b3eb02 `#10017 minor changes to docstring`

### comment:13 Changed 4 years ago by jdemeyer

In this block of code:

```            M = self.Minkowski_embedding(self.integral_basis(), prec=prec)
T = pari(M).qflll().sage()
ZK = self.ring_of_integers()
self.__reduced_basis = [ ZK(v.list()) for v in T.columns() ]
```

you are assuming that `self.integral_basis()` corresponds to the basis chosen in `ZK(v.list())`. After all, that's the bug on this ticket, right? To me, this looks way too implicit. I would prefer to see the code written in such a way that you are actually using `self.integral_basis()` in this computation, instead of the `ZK()` call. I.e. something like

```ZK = self.integral_basis()
reduced_basis = [sum(x * g for x, g in zip(v.list(), ZK)) for v in T.columns()]
```

### comment:14 follow-up: ↓ 16 Changed 4 years ago by cremona

Yes, you are right, and perhaps being explicit like this is safest. I still think (do you agree?) that being able to do ZK(list) is useful and must give the linear combination of ZK's basis with coefficients from the list.

A related different point is that after computing the reduced basis, it is stored as `self.__reduced_basis`, but not actually used anywhere. Perhaps we should implement the ability to change the integral basis itself; but that might be too dangerous since various functions using pari directly will assume that the integral basis does not change.

Last edited 4 years ago by cremona (previous) (diff)

### comment:15 Changed 4 years ago by vdelecroix

• Description modified (diff)

(fixing exponentiation in ticket description)

### comment:16 in reply to: ↑ 14 Changed 4 years ago by jdemeyer

Yes, you are right, and perhaps being explicit like this is safest. I still think (do you agree?) that being able to do ZK(list) is useful and must give the linear combination of ZK's basis with coefficients from the list.

Sure, why not.

A related different point is that after computing the reduced basis, it is stored as `self.__reduced_basis`, but not actually used anywhere.

I would be fine to not store it.

Perhaps we should implement the ability to change the integral basis itself; but that might be too dangerous since various functions using pari directly will assume that the integral basis does not change.

I agree with the latter. Better not do this.

### comment:17 Changed 4 years ago by cremona

The problems I was having with testing were due to a cut-and-paste error: in the patch the coefficient of x starts 3666 instead of 36066 which makes the polynomial rather harder to work with. I will fix that and the other issues from the previous comment now.

### comment:18 Changed 4 years ago by git

• Commit changed from 8b3eb02ede88e3de6bf03cc54f8565b2a0ef3bcf to 6ee081a5d2052f1adcfb44fbba9afb29c2b7be56

Branch pushed to git repo; I updated commit sha1. New commits:

 ​6ee081a `Merge branch 'develop' into 10017`

### comment:19 Changed 4 years ago by git

• Commit changed from 6ee081a5d2052f1adcfb44fbba9afb29c2b7be56 to a0448a21aea340aa6dc4be2cbca17fc6cf756ab5

Branch pushed to git repo; I updated commit sha1. New commits:

 ​a0448a2 `#10017: fix new example and simplify code in line with review comments`

### comment:20 Changed 4 years ago by cremona

New changes: code simplified since there's really only 2 lines different between the two cases (totally real or not). Removed caching of reduced basis (was never used). Follow suggestion of comments 13 and following.

Needs review again!

### comment:21 Changed 4 years ago by jdemeyer

Style issue: better write

```if isinstance(x, (tuple, list)):
```

```if isinstance(x,list) or isinstance(x,tuple):
```

### comment:22 Changed 4 years ago by git

• Commit changed from a0448a21aea340aa6dc4be2cbca17fc6cf756ab5 to 39fc2278ce182b46715ffd4fb2071178a4c11092

Branch pushed to git repo; I updated commit sha1. New commits:

 ​39fc227 `#10017: small stylistic improvement`

### comment:23 Changed 4 years ago by jdemeyer

• Reviewers set to Jeroen Demeyer
• Status changed from needs_review to needs_work

The patchbot reports a failure on 32-bit:

```sage -t --long src/sage/rings/number_field/number_field.py
**********************************************************************
File "src/sage/rings/number_field/number_field.py", line 5599, in sage.rings.number_field.number_field.NumberField_generic.reduced_basis
Failed example:
[c.minpoly() for c in new_basis]
Expected:
[x - 1,
x^6 + 3*x^5 + 258*x^4 + 511*x^3 + 3564*x^2 + 3309*x + 2347,
x^6 - 24*x^5 + 6126*x^4 - 312664*x^3 + 5407566*x^2 - 33643572*x + 95921443,
x^6 + 18*x^5 + 3366*x^4 + 82008*x^3 + 886962*x^2 - 840726*x + 5521647,
x^6 + 27*x^5 + 9579*x^4 + 623358*x^3 + 5060091*x^2 - 139224285*x + 880944177,
x^6 - 72*x^5 + 65286*x^4 - 10762768*x^3 + 473072922*x^2 - 2502686322*x + 54227921641]
Got:
[x - 1,
x^6 - 45*x^5 + 6993*x^4 - 519048*x^3 + 8139204*x^2 + 64936080*x + 394386192,
x^6 + 3*x^5 + 258*x^4 + 511*x^3 + 3564*x^2 + 3309*x + 2347,
x^6 - 12*x^5 + 459*x^4 - 13068*x^3 + 171423*x^2 - 1008720*x + 2218941,
x^6 - 72*x^5 + 21618*x^4 - 2725154*x^3 + 96844941*x^2 - 367602102*x + 5299367596,
x^6 - 36*x^5 + 4104*x^4 - 225746*x^3 + 4831632*x^2 - 39603546*x + 130361299]
**********************************************************************
```

### comment:24 Changed 4 years ago by cremona

Would it be OK to tag the current output as 64-bit and this alternative as 32-bit? Of course this "reduced basis" is in no way canonical.

### comment:25 Changed 4 years ago by jdemeyer

let me check if the `prec` argument makes a difference.

### comment:26 Changed 4 years ago by jdemeyer

On 64-bit, the output stabilizes for `prec >= 117`; the result it

```[x - 1,
x^2 - x + 1,
x^6 + 3*x^5 - 102*x^4 - 103*x^3 + 10572*x^2 - 59919*x + 127657,
x^6 - 3*x^5 - 102*x^4 + 315*x^3 + 10254*x^2 - 80955*x + 198147,
x^3 - 171*x + 848,
x^6 + 171*x^4 + 1696*x^3 + 29241*x^2 + 145008*x + 719104]
```

### comment:27 Changed 4 years ago by cremona

That's a lot better -- can you check if it's the same on 32-bit using the same prec (say 120)?

### comment:28 Changed 4 years ago by jdemeyer

Yes, on 32-bit I get the same result for `prec >= 68`.

So the best solution is increasing the `prec` in the doctest. While you're at it, could you please replace `k = NumberField(...,'a')` by `k.<a> = NumberField(...)`?

### comment:29 Changed 4 years ago by git

• Commit changed from 39fc2278ce182b46715ffd4fb2071178a4c11092 to 22746f8e5bd59a1f1f53189dc8e8281cc83933c4

Branch pushed to git repo; I updated commit sha1. New commits:

 ​22746f8 `#10017: adjust doctest for 32-bit`

### comment:30 Changed 4 years ago by cremona

• Status changed from needs_work to needs_review

Apologies for the delay -- I pushed to u/cremona/20017 by mistake 4 days ago. It's now where it should be.

### comment:31 Changed 4 years ago by jdemeyer

OK, looks good but let's wait for the patchbot.

### comment:32 follow-up: ↓ 33 Changed 4 years ago by jdemeyer

• Status changed from needs_review to needs_work

Sorry to say this, but the branch doesn't merge cleanly...

### comment:33 in reply to: ↑ 32 Changed 4 years ago by cremona

Sorry to say this, but the branch doesn't merge cleanly...

Not your fault! I am just building the latest beta and will merge it with my branch shortly.

### comment:34 Changed 4 years ago by git

• Commit changed from 22746f8e5bd59a1f1f53189dc8e8281cc83933c4 to 8f1f51d83fc198a31683740001bcffe9a3d30bbc

Branch pushed to git repo; I updated commit sha1. New commits:

 ​f56a7e8 `Merge branch 'develop' into 10017` ​8f1f51d `#10017: fix deprecation warning for Minkowski_embedding`

### comment:35 Changed 4 years ago by cremona

• Status changed from needs_work to needs_review

OK, I merged, fixed the small conflict and a deprecation warning. Let's see if the patchbots like this one better.

### comment:36 Changed 4 years ago by jdemeyer

• Status changed from needs_review to positive_review

### comment:37 Changed 4 years ago by vbraun

• Branch changed from u/cremona/10017 to 8f1f51d83fc198a31683740001bcffe9a3d30bbc
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.