# closest_vector for IntegerLattice is broken

Reported by: Owned by: dimpase blocker sage-9.2 number theory round malb, slelievre, chapoton Samuel Lelièvre, Dima Pasechnik Matthias Koeppe, Samuel Lelièvre N/A 99f0df7 99f0df74fa4b4b1c836b8b850ceeb516ef012d22

In this example we compute `v`, the closest vector to the vector `u` that lies in the integer lattice `L`. Problem: `v` is not in `L`.

```sage: from sage.modules.free_module_integer import IntegerLattice
sage: M = matrix(ZZ, [[20957228, -4966110], [9411844, 19625639]])
sage: L = IntegerLattice(M)
sage: u = vector([-423434678248195, -18882583298608161305227077482])
sage: v = L.closest_vector(u)
sage: v in L
False
sage: MGAP = libgap(M)
sage: vGAP = libgap(v)
sage: libgap.SolutionIntMat(MGAP, vGAP)  # independent test that v is not in L
fail
```

This is a simplified version of Taylor Huang's example posted on sage-devel.

This used to work in Sage 8.4, and was broken in 8.5.beta3 by the change from `round(x, 6)` to `x.n(digits=6)` in #26648.

The sage-devel thread also reports the membership test to be broken:

```from sage.modules.free_module_integer import IntegerLattice
coef = Matrix([-44429982080874270968379672793605458,
98931650854481334735580708522902113])
bMat = Matrix([[20957228, -4966110], [9411844, 19625639]])
L = IntegerLattice(bMat)
coef*bMat in L  # should be true, but returns False
bMatGAP = libgap(bMat)
v = libgap(coef*bMat)[0]
sol = libgap.SolutionIntMat(bMatGAP, v); sol  # prints
# [ -423434671769860, -18882583298608161305225815339 ]
vector(sol.sage())*bMat == vector(coef*bMat)  # sanity check - prints True
```

but this is merely a type "error", in the sense that `coef` is not a vector, but a 1x2 matrix (a different type, and automatic conversion does not happen here). Not sure whether the latter warrants a fix.

### comment:1 Changed 10 months ago by mkoeppe

• Description modified (diff)

### comment:2 Changed 10 months ago by slelievre

• Description modified (diff)

### comment:3 Changed 10 months ago by dimpase

this might be just a bug in membership testing

```from sage.modules.free_module_integer import IntegerLattice
coef = Matrix([-44429982080874270968379672793605458, 98931650854481334735580708522902113])
bMat = Matrix([[20957228, -4966110],[ 9411844, 19625639]])
L = IntegerLattice(bMat)
coef*bMat in L
```

prints `False`. (cf. continutiation of the thread on sage-devel in ticket description)

Last edited 10 months ago by dimpase (previous) (diff)

### comment:4 Changed 10 months ago by gh-asd00012334

Hi, I'm the person who came up with these examples. I can be sure that it is not only about membership testing. I've checked my examples with Smith normal form and indeed some of the closest vectors returned is not in L. So both membership and closest_vector are broken.

Last edited 10 months ago by gh-asd00012334 (previous) (diff)

### comment:5 Changed 10 months ago by dimpase

we confirm that closest vector is broken, nothing to do with membership test bug

```sage: MGAP=libgap(M)
sage: vGAP=libgap(v)
sage: libgap.SolutionIntMat(MGAP,vGAP)
fail
```

### comment:6 Changed 10 months ago by dimpase

• Description modified (diff)
• Priority changed from critical to blocker
• Summary changed from closest_vector for IntegerLattice is broken to closest_vector for IntegerLattice and membership test are broken

### comment:7 Changed 10 months ago by dimpase

• Description modified (diff)

### comment:8 Changed 10 months ago by dimpase

• Description modified (diff)

### comment:9 Changed 10 months ago by dimpase

• Description modified (diff)

### comment:10 Changed 10 months ago by dimpase

closest vector got broken somewhere between Sage 8.4 (works) and 8.5 (does not). Membership test already broken in 8.2.

Last edited 10 months ago by dimpase (previous) (diff)

### comment:11 Changed 10 months ago by dimpase

the heavy lifting is done by fplll/fpylll.

```Python 3.7.7 (default, May 30 2020, 01:27:43)
[GCC 8.3.1 20190518] on linux
>>> from fpylll import CVP, IntegerMatrix
>>> A = IntegerMatrix.from_matrix([[20957228, -4966110], [9411844, 19625639]])
>>> t=(-423434678248195, -18882583298608161305227077482)
>>> v0 = CVP.closest_vector(A, t)
>>> v0
(-423434681577544, -18882583298608161305224294058)
```

this is the same vector as computed by Sage code in the ticket description, oops. So the bug may be narrowed to fp(y)lll.

This vector is different from what Sage outputs, namely

```(-423434671769860, -18882583298608161305225815339)
```

and may be verified to be in L by calling libgap, say, or just directly:

