Category: Handbook of Applied Cryptography

// Algorithm found in the "Handbook of
// Applied Cryptography" (c) 1997 by
// Alfred J. Menezes, Paul C van Oorschot,
// and Scott Vanstone 3.105 Algorithm
// Chapter 3 pages 120 - 121

#pragma once
class LLL_Lattice
{

private:

	static double Scalar(
		int n,
		std::vector<double> u,
		std::vector<double> v);
	static void RED(
		int k, int l, int n,
		std::vector<std::vector<double>>& b,
		std::vector<std::vector<double>>& h,
		std::vector<std::vector<double>>& mu);
	static void SWAP(
		int k, int k1, int kmax, int n,
		double m, std::vector<double>& B,
		std::vector<std::vector<double>>& b,
		std::vector<std::vector<double>>& bs,
		std::vector<std::vector<double>>& h,
		std::vector<std::vector<double>>& mu);

public:

	static bool LLL(
		int n,
		std::vector<std::vector<double>>& b,
		std::vector<std::vector<double>>& h);
};

#include "pch.h"
#include "LLL_Lattice.h"

double LLL_Lattice::Scalar(
    int n,
    std::vector<double> u,
    std::vector<double> v)
{
    // Calculate the scalar product of two vectors [1..n]
    double sum = 0.0;

    for (int i = 1; i <= n; i++) sum += u[i] * v[i];
    return sum;
}

void LLL_Lattice::RED(
    int k, int l, int n,
    std::vector<std::vector<double>>& b,
    std::vector<std::vector<double>>& h,
    std::vector<std::vector<double>>& mu)
{
    int i, q = (int)(0.5 + mu[k][l]);

    if (fabs(mu[k][l]) > 0.5)
    {
        for (i = 1; i <= n; i++)
        {
            b[k][i] -= q * b[l][i];
            h[i][k] -= q * h[i][l];
        }

        mu[k][l] -= q;

        for (i = 1; i <= l - 1; i++)
        {
            mu[k][i] -= q * mu[l][i];
        }
    }
}

void LLL_Lattice::SWAP(
    int k, int k1, int kmax, int n,
    double m, std::vector<double>& B,
    std::vector<std::vector<double>>& b,
    std::vector<std::vector<double>>& bs,
    std::vector<std::vector<double>>& h,
    std::vector<std::vector<double>>& mu)
{
    double C, t;
    int i, j;
    std::vector<double> c(n + 1);

    for (j = 1; j <= n; j++)
    {
        t = b[k][j];
        b[k][j] = b[k1][j];
        b[k1][j] = t;
        t = h[j][k];
        h[j][k] = h[j][k1];
        h[j][k1] = t;
    }

    if (k > 2)
    {
        for (j = 1; j <= k - 2; j++)
        {
            t = mu[k][j];
            mu[k][j] = mu[k1][j];
            mu[k1][j] = t;
        }
    }

    C = B[k] + m * m * B[k1];
    mu[k][k1] = m * B[k1] / C;

    for (i = 1; i <= n; i++)
    {
        c[i] = bs[k1][i];
    }

    for (i = 1; i <= n; i++)
    {
        bs[k1][i] = bs[k][i] + m * c[i];
    }
    
    for (i = 1; i <= n; i++)
    {
        bs[k][i] = -m * bs[k][i] + B[k] * c[i] / C;
    }

    B[k] *= B[k1] / C;
    B[k1] = C;

    for (i = k + 1; i <= kmax; i++)
    {
        t = mu[i][k];
        mu[i][k] = mu[i][k1] - m * t;
        mu[i][k1] = t + mu[k][k1] * mu[i][k];
    }
}

bool LLL_Lattice::LLL(
    int n,
    std::vector<std::vector<double>>& b,
    std::vector<std::vector<double>>& h)
{
    // Lattice reduction algorithm
 
    double m;
    std::vector<double> B(n + 1);
    std::vector<double> bv(n + 1);
    std::vector<double> bw(n + 1);
    std::vector<std::vector<double>> bs(n + 1,
        std::vector<double>(n + 1));
    std::vector<std::vector<double>> mu(n + 1,
        std::vector<double>(n + 1));
    int i, j, k, k1, kmax = 1, l;

    for (i = 1; i <= n; i++)
        bv[i] = bs[1][i] = b[1][i];

    B[1] = Scalar(n, bv, bv);

    for (i = 1; i <= n; i++)
    {
        for (j = 1; j <= n; j++)
        {
            h[i][j] = 0.0;
        }
        
        h[i][i] = 1.0;
    }

    for (k = 2; k <= n; k++)
    {
        if (k <= kmax)
            goto Label3;

        kmax = k;
        for (i = 1; i <= n; i++)
        {
            bs[k][i] = b[k][i];
        }

        for (j = 1; j <= k - 1; j++)
        {
            for (l = 1; l <= n; l++)
            {
                bv[l] = b[k][l];
                bw[l] = bs[j][l];
            }

            mu[k][j] = Scalar(n, bv, bw) / B[j];

            for (i = 1; i <= n; i++)
            {
                bs[k][i] -= mu[k][j] * bs[j][i];
            }
        }

        for (i = 1; i <= n; i++)
        {
            bv[i] = bs[k][i];
        }

        B[k] = Scalar(n, bv, bv);

        if (B[k] == 0.0)
            return false;

    Label3:

        k1 = k - 1;
        m = mu[k][k1];
        RED(k, k1, n, b, h, mu);

        if (B[k] < (0.75 - m * m) * B[k1])
        {
            SWAP(k, k1, kmax, n, m, B, b, bs, h, mu);
            k = max(2, k1);
            goto Label3;
        }

        for (l = k - 2; l >= 1; l--)
        {
            RED(k, l, n, b, h, mu);
        }
    }

    return true;
}