From 0ea35105753c509587285989a79fb604e0866968 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 25 Aug 2019 18:58:05 +0300 Subject: [PATCH] =?UTF-8?q?Lenstra-Lenstra-Lov=C3=A1sz=20lattice=20basis?= =?UTF-8?q?=20reduction.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Former-commit-id: 3ada9f1ebf783d07c31671fd81cb6f7d6fe6187b --- qpms/CMakeLists.txt | 4 +- qpms/__init__.py | 4 +- qpms/lattices.h | 12 ++++++ qpms/lll.c | 99 +++++++++++++++++++++++++++++++++++++++++++++ qpms/qpms_c.pyx | 27 +++++++++++++ qpms/qpms_cdefs.pxd | 1 + 6 files changed, 144 insertions(+), 3 deletions(-) create mode 100644 qpms/lll.c diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt index 3dffad0..d335f93 100644 --- a/qpms/CMakeLists.txt +++ b/qpms/CMakeLists.txt @@ -13,7 +13,9 @@ include_directories(${DIRS}) add_library (qpms SHARED translations.c tmatrices.c vecprint.c vswf.c wigner.c ewald.c ewaldsf.c pointgroups.c latticegens.c lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c - bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c) + bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c + lll.c + ) use_c99() set(LIBS ${LIBS} ${GSL_LIBRARIES} ${GSLCBLAS_LIBRARIES}) diff --git a/qpms/__init__.py b/qpms/__init__.py index 735021f..3b2c737 100644 --- a/qpms/__init__.py +++ b/qpms/__init__.py @@ -7,14 +7,14 @@ from sys import platform as __platform import warnings as __warnings try: - from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix, pitau, set_gsl_pythonic_error_handling, pgsl_ignore_error, gamma_inc + from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix, pitau, set_gsl_pythonic_error_handling, pgsl_ignore_error, gamma_inc, lll_reduce except ImportError as ex: if __platform == "linux" or __platform == "linux2": if 'LD_LIBRARY_PATH' not in __os.environ.keys(): __warnings.warn("Environment variable LD_LIBRARY_PATH has not been set. Make it point to a directory where you installed libqpms and run python again") else: __warnings.warn("Does your LD_LIBRARY_PATH include a directory where you installed libqpms? Check and run python again." - 'Currently, I see LD_LIBRARY_PATH="%s"' % __os.environ['LD_LIBRARY_PATH']) + '\nCurrently, I see LD_LIBRARY_PATH="%s"' % __os.environ['LD_LIBRARY_PATH']) raise ex from .qpms_p import * from .cyquaternions import CQuat, IRot3 diff --git a/qpms/lattices.h b/qpms/lattices.h index 1798e93..b9c28dd 100644 --- a/qpms/lattices.h +++ b/qpms/lattices.h @@ -63,6 +63,18 @@ static inline point2d point2d_fromxy(const double x, const double y) { return p; } +/// Lattice basis reduction. +/** This is currenty a bit naïve implementation of + * Lenstra-Lenstra-Lovász algorithm. + * + * The reduction happens in-place, i.e. the basis vectors in \a b are + * replaced with the reduced basis. + */ +int qpms_reduce_lattice_basis(double *b, ///< Array of dimension [bsize][ndim]. + const size_t bsize, ///< Number of the basis vectors (dimensionality of the lattice). + const size_t ndim ///< Dimension of the space into which the lattice is embedded. + ); + /// Generic lattice point generator type. /** * A bit of OOP-in-C brainfuck here. diff --git a/qpms/lll.c b/qpms/lll.c new file mode 100644 index 0000000..02fc0da --- /dev/null +++ b/qpms/lll.c @@ -0,0 +1,99 @@ +#include +#include +#include +#include + +static inline size_t mu_index(size_t k, size_t j) { + return k * (k - 1) / 2 + j; +} + +/// Gram-Schmidt orthogonalisation. +/** Does not return the actual orthogonal basis (as it is not needed + * for the LLL algorithm as such) but rather only + * the mu(i,j) coefficients and squared norms of the orthogonal vectors + */ +static void gram_schmidt( + double *mu, ///< Array of \f[ \mu_{k,j} = \frac{\vect{v}_i \cdot\vect{v}_i^*}{|\vect v_j^*|^2}\f] of length mu_index(bsize, 0), + double *vstar_sqnorm, ///< Array of \f$ \abs{\vect v_i^*}^2 \f$ of length bsize. + const double *v, ///< Vectors to orthogonalise, size [bsize][ndim], + const size_t bsize, ///< Size of the basis ( = dimensionality of the lattice) + const size_t ndim ///< Dimensionality of the space. + ) +{ + double v_star[bsize][ndim]; + for (size_t i = 0; i < bsize; ++i) { + memcpy(v_star[i], v+i*ndim, ndim*sizeof(double)); + double parallel_part[ndim /*???*/]; + memset(parallel_part, 0, sizeof(parallel_part)); + for (size_t j = 0; j < i; ++j) { + double mu_numerator = cblas_ddot(ndim, v + i*ndim, 1, v_star[j], 1); + mu[mu_index(i, j)] = mu_numerator / vstar_sqnorm[j]; + cblas_daxpy(ndim, mu[mu_index(i, j)], v_star[j], 1, parallel_part, 1); + } + cblas_daxpy(ndim, -1, parallel_part, 1, v_star[i], 1); + vstar_sqnorm[i] = cblas_ddot(ndim, v_star[i], 1, v_star[i], 1); + } +} + +static inline double fsq(double x) { return x * x; }; + +// A naïve implementation of Lenstra-Lenstra-Lovász algorithm. +int qpms_reduce_lattice_basis(double *b, const size_t bsize, const size_t ndim) +{ + QPMS_ENSURE(bsize <= ndim, + "The basis must have less elements (%zd) than the space dimension (%zd).", + bsize, ndim); + double mu[mu_index(bsize,0)]; + double bstar_sqnorm[bsize]; + gram_schmidt(mu, bstar_sqnorm, b, bsize, ndim); + size_t k = 1; + while (k < bsize) { + // Step 1 of LLL, achieve mu(k, k-1) <= 0.5 + if (fabs(mu[mu_index(k, k-1)]) > 0.5) { + double r = round(mu[mu_index(k, k-1)]); + // "normalize" b(k), replacing it with b(k) - r b(k-1) + cblas_daxpy(ndim, -r, b+(k-1)*ndim, 1, b+k*ndim, 1); + // update mu to correspond to the new b(k) + for(size_t j = 0; j < bsize; ++j) + mu[mu_index(k, j)] -= r*mu[mu_index(k-1, j)]; + } + // Step 2 + if (k > 0 && // Case 1 + bstar_sqnorm[k] < (0.75 - fsq(mu[mu_index(k, k-1)])) * bstar_sqnorm[k-1]) { + // swap b(k) and b(k-1) + cblas_dswap(ndim, &b[k*ndim], 1, &b[(k-1)*ndim], 1); + double B_k = bstar_sqnorm[k]; + double mu_kkm1_old = mu[mu_index(k, k-1)]; + double C = B_k + fsq(mu_kkm1_old) * bstar_sqnorm[k-1]; + mu[mu_index(k, k-1)] *= bstar_sqnorm[k-1] / C; + bstar_sqnorm[k] *= bstar_sqnorm[k-1] / C; + bstar_sqnorm[k-1] = C; + for(size_t j = k+1; j < bsize; ++j) { + double m = mu[mu_index(j, k-1)]; + mu[mu_index(j, k-1)] = m*m + mu[mu_index(j, k)] * B_k / C; + mu[mu_index(j, k)] = m - mu[mu_index(j, k)] * mu_kkm1_old; + } + for(size_t j = 0; j < k-1; ++j) { + double m = mu[mu_index(k-1, j)]; + mu[mu_index(k-1, j)] = mu[mu_index(k, j)]; + mu[mu_index(k, j)] = m; + } + --k; + } else { // Case 2 + size_t l = k; + while(l > 0) { + --l; + if(fabs(mu[mu_index(k, l)] > 0.5)) { + double r = round(mu[mu_index(k, l)]); + cblas_daxpy(ndim, -r, b+l*ndim, 1, b+k*ndim, 1); + for (size_t j = 0; j < l; ++j) + mu[mu_index(k, j)] -= r * mu[mu_index(l, j)]; + mu[mu_index(k, l)] -= r; + l = k; + } + } + ++k; + } + } + return 0; +} diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index d539f2b..be09550 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -609,3 +609,30 @@ def gamma_inc_CF(double a, cdouble x): with pgsl_ignore_error(15): #15 is underflow cx_gamma_inc_CF_e(a, x, &res) return (res.val, res.err) + +def lll_reduce(basis): + """ + Lattice basis reduction with the Lenstra-Lenstra-Lovász algorithm. + + basis is array_like with dimensions (n, d), where + n is the size of the basis (dimensionality of the lattice) + and d is the dimensionality of the space into which the lattice + is embedded. + """ + basis = np.array(basis, copy=True, order='C') + if len(basis.shape) != 2: + raise ValueError("Expected two-dimensional array (got %d-dimensional)" + % len(basis.shape)) + cdef size_t n, d + n, d = basis.shape + if n > d: + raise ValueError("Real space dimensionality (%d) cannot be smaller than" + "the dimensionality of the lattice (%d) embedded into it." + % (d, n)) + cdef double [:,:] basis_view = basis + if 0 != qpms_reduce_lattice_basis(&basis_view[0,0], n, d): + raise RuntimeError("Something weird happened") + return basis + + + diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 2efe72c..9987119 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -209,6 +209,7 @@ cdef extern from "lattices.h": PGEN_1D_INC_TOWARDS_ORIGIN PGen PGen_1D_new_minMaxR(double period, double offset, double minR, bint inc_minR, double maxR, bint inc_maxR, PGen_1D_incrementDirection incdir) + int qpms_reduce_lattice_basis(double *b, size_t bsize, size_t ndim) cdef extern from "quaternions.h":