```sage: vector((-423434681577544, -18882583298608161305224294058)) in L
True
```

Note a comment on https://github.com/fplll/fpylll/issues/124 which says that the basis should be LLL-reduced for this to work correctly.

Last edited 10 months ago by dimpase (previous) (diff)

### comment:12 Changed 10 months ago by dimpase

Thus, in view of comment:11 (corrected) the closest_vector bug is a Sage interface to fp(y)lll bug.

rather, Sage has a different implementation --- which seems to be unchanged since 2015 (i.e. much before Sage 8.5...)

Last edited 10 months ago by dimpase (previous) (diff)

### comment:13 Changed 10 months ago by dimpase

• Description modified (diff)

### comment:14 Changed 10 months ago by dimpase

• Description modified (diff)
• Summary changed from closest_vector for IntegerLattice and membership test are broken to closest_vector for IntegerLattice is broken

The containment is "broken", as the 1x2 matrix `coef*bMat` fails to be coersed to a vector. Namely, with the example in the

```sage: vector(coef)*bMat in L
True
sage: vector(coef*bMat) in L
True
sage: coef*bMat in L
False
sage: M=matrix([1,1])
sage: M*bMat in L
False
```

so this is not really a bug, it's just a conversion issue, absence of automatic conversion.

### comment:15 Changed 10 months ago by slelievre

This seems suspicious.

```sage: from sage.modules.free_module_integer import IntegerLattice
sage: a = vector(ZZ, [20957228, -4966110])
sage: b = vector(ZZ, [9411844, 19625639])
sage: L = IntegerLattice([a, b])
sage: R = L.voronoi_relevant_vectors()
sage: S = [-b, -a, b - a, b, a, a - b]
sage: print([r - s for r, s in zip(R, S)])
[(0, 0), (0, 0), (0, -1), (0, 0), (0, 0), (0, 1)]
sage: print(R[-1]); print(a - b)
(11545384, -24591748)
(11545384, -24591749)
```

### comment:16 Changed 10 months ago by dimpase

note that on Sage 8.4 (where closest_vector() still worked) this is

```...
sage: print([r - s for r, s in zip(R, S)])
[(-2133540, 44217388), (11545384, -24591749), (-9411844, -19625639), (2133540, -44217388), (0, 0), (-2133540, 44217388)]
sage: print(R[-1]); print(a - b)
(9411844, 19625639)
(11545384, -24591749)
```

### comment:17 Changed 10 months ago by dimpase

If fact, R computed in comment:15 is the same as the one in comment:16 (with Sage 8.4), but in a different order.

### comment:18 Changed 10 months ago by slelievre

The voronoi cell is slightly off. It should be exactly `P` below.

```sage: V = L.voronoi_cell()
sage: ieqs = [(x^2 + y^2, 2*x, 2*y) for x, y in (-b, -a, a-b, b, a, b-a)]
sage: P = Polyhedron(ieqs=ieqs)
sage: VP = V.plot(fill='cyan', alpha=0.5) + P.plot(fill='magenta', alpha=0.5)
sage: for v in V.vertices():
....:     x, y = v.vector()
....:     VP.show(xmin=x-1, xmax=x+1, ymin=y-1, ymax=y+1, figsize=3)
```
Last edited 10 months ago by slelievre (previous) (diff)

### comment:19 Changed 10 months ago by dimpase

indeed, on Sage 9.2.beta1

```sage: V.vertices()
(A vertex at (-483879947381762629327/57254902852288, 315998472693125499957/28627451426144),
A vertex at (483879947381762629327/57254902852288, -315998472693125499957/28627451426144),
A vertex at (5728193175888805065973/458039243775532, 1390651404779390505975/229019621887766),
A vertex at (-708599510920502622251/229019617181844, 1552002451289776620941/114509808590922),
A vertex at (708599510920502622251/229019617181844, -1552002451289776620941/114509808590922),
A vertex at (-5728193175888805065973/458039243775532, -1390651404779390505975/229019621887766))
```

but on Sage 8.4

```sage: V.vertices()
(A vertex at (5728193175888805065973/458039243775532, 1390651404779390505975/229019621887766),
A vertex at (3871039688862599879323/458039243775532, -2527988039232444116235/229019621887766),
A vertex at (1417199267595526864965/458039243775532, -3104005018306403526499/229019621887766),
A vertex at (-5728193175888805065973/458039243775532, -1390651404779390505975/229019621887766),
A vertex at (-3871039688862599879323/458039243775532, 2527988039232444116235/229019621887766),
A vertex at (-1417199267595526864965/458039243775532, 3104005018306403526499/229019621887766))
```

quite a bit different

### comment:20 Changed 10 months ago by dimpase

I'll try bisecting Sage 8.4 and 8.5. Might take a while...

### comment:21 Changed 10 months ago by slelievre

Could it be #26648 with this change to the function `diamond_cut` in `src/sage/modules/diamond_cutting.py`?

