Repository: fplll/fpylll (language: python)

Issue:
Babai failure and precision issues
The attached file (which is really a Python file, but can't be attached as one) gives an example of Babai's algorithm in fpylll failing when it shouldn't due presumably to precision issues.

[issue.txt](https://github.com/fplll/fpylll/files/8674033/issue.txt)


Discussion / hints:
To make it easier: 

```python
import numpy as np
from fpylll import FPLLL, LLL, GSO, IntegerMatrix
from math import ceil, sqrt

# Global set-up.
FPLLL.set_random_seed(1337)
n = 10

# Just to emphasize that it doesn't help.
FPLLL.float_type="mpfr"
FPLLL.set_precision(1000)

# Set up basis and Gram-Schmidt basis.
B = IntegerMatrix(n, n + 1)
B.randomize("intrel", bits=100)

# Set up target.
x = [1, 0, 0, 1, 0, 0, 1, 0, 0, 0]
v_opt = B.multiply_left(x)
s = v_opt[0] # s = <a, x>, where a is the vector of knapsack values.
t = [s] + (n * [0])

# Reduce basis. 
LLL.reduction(B, delta=0.75)
M = GSO.Mat(B)
M.update_gso()

# Get a coefficient vector of a close vector to t from Babai.
y = M.babai(t)
v = B.multiply_left(y)

# These are not equal.
print("t[0]:", t[0])
print("v[0]:", v[0])
print("v[0] = t[0]?", v[0] == t[0])

# Upper bound on Babai approx. factor using a 0.75-LLL-reduced basis.
gamma = ceil(sqrt((n + 1) * (2 ** n)))

# Using NumPy for linear algebra.
# Type numbers in NumPy as Python objects to avoid precision issues.
np_array_obj = lambda a: np.array(a, "object")

# Babai does not work correctly.
assert(np.linalg.norm(np_array_obj(v) - np_array_obj(t)) <= \
       gamma * np.linalg.norm(np_array_obj(v_opt) - np_array_obj(t)))
```
I am hit by the same problem. This other Sagemath script illustrates the problem.

```python
import random
import fpylll

q = random_prime(n=2**256)
bound = 10
n = 10         # variables
m = 5          # linear equations

# build ISIS instance
K = GF(q)
A = matrix.random(GF(q), n, m)
x = vector(K, [random.randrange(bound) for _ in range(n)])
y = x * A

# assemble CVP instance
mu = 1
M = matrix.block(ZZ, [[mu*A, 1], [mu*q, 0]])
t = vector(ZZ, list(mu * vector(ZZ, y)) + [0 for _ in range(n)])   # CVP target
s = vector(ZZ, list(mu * vector(ZZ, y)) + list(x))                 # lattice vector we hope to find

# check that s actually belongs to the lattice
u = vector(ZZ, M.solve_left(s))
assert u*M == s

# check that the distance between s and t is small
dist = (s - t).norm().n()
assert dist < sqrt(n) * bound

# solve CVP exactly using enumeration
G = fpylll.IntegerMatrix.from_matrix(M)
fpylll.LLL.reduction(G)
cvp = fpylll.CVP.closest_vector(G, t)             # finds s in principle
assert cvp == tuple(s)


# now try again using Babai's nearest plane
H = fpylll.IntegerMatrix.from_matrix(M)
_ = fpylll.LLL.reduction(H)
S = fpylll.GSO.Mat(H)
_ = S.update_gso()
w = S.babai(t)
cvp_bis = vector(ZZ, H.multiply_left(w))          # ERRROR ! Very far away from target


# check that CVP approximation ratio is less than sqrt(4/3)**(n+m)
actual_dist = (cvp_bis - t).norm().n()            # ERROR ! astoundingly large
assert actual_dist < 1.1547**(n+m) * dist
```

Fix patch:
diff --git a/src/fpylll/algorithms/babai.py b/src/fpylll/algorithms/babai.py
new file mode 100644
index 00000000..b7942d86
--- /dev/null
+++ b/src/fpylll/algorithms/babai.py
@@ -0,0 +1,67 @@
+# -*- coding: utf-8 -*-
+"""
+Babai's Nearest Plane algorithm
+
+
+..  moduleauthor:: Martin R.  Albrecht <martinralbrecht+fpylll@googlemail.com>
+
+"""
+
+from fpylll import IntegerMatrix, LLL
+from math import ceil
+
+
+def babai(B, t):
+    """
+    Run Babai's Nearest Plane algorithm by running LLL.
+
+    :param B: Input lattice basis.
+    :param target: Target point (∈ ZZ^n)
+    :returns coordinates of the solution vector:
+
+    This implementation is more numerically stable compared to the one offered by `MatGSO.babai()`.
+    On the other hand, this implementation will only accept tatgets with Integer coefficients.
+
+    EXAMPLE ::
+
+       >>> from fpylll import *
+       >>> n = 10
+       >>> B = IntegerMatrix(n, n + 1)
+       >>> B.randomize("intrel", bits=100)
+       >>> v_opt = B.multiply_left([1,0,1,0,1,1,0,0,1,1])
+       >>> s = v_opt[0] # s = <a, x>, where a is vector of knapsack values.
+       >>> t = [s] + (n * [0])
+       >>> v = CVP.babai(B, t)
+       >>> v[0] == t[0]
+       True
+       >>> v[1:]
+       (1, 0, 1, 0, 1, 1, 0, 0, 1, 1)
+       >>> _ = LLL.reduction(B)
+       >>> v == CVP.closest_vector(B, t)
+       True
+    """
+    A = IntegerMatrix(B.nrows + 1, B.ncols + 1)
+    for i in range(B.nrows):
+        for j in range(B.ncols):
+            A[i, j] = B[i, j]
+
+    # make sure the input is LLL reduced before reading the norm of the last vector
+    LLL.reduction(A)
+    # zero vector at the end
+    A.swap_rows(0, B.nrows)
+
+    for j in range(B.ncols):
+        A[-1, j] = t[j]
+    A[-1, -1] = ceil(A[-2].norm())
+
+    LLL.reduction(A)  # now call LLL to run Babai
+
+    v = [0] * len(t)
+    if A[-1, -1] > 0:
+        for i in range(len(t)):
+            v[i] = t[i] - A[-1][i]
+    else:
+        for i in range(len(t)):
+            v[i] = t[i] + A[-1][i]
+
+    return tuple(v)
diff --git a/src/fpylll/algorithms/bkz2.py b/src/fpylll/algorithms/bkz2.py
index 3228bb0c..ccfcb724 100644
--- a/src/fpylll/algorithms/bkz2.py
+++ b/src/fpylll/algorithms/bkz2.py
@@ -1,4 +1,10 @@
 # -*- coding: utf-8 -*-
+"""
+BKZ 2.0 algorithm in Python.
+
+..  moduleauthor:: Martin R.  Albrecht <martinralbrecht+fpylll@googlemail.com>
+
+"""
 
 from fpylll import BKZ, Enumeration, EnumerationError
 from fpylll.algorithms.bkz import BKZReduction as BKZBase
diff --git a/src/fpylll/fplll/decl.pxd b/src/fpylll/fplll/decl.pxd
index f199d984..b3be8d1c 100644
--- a/src/fpylll/fplll/decl.pxd
+++ b/src/fpylll/fplll/decl.pxd
@@ -460,3 +460,7 @@ ELSE:
             vector[FP_NR[d_t]] d
             vector[FP_NR[dpe_t]] dpe
             vector[FP_NR[mpfr_t]] mpfr
+
+ctypedef struct vector_z_nr_t:
+    vector[Z_NR[mpz_t]] mpz
+    vector[Z_NR[long]] long
diff --git a/src/fpylll/fplll/fplll.pxd b/src/fpylll/fplll/fplll.pxd
index 2635c250..4b5207ec 100644
--- a/src/fpylll/fplll/fplll.pxd
+++ b/src/fpylll/fplll/fplll.pxd
@@ -124,6 +124,7 @@ cdef extern from "fplll/nr/nr.h" namespace "fplll":
         double get_d() nogil
         inline void operator=(const FP_NR[T]& a) nogil
         inline void operator=(double a) nogil
+        inline void operator=(const char *s) nogil
 
         @staticmethod
         unsigned int get_prec() nogil
@@ -464,6 +465,10 @@ cdef extern from "fplll/gso.h" namespace "fplll":
         FT get_log_det(int start_row, int stop_row) nogil
         FT get_slide_potential(int start_row, int stop_row, int block_size) nogil
 
+        void to_canonical(vector[FT] &w, const vector[FT] &v, long start) nogil
+        void from_canonical(vector[FT] &v, const vector[FT] &w, long start, long dimension) nogil
+        int babai(vector[ZT] w, vector[FT] v, int start, int dimension, bool gsa) nogil
+
         const int enable_int_gram
         const int enable_row_expo
         const int enable_transform
@@ -565,6 +570,8 @@ cdef extern from "fplll/gso_interface.h" namespace "fplll":
         FT get_log_det(int start_row, int end_row) nogil
         FT get_slide_potential(int start_row, int end_row, int block_size) nogil
 
+        int babai(vector[ZT] w, vector[FT] v, int start, int dimension) nogil
+
         const bool enable_int_gram
         const bool enable_row_expo
         const bool enable_transform
@@ -989,5 +996,3 @@ cdef extern from "fplll/fplll.h" namespace "fplll":
     int hkz_reduction(ZZ_mat[mpz_t] b) nogil
 
     const char* get_red_status_str(int status) nogil
-
-
diff --git a/src/fpylll/fplll/gso.pyx b/src/fpylll/fplll/gso.pyx
index 57d6e812..ca606e55 100644
--- a/src/fpylll/fplll/gso.pyx
+++ b/src/fpylll/fplll/gso.pyx
@@ -10,17 +10,19 @@ A ``MatGSO`` object stores the following information:
 
     - the Gram-Schmidt coefficients `μ_{i,j} = `⟨b_i, b^*_j⟩ / ||b^*_j||^2` for `i>j`, and
 
-    - the coefficients `r_{i,j} = ⟨b_i, b^*_j⟩` for `i>j`
+    - the coefficients `r_{i,j} = ⟨b_i, b^*_j⟩` for `i≥j`
 
 It holds that: `B = R × Q = (μ × D) × (D^{-1} × B^*)` where `Q` is orthonormal and `R` is lower
 triangular.
+
 """
 
 include "fpylll/config.pxi"
 
 from cysignals.signals cimport sig_on, sig_off
 
-from .decl cimport mat_gso_mpz_d, mat_gso_mpz_ld, mat_gso_mpz_dpe, mat_gso_mpz_mpfr, fp_nr_t, zz_mat_core_t, z_nr_t
+from .decl cimport mat_gso_mpz_d, mat_gso_mpz_ld, mat_gso_mpz_dpe, mat_gso_mpz_mpfr
+from .decl cimport vector_fp_nr_t, vector_z_nr_t, z_nr_t, fp_nr_t, zz_mat_core_t
 from .decl cimport mat_gso_long_d, mat_gso_long_ld, mat_gso_long_dpe, mat_gso_long_mpfr
 from .decl cimport d_t
 from .decl cimport mat_gso_gso_t, mat_gso_gram_t
@@ -37,6 +39,7 @@ from .fplll cimport get_current_slope
 from fpylll.gmp.mpz cimport mpz_t
 from fpylll.mpfr.mpfr cimport mpfr_t
 from fpylll.util cimport preprocess_indices, check_float_type
+from fpylll.io cimport vector_fp_nr_barf, vector_fp_nr_slurp, vector_z_nr_slurp
 from fpylll.io cimport mpz_get_python
 from .integer_matrix cimport IntegerMatrix
 
@@ -2032,12 +2035,12 @@ cdef class MatGSO:
         raise RuntimeError("MatGSO object '%s' has no core."%self)
 
 
-    def from_canonical(self, v, int start=0, int dimension=-1):
-        """Given a vector `v` wrt the canonical basis `\mathbb{Z}^n` return a vector wrt the
+    def from_canonical(self, w, int start=0, int dimension=-1):
+        """Given a vector `w` wrt the canonical basis `\ZZ^n` return a vector `v` wrt the
         Gram-Schmidt basis `B^*`
 
         :param v: a tuple-like object of dimension ``M.B.ncols``
-        :param start: only consider subbasis starting at ``start```
+        :param start: only consider subbasis starting at ``start``
         :param dimension: only consider ``dimension`` vectors or all if ``-1``
 
         :returns: a tuple of dimension ``dimension``` or ``M.d``` when ``dimension`` is ``None``
@@ -2058,67 +2061,244 @@ cdef class MatGSO:
             True
 
         """
-        cdef Py_ssize_t i, j, d
-
         if self._alg != mat_gso_gso_t:
             raise TypeError("This function is only defined for GSO objects over a basis")
 
-        if dimension == -1:
-            d = self.d - start
-        else:
-            d = dimension
+        cdef vector_fp_nr_t cv, cw
+        cdef fp_nr_t tmp
+
+        if self._type == mat_gso_mpz_d:
+            vector_fp_nr_barf(cw, w, FT_DOUBLE)
+            sig_on()
+            (<MatGSO_c[Z_NR[mpz_t],FP_NR[d_t]]*>self._core.mpz_d).from_canonical(cv.d, cw.d, start, dimension)
+            sig_off()
+            return vector_fp_nr_slurp(cv, FT_DOUBLE)
+        elif self._type == mat_gso_mpz_ld:
+            vector_fp_nr_barf(cw, w, FT_LONG_DOUBLE)
+            sig_on()
+            (<MatGSO_c[Z_NR[mpz_t],FP_NR[ld_t]]*>self._core.mpz_ld).from_canonical(cv.ld, cw.ld, start, dimension)
+            sig_off()
+            return vector_fp_nr_slurp(cv, FT_LONG_DOUBLE)
+        # elif self._type == mat_gso_mpz_dpe:
+        #     vector_fp_nr_barf(cw, w, FT_DPE)
+        #     sig_on()
+        #     (<MatGSO_c[Z_NR[mpz_t],FP_NR[dpe_t]]*>self._core.mpz_dpe).from_canonical(cv.dpe, cw.dpe, start, dimension)
+        #     sig_off()
+        #     return vector_fp_nr_slurp(cv, FT_DPE)
+        elif self._type == mat_gso_mpz_mpfr:
+            vector_fp_nr_barf(cw, w, FT_MPFR)
+            sig_on()
+            (<MatGSO_c[Z_NR[mpz_t],FP_NR[mpfr_t]]*>self._core.mpz_mpfr).from_canonical(cv.mpfr, cw.mpfr, start, dimension)
+            sig_off()
+            return vector_fp_nr_slurp(cv, FT_MPFR)
+        elif self._type == mat_gso_long_d:
+            vector_fp_nr_barf(cw, w, FT_DOUBLE)
+            sig_on()
+            (<MatGSO_c[Z_NR[long],FP_NR[d_t]]*>self._core.long_d).from_canonical(cv.d, cw.d, start, dimension)
+            sig_off()
+            return vector_fp_nr_slurp(cv, FT_DOUBLE)
+        elif self._type == mat_gso_long_ld:
+            vector_fp_nr_barf(cw, w, FT_LONG_DOUBLE)
+            sig_on()
+            (<MatGSO_c[Z_NR[long],FP_NR[ld_t]]*>self._core.long_ld).from_canonical(cv.ld, cw.ld, start, dimension)
+            sig_off()
+            return vector_fp_nr_slurp(cv, FT_LONG_DOUBLE)
+        # elif self._type == mat_gso_long_dpe:
+        #     vector_fp_nr_barf(cw, w, FT_DPE)
+        #     sig_on()
+        #     (<MatGSO_c[Z_NR[long],FP_NR[dpe_t]]*>self._core.long_dpe).from_canonical(cv.dpe, cw.dpe, start, dimension)
+        #     sig_off()
+        #     return vector_fp_nr_slurp(cv, FT_DPE)
+        # elif self._type == mat_gso_long_mpfr:
+        #     vector_fp_nr_barf(cw, w, FT_MPFR)
+        #     sig_on()
+        #     (<MatGSO_c[Z_NR[long],FP_NR[mpfr_t]]*>self._core.long_mpfr).from_canonical(cv.mpfr, cw.mpfr, start, dimension)
+        #     sig_off()
+        #     return vector_fp_nr_slurp(cv, FT_MPFR)
+        IF HAVE_QD:
+            if self._type == mat_gso_mpz_dd:
+                vector_fp_nr_barf(cw, w, FT_DD)
+                sig_on()
+                (<MatGSO_c[Z_NR[mpz_t],FP_NR[dd_t]]*>self._core.mpz_dd).from_canonical(cv.dd, cw.dd, start, dimension)
+                sig_off()
+                return vector_fp_nr_slurp(cv, FT_DD)
+            elif self._type == mat_gso_mpz_qd:
+                vector_fp_nr_barf(cw, w, FT_QD)
+                sig_on()
+                (<MatGSO_c[Z_NR[mpz_t],FP_NR[qd_t]]*>self._core.mpz_ld).from_canonical(cv.qd, cw.qd, start, dimension)
+                sig_off()
+                return vector_fp_nr_slurp(cv, FT_QD)
+            elif self._type == mat_gso_long_dd:
+                vector_fp_nr_barf(cw, w, FT_DD)
+                sig_on()
+                (<MatGSO_c[Z_NR[long],FP_NR[dd_t]]*>self._core.long_dd).from_canonical(cv.dd, cw.dd, start, dimension)
+                sig_off()
+                return vector_fp_nr_slurp(cv, FT_DD)
+            # elif self._type == mat_gso_long_qd:
+            #     vector_fp_nr_barf(cw, w, FT_QD)
+            #     sig_on()
+            #     (<MatGSO_c[Z_NR[long],FP_NR[qd_t]]*>self._core.long_ld).from_canonical(cv.qd, cw.qd, start, dimension)
+            #     sig_off()
+            #     return vector_fp_nr_slurp(cw, FT_QD)
+
+        raise NotImplementedError
+
+
+    # def _from_canonical_old(self, w, int start=0, int dimension=-1):
+    #     cdef Py_ssize_t i, j, d
 
-        cdef list ret = [0]*(start+d)
-        for i in range(start+d):
-            for j in range(self.B.ncols):
-                ret[i] += self.B[i, j] * v[j]
+    #     if self._alg != mat_gso_gso_t:
+    #         raise TypeError("This function is only defined for GSO objects over a basis")
 
-            for j in range(i):
-                ret[i] -= self.get_mu(i, j) * ret[j]
+    #     if dimension == -1:
+    #         d = self.d - start
+    #     else:
+    #         d = dimension
 
-        # we drop the first ``start`` entries anyway, so no need to update
-        for i in range(d):
-            ret[start+i] /= self.get_r(start+i, start+i)
+    #     cdef list ret = [0]*(start+d)
+    #     for i in range(start+d):
+    #         for j in range(self.B.ncols):
+    #             ret[i] += self.B[i, j] * w[j]
 
-        return tuple(ret)[start:]
+    #         for j in range(i):
+    #             ret[i] -= self.get_mu(i, j) * ret[j]
 
+    #     # we drop the first ``start`` entries anyway, so no need to update
+    #     for i in range(d):
+    #         ret[start+i] /= self.get_r(start+i, start+i)
+
+    #     return tuple(ret)[start:]
 
     def to_canonical(self, v, int start=0):
         """
-        Given a vector `v` wrt the Gram-Schmidt basis `B^*` return a vector wrt the canonical basis
-        `\mathbb{Z}^n`
+        Given a vector `v` wrt the Gram-Schmidt basis `B^*` return a vector `w` wrt the
+        canonical basis `ZZ^n`, i.e. solve `w = v⋅B^*`.
 
         :param v: a tuple-like object of dimension ``M.d``
-        :param start: only consider subbasis starting at ``start```
+        :param start: only consider subbasis starting at ``start``
 
         :returns: a tuple of dimension ``M.B.ncols``
-
         """
         if self._alg != mat_gso_gso_t:
             raise TypeError("This function is only defined for GSO objects over a basis")
 
-        cdef list vv = list(v)
-        cdef Py_ssize_t i, j
-        cdef Py_ssize_t d = min(len(vv), self.d-start)
-        for i in range(d)[::-1]:
-            for j in range(i+1, d):
-                vv[i] -= self.get_mu(start+j, start+i) * vv[j]
+        cdef vector_fp_nr_t cv, cw
+        cdef fp_nr_t tmp
 
-        ret = [0]*self.B.ncols
-        for i in range(d):
-            for j in range(self.B.ncols):
-                ret[j] += vv[i] * self.B[start+i,j]
+        if self._type == mat_gso_mpz_d:
+            vector_fp_nr_barf(cv, v, FT_DOUBLE)
+            sig_on()
+            (<MatGSO_c[Z_NR[mpz_t],FP_NR[d_t]]*>self._core.mpz_d).to_canonical(cw.d, cv.d, start)
+            sig_off()
+            return vector_fp_nr_slurp(cw, FT_DOUBLE)
+        elif self._type == mat_gso_mpz_ld:
+            vector_fp_nr_barf(cv, v, FT_LONG_DOUBLE)
+            sig_on()
+            (<MatGSO_c[Z_NR[mpz_t],FP_NR[ld_t]]*>self._core.mpz_ld).to_canonical(cw.ld, cv.ld, start)
+            sig_off()
+            return vector_fp_nr_slurp(cw, FT_LONG_DOUBLE)
+        # # https://github.com/fplll/fplll/issues/493
+        # elif self._type == mat_gso_mpz_dpe:
+        #     vector_fp_nr_barf(cv, v, FT_DPE)
+        #     sig_on()
+        #     (<MatGSO_c[Z_NR[mpz_t],FP_NR[dpe_t]]*>self._core.mpz_dpe).to_canonical(cw.dpe, cv.dpe, start)
+        #     sig_off()
+        #     return vector_fp_nr_slurp(cw, FT_DPE)
+        elif self._type == mat_gso_mpz_mpfr:
+            vector_fp_nr_barf(cv, v, FT_MPFR)
+            sig_on()
+            (<MatGSO_c[Z_NR[mpz_t],FP_NR[mpfr_t]]*>self._core.mpz_mpfr).to_canonical(cw.mpfr, cv.mpfr, start)
+            sig_off()
+            return vector_fp_nr_slurp(cw, FT_MPFR)
+        elif self._type == mat_gso_long_d:
+            vector_fp_nr_barf(cv, v, FT_DOUBLE)
+            sig_on()
+            (<MatGSO_c[Z_NR[long],FP_NR[d_t]]*>self._core.long_d).to_canonical(cw.d, cv.d, start)
+            sig_off()
+            return vector_fp_nr_slurp(cw, FT_DOUBLE)
+        elif self._type == mat_gso_long_ld:
+            vector_fp_nr_barf(cv, v, FT_LONG_DOUBLE)
+            sig_on()
+            (<MatGSO_c[Z_NR[long],FP_NR[ld_t]]*>self._core.long_ld).to_canonical(cw.ld, cv.ld, start)
+            sig_off()
+            return vector_fp_nr_slurp(cw, FT_LONG_DOUBLE)
+        # # https://github.com/fplll/fplll/issues/493
+        # elif self._type == mat_gso_long_dpe:
+        #     vector_fp_nr_barf(cv, v, FT_DPE)
+        #     sig_on()
+        #     (<MatGSO_c[Z_NR[long],FP_NR[dpe_t]]*>self._core.long_dpe).to_canonical(cw.dpe, cv.dpe, start)
+        #     sig_off()
+        #     return vector_fp_nr_slurp(cw, FT_DPE)
+        # elif self._type == mat_gso_long_mpfr:
+        #     vector_fp_nr_barf(cv, v, FT_MPFR)
+        #     sig_on()
+        #     (<MatGSO_c[Z_NR[long],FP_NR[mpfr_t]]*>self._core.long_mpfr).to_canonical(cw.mpfr, cv.mpfr, start)
+        #     sig_off()
+        #     return vector_fp_nr_slurp(cw, FT_MPFR)
+        IF HAVE_QD:
+            if self._type == mat_gso_mpz_dd:
+                vector_fp_nr_barf(cv, v, FT_DD)
+                sig_on()
+                (<MatGSO_c[Z_NR[mpz_t],FP_NR[dd_t]]*>self._core.mpz_dd).to_canonical(cw.dd, cv.dd, start)
+                sig_off()
+                return vector_fp_nr_slurp(cw, FT_DD)
+            # # https://github.com/fplll/fplll/issues/493
+            # elif self._type == mat_gso_mpz_qd:
+            #     vector_fp_nr_barf(cv, v, FT_QD)
+            #     sig_on()
+            #     (<MatGSO_c[Z_NR[mpz_t],FP_NR[qd_t]]*>self._core.mpz_ld).to_canonical(cw.qd, cv.qd, start)
+            #     sig_off()
+            #     return vector_fp_nr_slurp(cw, FT_QD)
+            elif self._type == mat_gso_long_dd:
+                vector_fp_nr_barf(cv, v, FT_DD)
+                sig_on()
+                (<MatGSO_c[Z_NR[long],FP_NR[dd_t]]*>self._core.long_dd).to_canonical(cw.dd, cv.dd, start)
+                sig_off()
+                return vector_fp_nr_slurp(cw, FT_DD)
+            # # https://github.com/fplll/fplll/issues/493
+            # elif self._type == mat_gso_long_qd:
+            #     vector_fp_nr_barf(cv, v, FT_QD)
+            #     sig_on()
+            #     (<MatGSO_c[Z_NR[long],FP_NR[qd_t]]*>self._core.long_ld).to_canonical(cw.qd, cv.qd, start)
+            #     sig_off()
+            #     return vector_fp_nr_slurp(cw, FT_QD)
+
+        raise NotImplementedError
 
-        return tuple(ret)
+    # def _to_canonical_old(self, v, int start=0):
 
+    #     # We have some `v` s.t. `w = v ⋅ B^*` and we want to recover `w`. We do not store `B^*`, but
+    #     # we store `(B, μ)` s.t. `μ ⋅ B^* = B`:
+    #     # - `w = v ⋅ B^*`
+    #     # - `w = v ⋅ μ^{-1} ⋅ B`
+    #     # - let `x` s.t. `x⋅μ = v`
+    #     # - `w = x ⋅ B`
+
+    #     # 1. triangular system solving
+
+    #     cdef list x = list(v)
+    #     cdef Py_ssize_t i, j
+    #     cdef Py_ssize_t d = min(len(x), self.d-start)
+    #     for i in range(d)[::-1]:
+    #         for j in range(i+1, d):
+    #             x[i] -= self.get_mu(start+j, start+i) * x[j]
+
+    #     # 2. multiply by `B`
+
+    #     w = [0]*self.B.ncols
+    #     for i in range(d):
+    #         for j in range(self.B.ncols):
+    #             w[j] += x[i] * self.B[start+i,j]
+
+    #     return tuple(w)
 
     def babai(self, v, int start=0, int dimension=-1, gso=False):
         """
         Return vector `w` s.t. `‖w⋅B - v‖` is small using Babai's nearest plane algorithm.
 
         :param v: a tuple-like object
-        :param start: only consider subbasis starting at ``start```
-        :param dimension: only consider ``dimension`` vectors or all if ``-1```
+        :param start: only consider subbasis starting at ``start``
+        :param dimension: only consider ``dimension`` vectors or all if ``-1``
         :param gso: if ``True`` vector is represented wrt to the Gram-Schmidt basis, otherwise
             canonical basis is assumed.
 
@@ -2143,19 +2323,214 @@ cdef class MatGSO:
             >>> vector_norm(A.multiply_left(w), (100,)*50)
             58851
 
+        We compute a more interesting example and solve a simple Knapsack::
+
+            >>> from fpylll import *
+            >>> _ = FPLLL.set_precision(500)
+            >>> n = 10
+            >>> B = IntegerMatrix(n, n + 1)
+            >>> B.randomize("intrel", bits=500)
+            >>> v_opt = B.multiply_left([1,0,1,0,1,1,0,0,1,1])
+            >>> s = v_opt[0] # s = <a, x>, where a is vector of knapsack values.
+            >>> t = [s] + (n * [0])
+
+            >>> _ = LLL.reduction(B)
+            >>> M = GSO.Mat(B, update=True, float_type="mpfr")
+
+            >>> y = M.babai(t)
+            >>> v = B.multiply_left(y)
+            >>> t[0] == v[0]
+            True
+            >>> v_ = CVP.babai(B, t)
+            >>> v_ == v
+            True
+
+        .. note :: A separate implementation is available at `CVP.babai()`. That implementation is
+           more numerically stable than this one (by repeatedly calling the Nearest Plane algorithm
+           at a given precision to improve the solution). On the other hand, this implementation
+           supports floating point target vectors and `CVP.babai()` does not.
+
         """
         if dimension == -1:
             dimension = self.d - start
+
+        cdef vector_fp_nr_t cv
+        cdef vector_z_nr_t cw
+
+        if gso is False and self._alg != mat_gso_gso_t:
+            raise TypeError("Can only convert to GSO representation with a basis.")
+
         if not gso:
-            v = self.from_canonical(v, start, dimension)
-
-        cdef Py_ssize_t i, j
-        cdef list vv = list(v)
-        for i in range(dimension)[::-1]:
-            vv[i] = int(round(vv[i]))
-            for j in range(i):
-                vv[j] -= self.get_mu(start+i, start+j) * vv[i]
-        return tuple(vv)
+            if self._type == mat_gso_mpz_d:
+                vector_fp_nr_barf(cv, v, FT_DOUBLE)
+                sig_on()
+                (<MatGSO_c[Z_NR[mpz_t],FP_NR[d_t]]*>self._core.mpz_d).babai(cw.mpz, cv.d, start, dimension, gso)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_MPZ)
+            elif self._type == mat_gso_mpz_ld:
+                vector_fp_nr_barf(cv, v, FT_LONG_DOUBLE)
+                sig_on()
+                (<MatGSO_c[Z_NR[mpz_t],FP_NR[ld_t]]*>self._core.mpz_ld).babai(cw.mpz, cv.ld, start, dimension, gso)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_MPZ)
+            # # https://github.com/fplll/fplll/issues/493
+            # elif self._type == mat_gso_mpz_dpe:
+            #     vector_fp_nr_barf(cv, v, FT_DPE)
+            #     sig_on()
+            #     (<MatGSO_c[Z_NR[mpz_t],FP_NR[dpe_t]]*>self._core.mpz_dpe).babai(cw.mpz, cv.dpe, start, dimension, gso)
+            #     sig_off()
+            #     return vector_z_nr_slurp(cw, ZT_MPZ)
+            elif self._type == mat_gso_mpz_mpfr:
+                vector_fp_nr_barf(cv, v, FT_MPFR)
+                sig_on()
+                (<MatGSO_c[Z_NR[mpz_t],FP_NR[mpfr_t]]*>self._core.mpz_mpfr).babai(cw.mpz, cv.mpfr, start, dimension, gso)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_MPZ)
+            elif self._type == mat_gso_long_d:
+                vector_fp_nr_barf(cv, v, FT_DOUBLE)
+                sig_on()
+                (<MatGSO_c[Z_NR[long],FP_NR[d_t]]*>self._core.long_d).babai(cw.long, cv.d, start, dimension, gso)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_LONG)
+            elif self._type == mat_gso_long_ld:
+                vector_fp_nr_barf(cv, v, FT_LONG_DOUBLE)
+                sig_on()
+                (<MatGSO_c[Z_NR[long],FP_NR[ld_t]]*>self._core.long_ld).babai(cw.long, cv.ld, start, dimension, gso)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_LONG)
+            # # https://github.com/fplll/fplll/issues/493
+            # elif self._type == mat_gso_long_dpe:
+            #     vector_fp_nr_barf(cv, v, FT_DPE)
+            #     sig_on()
+            #     (<MatGSO_c[Z_NR[long],FP_NR[dpe_t]]*>self._core.long_dpe).babai(cw.long, cv.dpe, start, dimension, gso)
+            #     sig_off()
+            #     return vector_z_nr_slurp(cw, ZT_LONG)
+            # elif self._type == mat_gso_long_mpfr:
+            #     vector_fp_nr_barf(cv, v, FT_MPFR)
+            #     sig_on()
+            #     (<MatGSO_c[Z_NR[long],FP_NR[mpfr_t]]*>self._core.long_mpfr).babai(cw.long, cv.mpfr, start, dimension, gso)
+            #     sig_off()
+            #     return vector_z_nr_slurp(cw, ZT_LONG)
+            IF HAVE_QD:
+                if self._type == mat_gso_mpz_dd:
+                    vector_fp_nr_barf(cv, v, FT_DD)
+                    sig_on()
+                    (<MatGSO_c[Z_NR[mpz_t],FP_NR[dd_t]]*>self._core.mpz_dd).babai(cw.mpz, cv.dd, start, dimension, gso)
+                    sig_off()
+                    return vector_z_nr_slurp(cw, ZT_MPZ)
+                # # https://github.com/fplll/fplll/issues/493
+                # elif self._type == mat_gso_mpz_qd:
+                #     vector_fp_nr_barf(cv, v, FT_QD)
+                #     sig_on()
+                #     (<MatGSO_c[Z_NR[mpz_t],FP_NR[qd_t]]*>self._core.mpz_qd).babai(cw.mpz, cv.qd, start, dimension, gso)
+                #     sig_off()
+                #     return vector_z_nr_slurp(cw, ZT_MPZ)
+                elif self._type == mat_gso_long_dd:
+                    vector_fp_nr_barf(cv, v, FT_DD)
+                    sig_on()
+                    (<MatGSO_c[Z_NR[long],FP_NR[dd_t]]*>self._core.long_dd).babai(cw.long, cv.dd, start, dimension, gso)
+                    sig_off()
+                    return vector_z_nr_slurp(cw, ZT_LONG)
+                # # https://github.com/fplll/fplll/issues/493
+                # elif self._type == mat_gso_long_qd:
+                #     vector_fp_nr_barf(cv, v, FT_QD)
+                #     sig_on()
+                #     (<MatGSO_c[Z_NR[long],FP_NR[qd_t]]*>self._core.long_qd).babai(cw.long, cv.qd, start, dimension, gso)
+                #     sig_off()
+                #     return vector_z_nr_slurp(cw, ZT_LONG)
+
+            raise NotImplementedError
+        else:
+            if self._type == mat_gso_mpz_d:
+                vector_fp_nr_barf(cv, v, FT_DOUBLE)
+                sig_on()
+                self._core.mpz_d.babai(cw.mpz, cv.d, start, dimension)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_MPZ)
+            elif self._type == mat_gso_mpz_ld:
+                vector_fp_nr_barf(cv, v, FT_LONG_DOUBLE)
+                sig_on()
+                self._core.mpz_ld.babai(cw.mpz, cv.ld, start, dimension)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_MPZ)
+            elif self._type == mat_gso_mpz_dpe:
+                vector_fp_nr_barf(cv, v, FT_DPE)
+                sig_on()
+                self._core.mpz_dpe.babai(cw.mpz, cv.dpe, start, dimension)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_MPZ)
+            elif self._type == mat_gso_mpz_mpfr:
+                vector_fp_nr_barf(cv, v, FT_MPFR)
+                sig_on()
+                self._core.mpz_mpfr.babai(cw.mpz, cv.mpfr, start, dimension)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_MPZ)
+            elif self._type == mat_gso_long_d:
+                vector_fp_nr_barf(cv, v, FT_DOUBLE)
+                sig_on()
+                self._core.long_d.babai(cw.long, cv.d, start, dimension)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_LONG)
+            elif self._type == mat_gso_long_ld:
+                vector_fp_nr_barf(cv, v, FT_LONG_DOUBLE)
+                sig_on()
+                self._core.long_ld.babai(cw.long, cv.ld, start, dimension)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_LONG)
+            elif self._type == mat_gso_long_dpe:
+                vector_fp_nr_barf(cv, v, FT_DPE)
+                sig_on()
+                self._core.long_dpe.babai(cw.long, cv.dpe, start, dimension)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_LONG)
+            elif self._type == mat_gso_long_mpfr:
+                vector_fp_nr_barf(cv, v, FT_MPFR)
+                sig_on()
+                self._core.long_mpfr.babai(cw.long, cv.mpfr, start, dimension)
+                sig_off()
+                return vector_z_nr_slurp(cw, ZT_LONG)
+
+            IF HAVE_QD:
+                if self._type == mat_gso_mpz_dd:
+                    vector_fp_nr_barf(cv, v, FT_DD)
+                    sig_on()
+                    self._core.mpz_dd.babai(cw.mpz, cv.dd, start, dimension)
+                    sig_off()
+                    return vector_z_nr_slurp(cw, ZT_MPZ)
+                elif self._type == mat_gso_mpz_qd:
+                    vector_fp_nr_barf(cv, v, FT_QD)
+                    sig_on()
+                    self._core.mpz_qd.babai(cw.mpz, cv.qd, start, dimension)
+                    sig_off()
+                    return vector_z_nr_slurp(cw, ZT_MPZ)
+                elif self._type == mat_gso_long_dd:
+                    vector_fp_nr_barf(cv, v, FT_DD)
+                    sig_on()
+                    self._core.long_dd.babai(cw.long, cv.dd, start, dimension)
+                    sig_off()
+                    return vector_z_nr_slurp(cw, ZT_LONG)
+                elif self._type == mat_gso_long_qd:
+                    vector_fp_nr_barf(cv, v, FT_QD)
+                    sig_on()
+                    self._core.long_qd.babai(cw.long, cv.qd, start, dimension)
+                    sig_off()
+                    return vector_z_nr_slurp(cw, ZT_LONG)
+
+            raise NotImplementedError
+
+    # def _babai_old(self, v, int start=0, int dimension=-1, gso=False):
+    #     if dimension == -1:
+    #         dimension = self.d - start
+    #     if not gso:
+    #         v = self.from_canonical(v, start, dimension)
+
+    #     cdef Py_ssize_t i, j
+    #     cdef list vv = list(v)
+    #     for i in range(dimension)[::-1]:
+    #         vv[i] = int(round(vv[i]))
+    #         for j in range(i):
+    #             vv[j] -= self.get_mu(start+i, start+j) * vv[i]
+    #     return tuple(vv)
 
     def r(self, start=0, end=-1):
         """
diff --git a/src/fpylll/fplll/svpcvp.pyx b/src/fpylll/fplll/svpcvp.pyx
index 126c2f36..6bd1e335 100644
--- a/src/fpylll/fplll/svpcvp.pyx
+++ b/src/fpylll/fplll/svpcvp.pyx
@@ -34,6 +34,7 @@ from .pruner import Pruning
 from fpylll.io cimport assign_Z_NR_mpz, mpz_get_python
 from fpylll.io import SuppressStream
 from fpylll.util import ReductionError
+from fpylll.algorithms.babai import babai
 
 from .integer_matrix cimport IntegerMatrix
 
@@ -161,18 +162,17 @@ class SVP:
     VERBOSE = SVP_VERBOSE
     OVERRIDE_BND = SVP_OVERRIDE_BND
 
-def closest_vector(IntegerMatrix B, target, method="fast", int flags=CVP_DEFAULT):
+def closest_vector(IntegerMatrix B, t, method="fast", int flags=CVP_DEFAULT):
     """Return a closest vector.
 
     The basis must be LLL-reduced with delta=``LLL.DEFAULT_DELTA`` and eta=``LLL.DEFAULT_ETA``.  The
     result is guaranteed if method = "proved", default is "fast".
 
-    :param IntegerMatrix B:
-    :param vector[Z_NR[mpz_t]] target:
+    :param IntegerMatrix B: Input lattice basis.
+    :param t: Target point (∈ ZZ^n)
     :param method: One of "fast" or "proved".
     :param int flags: Either ``CVP.DEFAULT`` or ``CVP.VERBOSE``.
     :returns coordinates of the solution vector:
-    :rtype tuple:
 
     EXAMPLE ::
 
@@ -184,6 +184,20 @@ def closest_vector(IntegerMatrix B, target, method="fast", int flags=CVP_DEFAULT
         >>> print (CVP.closest_vector(A, t))
         (-34, 109, 204, 360, -1548)
 
+       >>> from fpylll import *
+       >>> n = 10
+       >>> B = IntegerMatrix(n, n + 1)
+       >>> B.randomize("intrel", bits=100)
+       >>> v_opt = B.multiply_left([1,0,1,0,1,1,0,0,1,1])
+       >>> s = v_opt[0] # s = <a, x>, where a is vector of knapsack values.
+       >>> t = [s] + (n * [0])
+       >>> _ = LLL.reduction(B)
+       >>> v = CVP.closest_vector(B, t)
+       >>> v[0] == t[0]
+       True
+       >>> v[1:]
+       (1, 0, 1, 0, 1, 1, 0, 0, 1, 1)
+
     """
 
     if B.nrows > MAX_ENUM_DIM:
@@ -206,11 +220,10 @@ def closest_vector(IntegerMatrix B, target, method="fast", int flags=CVP_DEFAULT
     cdef vector[Z_NR[mpz_t]] sol_coord
     cdef vector[Z_NR[mpz_t]] solution
 
-    int_target.resize(len(target))
-    cdef Z_NR[mpz_t] t
+    int_target.resize(len(t))
 
-    for i in range(len(target)):
-        assign_Z_NR_mpz(int_target[i], target[i])
+    for i in range(len(t)):
+        assign_Z_NR_mpz(int_target[i], t[i])
 
     sig_on()
     r = closest_vector_c(B._core.mpz[0], int_target, sol_coord, method_, flags)
@@ -229,6 +242,7 @@ def closest_vector(IntegerMatrix B, target, method="fast", int flags=CVP_DEFAULT
     return tuple(v)
 
 class CVP:
+    babai = staticmethod(babai)
     closest_vector = staticmethod(closest_vector)
     DEFAULT = CVP_DEFAULT
     VERBOSE = CVP_VERBOSE
diff --git a/src/fpylll/io.pxd b/src/fpylll/io.pxd
index 321850c8..b5599fd0 100644
--- a/src/fpylll/io.pxd
+++ b/src/fpylll/io.pxd
@@ -3,7 +3,12 @@ include "fpylll/config.pxi"
 from .fplll.fplll cimport Z_NR
 from .gmp.mpz cimport mpz_t
 from .gmp.types cimport mpz_srcptr
+from .fplll.decl cimport fp_nr_t, vector_fp_nr_t, vector_z_nr_t
+from fplll.fplll cimport FloatType, IntType
 
 cdef int assign_Z_NR_mpz(Z_NR[mpz_t]& t, value) except -1
 cdef int assign_mpz(mpz_t& t, value) except -1
 cdef object mpz_get_python(mpz_srcptr z)
+cdef void vector_fp_nr_barf(vector_fp_nr_t &out, object inp, FloatType float_type)
+cdef object vector_fp_nr_slurp(vector_fp_nr_t &inp, FloatType float_type)
+cdef object vector_z_nr_slurp(vector_z_nr_t &inp, IntType int_type)
diff --git a/src/fpylll/io.pyx b/src/fpylll/io.pyx
index 5ac837d2..6b07502d 100644
--- a/src/fpylll/io.pyx
+++ b/src/fpylll/io.pyx
@@ -9,6 +9,11 @@ from fpylll.gmp.mpz cimport mpz_init, mpz_clear, mpz_set
 from fpylll.gmp.pylong cimport mpz_get_pyintlong, mpz_set_pylong
 from .gmp.mpz cimport mpz_t, mpz_set_si, mpz_set
 from cpython.version cimport PY_MAJOR_VERSION
+from fplll.fplll cimport FT_DEFAULT, FT_DOUBLE, FT_LONG_DOUBLE, FT_DPE, FT_MPFR
+from fplll.fplll cimport ZT_MPZ, ZT_LONG
+
+IF HAVE_QD:
+    from fpylll.fplll.fplll cimport FT_DD, FT_QD
 
 try:
     from sage.all import ZZ
@@ -55,6 +60,88 @@ cdef object mpz_get_python(mpz_srcptr z):
     else:
         return r
 
+cdef void vector_fp_nr_barf(vector_fp_nr_t &out, object inp, FloatType float_type):
+    cdef fp_nr_t tmp
+    cdef bytes py_bytes
+    if float_type == FT_DOUBLE:
+        for entry in inp:
+            # this is slow but we want to cover all kinds of Python types here
+            py_bytes = str(entry).encode()
+            tmp.d = <char*>py_bytes
+            out.d.push_back(tmp.d)
+    elif float_type == FT_LONG_DOUBLE:
+        for entry in inp:
+            py_bytes = str(entry).encode()
+            tmp.ld = <char*>py_bytes
+            out.ld.push_back(tmp.ld)
+    elif float_type == FT_DPE:
+        for entry in inp:
+            py_bytes = str(entry).encode()
+            tmp.dpe = <char*>py_bytes
+            out.dpe.push_back(tmp.dpe)
+    elif float_type == FT_MPFR:
+        for entry in inp:
+            py_bytes = str(entry).encode()
+            tmp.mpfr = <char*>py_bytes
+            out.mpfr.push_back(tmp.mpfr)
+    else:
+        IF HAVE_QD:
+            if float_type == FT_DD:
+                for entry in inp:
+                    py_bytes = str(entry).encode()
+                    tmp.dd = <char*>py_bytes
+                    out.dd.push_back(tmp.dd)
+            elif float_type == FT_QD:
+                for entry in inp:
+                    py_bytes = str(entry).encode()
+                    tmp.qd = <char*>py_bytes
+                    out.qd.push_back(tmp.qd)
+            else:
+                raise ValueError("Float type '%s' not understood."%float_type)
+        ELSE:
+            raise ValueError("Float type '%s' not understood."%float_type)
+
+cdef object vector_fp_nr_slurp(vector_fp_nr_t &inp, FloatType float_type):
+    out = []
+    if float_type == FT_DOUBLE:
+        for i in range(inp.d.size()):
+            out.append(inp.d[i].get_d())
+    elif float_type == FT_LONG_DOUBLE:
+        for i in range(inp.ld.size()):
+            out.append(inp.ld[i].get_d())
+    elif float_type == FT_DPE:
+        for i in range(inp.dpe.size()):
+            out.append(inp.dpe[i].get_d())
+    elif float_type == FT_MPFR:
+        for i in range(inp.mpfr.size()):
+            out.append(inp.mpfr[i].get_d())
+    else:
+        IF HAVE_QD:
+            if float_type == FT_DD:
+                for i in range(inp.dd.size()):
+                    out.append(inp.dd[i].get_d())
+            elif float_type == FT_QD:
+                for i in range(inp.qd.size()):
+                    out.append(inp.qd[i].get_d())
+            else:
+                raise ValueError("Float type '%s' not understood."%float_type)
+        ELSE:
+            raise ValueError("Float type '%s' not understood."%float_type)
+    return tuple(out)
+
+cdef object vector_z_nr_slurp(vector_z_nr_t &inp, IntType int_type):
+    out = []
+    if int_type == ZT_MPZ:
+        for i in range(inp.mpz.size()):
+            out.append(mpz_get_python(inp.mpz[i].get_data()))
+    elif int_type == ZT_LONG:
+        for i in range(inp.long.size()):
+            out.append(inp.long[i].get_data())
+    else:
+        raise ValueError("Int type '%s' not understood."%int_type)
+    return tuple(out)
+
+
 class SuppressStream(object):
     """
     Suppress errors (being printed by FPLLL, which are to be expected).
diff --git a/src/fpylll/util.pyx b/src/fpylll/util.pyx
index 22755bbd..40b03d87 100644
--- a/src/fpylll/util.pyx
+++ b/src/fpylll/util.pyx
@@ -3,7 +3,7 @@ include "fpylll/config.pxi"
 
 from contextlib import contextmanager
 
-from fpylll.fplll.decl cimport fp_nr_t
+
 from fpylll.fplll.fplll cimport FP_NR, RandGen, dpe_t
 from fpylll.fplll.fplll cimport FT_DEFAULT, FT_DOUBLE, FT_LONG_DOUBLE, FT_DPE, FT_MPFR
 from fpylll.fplll.fplll cimport IntType, ZT_LONG, ZT_MPZ
@@ -118,7 +118,6 @@ cdef int check_delta(float delta) except -1:
     elif delta > 1.0:
         raise TypeError("delta must be <= 1.0")
 
-
 def set_random_seed(unsigned long seed):
     """Set random seed.
 


Test patch:
diff --git a/tests/test_callback_enum.py b/tests/test_callback_enum.py
index c4340192..60ad854c 100644
--- a/tests/test_callback_enum.py
+++ b/tests/test_callback_enum.py
@@ -11,9 +11,9 @@ def test_callback_enum(d=40):
 
     # we are not imposing a constraint
     enum_obj = Enumeration(M)
-    solutions = enum_obj.enumerate(0, d, 0.99*M.get_r(0, 0), 0)
+    solutions = enum_obj.enumerate(0, d, 0.99 * M.get_r(0, 0), 0)
     max_dist, sol = solutions[0]
-    assert(A.multiply_left(sol)[0] != 2)
+    assert A.multiply_left(sol)[0] != 2
 
     # now we do
     def callback(new_sol_coord):
@@ -23,7 +23,7 @@ def callback(new_sol_coord):
             return False
 
     enum_obj = Enumeration(M, callbackf=callback)
-    solutions = enum_obj.enumerate(0, d, 0.99*M.get_r(0, 0), 0)
+    solutions = enum_obj.enumerate(0, d, 0.99 * M.get_r(0, 0), 0)
     max_dist, sol = solutions[0]
 
-    assert(A.multiply_left(sol)[0] == 2)
+    assert A.multiply_left(sol)[0] == 2
diff --git a/tests/test_gso.py b/tests/test_gso.py
index 3b0616c8..8f881e31 100644
--- a/tests/test_gso.py
+++ b/tests/test_gso.py
@@ -99,9 +99,9 @@ def test_gso_update_gso():
                 assert g00[0] ==  pytest.approx(g00[i], rel=EPSILON)
 
 
-def test_gso_io():
+def test_gso_babai():
     for int_type in int_types:
-        for m, n in dimensions:
+        for m, n in ((0, 0), (2, 2), (3, 3), (10, 10), (30, 30)):
             if m <= 2 or n <= 2:
                 continue
 
@@ -110,12 +110,34 @@ def test_gso_io():
             LLL.reduction(A)
 
             for float_type in float_types:
-                M = GSO.Mat(copy(A), float_type=float_type)
-                M.update_gso()
-                w = M.babai(v)
-                v_ = IntegerMatrix.from_iterable(1, m, w) * A
-                v_ = list(v_[0])
-                assert v == v_
+                M = GSO.Mat(copy(A), update=True, float_type=float_type)
+                try:
+                    w = M.babai(v)
+                    v_ = IntegerMatrix.from_iterable(1, m, w) * A
+                    v_ = list(v_[0])
+                    assert v == v_
+                except NotImplementedError:
+                    pass
+
+
+def test_gso_conversion():
+    for int_type in int_types:
+        for m, n in ((0, 0), (2, 2), (3, 3), (10, 10), (30, 30)):
+            if m <= 2 or n <= 2:
+                continue
+
+            A = make_integer_matrix(m, n, int_type=int_type)
+            v = list(A[0])
+            LLL.reduction(A)
+
+            for float_type in float_types:
+                M = GSO.Mat(copy(A), update=True, float_type=float_type)
+                try:
+                    w = M.from_canonical(v)
+                    v_ = [int(round(v__)) for v__ in M.to_canonical(w)]
+                    assert v == v_
+                except NotImplementedError:
+                    pass
 
 
 def test_gso_coherence_gram_matrix():