```-                hv = [QQ(round(elmt, 6)) for elmt in hv]
+                hv = [QQ(elmt.n(digits=6)) for elmt in hv]
```

### comment:22 Changed 10 months ago by dimpase

good catch - this change is certainly bad. cf.

```sage: r=5555.55555555555555555
sage: round(r,6)
5555.555556
sage: r.n(digits=6)
5555.56
```

why it's always 6 that suffices is another question, though

### comment:23 Changed 10 months ago by dimpase

yes, this fixes this error:

• ## src/sage/modules/diamond_cutting.py

 a def diamond_cut(V, GM, C, verbose=False): (A vertex at (2), A vertex at (0)) """ # coerce to floats from sage.misc.functional import round GM = GM.n() C = float(C) if verbose: def diamond_cut(V, GM, C, verbose=False): cut_count += 1 if verbose: print("\n%d) Cut using normal vector %s" % (cut_count, hv)) hv = [QQ(elmt.n(digits=6)) for elmt in hv] hv = [QQ(round(elmt,6)) for elmt in hv] # hv = [QQ(elmt.n(digits=6)) for elmt in hv] inequalities.append(plane_inequality(hv)) if verbose:

### comment:24 Changed 10 months ago by dimpase

• Authors set to Samuel Lelièvre, Dima Pasechnik
• Branch set to u/dimpase/modules/corrround
• Commit set to d734de1617b0ab38a60d7bc2e27f8e77107bede4
• Status changed from new to needs_review

I wonder how many more of these `round()` calls were butchered this way, or this was an one-off event.

New commits:

 ​d734de1 `the correct precision of 6 digits after dot`

### comment:25 Changed 10 months ago by dimpase

The line

```hv = [QQ(round(elmt, 6)) for elmt in hv]
```

still bothers me a lot. Why `6` there? Is it just waiting for an example with bigger numbers to break it again?

Last edited 10 months ago by dimpase (previous) (diff)

### comment:26 Changed 10 months ago by git

• Commit changed from d734de1617b0ab38a60d7bc2e27f8e77107bede4 to 1faa6aceeb6799ea31517c9b84b3e4164153e78d

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

 ​1faa6ac `get rid of conversion to floats and round()`

### comment:27 Changed 10 months ago by dimpase

this should be robust now

### comment:28 follow-up: ↓ 29 Changed 10 months ago by slelievre

• Description modified (diff)

Could it be that the conversion to `QQ` was there for inexact vectors, say vectors over `RDF` or `RR`?

Should we add the example in the ticket description as a doctest?

### comment:29 in reply to: ↑ 28 ; follow-up: ↓ 30 Changed 10 months ago by dimpase

Could it be that the conversion to `QQ` was there for inexact vectors, say vectors over `RDF` or `RR`?

for IntegerLattice this looks irrelevant.

probably it was a dodgy old workaround, from pre-PPL days.

Should we add the example in the ticket description as a doctest?

OK, I can do this.

### comment:30 in reply to: ↑ 29 ; follow-up: ↓ 32 Changed 10 months ago by slelievre

Could it be that the conversion to `QQ` was there for inexact vectors, say vectors over `RDF` or `RR`?

for IntegerLattice this looks irrelevant.

Even though we are working with an integer lattice, doesn't this offer to find the closest lattice point to a non-integer vector?

### comment:31 Changed 10 months ago by git

• Commit changed from 1faa6aceeb6799ea31517c9b84b3e4164153e78d to 685ac9bdb3b9451c97707add3118d63b1006f366

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

 ​685ac9b `adding test from trac #29866`

### comment:32 in reply to: ↑ 30 Changed 10 months ago by dimpase

Could it be that the conversion to `QQ` was there for inexact vectors, say vectors over `RDF` or `RR`?

for IntegerLattice this looks irrelevant.

Even though we are working with an integer lattice, doesn't this offer to find the closest lattice point to a non-integer vector?

the code we're changing only deals with the lattice, not with the vector we are finding a closest lattice point for.

### comment:33 Changed 10 months ago by slelievre

Right! One last thing, QQ is now imported but not used, can you remove the import?

### comment:34 Changed 10 months ago by git

• Commit changed from 685ac9bdb3b9451c97707add3118d63b1006f366 to 99f0df74fa4b4b1c836b8b850ceeb516ef012d22

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

 ​99f0df7 `remove unused import`

done

### comment:37 Changed 10 months ago by mkoeppe

• Reviewers set to Matthias Koeppe
• Status changed from needs_review to positive_review

### comment:38 Changed 10 months ago by slelievre

• Reviewers changed from Matthias Koeppe to Matthias Koeppe, Samuel Lelièvre

Thanks!

### comment:39 Changed 10 months ago by vbraun

• Branch changed from u/dimpase/modules/corrround to 99f0df74fa4b4b1c836b8b850ceeb516ef012d22
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.