色々な言語でLLL簡約を試してみる


# 目次
[1.C言語で実装](#c言語での実装)
[2.C++で実装](#cでの実装)
[3.Pythonで実装](#pythonでの実装)
[4.SageMathで実装](#sagemathでの実装)
[5.Risa/Asirで実装](#risaasirでの実装)
[6.実行してみた](#それぞれの実行時間の比較)
[7.比較してみた](#それぞれの結果の比較)
[8.GitHubにあげました](#最後に)

:::note warn
コメントでの指摘にもある通り、Cなどでヘッダに関数を直接書くのは危ないようです。そもそもこのなんちゃって実験においてわざわざヘッダに書く必要もないのでcファイル、cppファイルに直接関数を定義することにしました。
:::


# C言語での実装

```c: LLL簡約（C言語）
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* 内積 */
double dot_dbl_dbl(double *x, double *y, const int n){
    double s = 0.0;
    for(int i = 0; i < n; ++i) s += x[i] * y[i];
    return s;
}
double dot_int_dbl(int *x, double *y, const int n){
    double s = 0.0;
    for(int i = 0; i < n; ++i) s += y[i] * x[i];
    return s;
}


/* Gram-Schmidtの直交化 */
void GSO(int **b, double *B, double **mu, const int n, const int m){
    int i, j, k;
    double t, s, **GSOb;
    GSOb = (double **)malloc(n * sizeof(double *));
    for(i = 0; i < n; ++i) GSOb[i] = (double *)malloc(m * sizeof(double));

    for(i = 0; i < n; ++i){
        mu[i][i] = 1.0;
        for(j = 0; j < m; ++j) GSOb[i][j] = b[i][j];
        for(j = 0; j < i; ++j){
            mu[i][j] = dot_int_dbl(b[i], GSOb[j], m) / dot_dbl_dbl(GSOb[j], GSOb[j], m);
            for(k = 0; k < m; ++k) GSOb[i][k] -= mu[i][j] * GSOb[j][k];
        }
        B[i] = dot_dbl_dbl(GSOb[i], GSOb[i], m);
    }
}


/* 部分サイズ基底簡約 */
void SizeReduce(int **b, double **mu, const int i, const int j, const int m){
    int k;
    if(mu[i][j] > 0.5 || mu[i][j] < -0.5){
        const int q = round(mu[i][j]);
        for(k = 0; k < m; ++k) b[i][k] -= q * b[j][k];
        for(k = 0; k <= j; ++k) mu[i][k] -= mu[j][k] * q;
    }
}

/* LLL基底簡約 */
void LLLReduce(int **b, const double d, const int n, const int m){
    int j, i, h;
    double **mu, *B, nu, BB, C, t;
    mu = (double **)malloc(n * sizeof(double *));
    B = (double *)malloc(n * sizeof(double));
    for(i = 0; i < n; ++i) mu[i] = (double *)malloc(n * sizeof(double));
    GSO(b, B, mu, n, m);

    int tmp;
    for(int k = 1; k < n;){
        h = k - 1;
        for(j = h; j > -1; --j) SizeReduce(b, mu, k, j, m);

        if(k > 0 && B[k] < (d - mu[k][h] * mu[k][h]) * B[h]){
            for(i = 0; i < m; ++i){tmp = b[h][i]; b[h][i] = b[k][i]; b[k][i] = tmp;}
            
            nu = mu[k][k - 1]; BB = B[k] + nu * nu * B[k - 1]; C = 1.0 / BB;
            mu[k][k - 1] = nu * B[k - 1] * C; B[k] *= B[k - 1] * C; B[k - 1] = BB;

            for(i = 0; i <= k - 2; ++i){
                t = mu[k - 1][i]; mu[k - 1][i] = mu[k][i]; mu[k][i] = t;
            }
            for(i = k + 1; i < n; ++i){
                t = mu[i][k]; mu[i][k] = mu[i][k - 1] - nu * t;
                mu[i][k - 1] = t + mu[k][k - 1] * mu[i][k];
            }
            
            k = h;
        }else ++k;
    }
}
```

# C++での実装

## 外部ライブラリ不使用版

```cpp: LLL簡約（C++）
#include <iostream>
#include <vector>
#include <tuple>

/* 内積 */
double dot(const std::vector<int> x, const std::vector<double> y){
	double z = 0.0;
	const int n = x.size();
	for(int i = 0; i < n; ++i) z += x.at(i) * y.at(i);
	return z;
}
double dot(const std::vector<double> x, const std::vector<double> y){
	double z = 0.0;
	const int n = x.size();
	for(int i = 0; i < n; ++i) z += x.at(i) * y.at(i);
	return z;
}
double dot(const std::vector<int> x, const std::vector<int> y){
	double z = 0.0;
	const int n = x.size();
	for(int i = 0; i < n; ++i) z += x.at(i) * y.at(i);
	return z;
}


/* Gram-Schmidtの直交化 */
std::tuple<std::vector<double>, std::vector<std::vector<double>>> Gram_Schmidt_squared(const std::vector<std::vector<int>> b){
	const int n = b.size(), m = b.at(0).size(); int i, j, k;
    std::vector<double> B(n);
	std::vector<std::vector<double>> GSOb(n, std::vector<double>(m)), mu(n, std::vector<double>(n));
	for(i = 0; i < n; ++i){
		mu.at(i).at(i)= 1.0;
		for(j = 0; j < m; ++j) GSOb.at(i).at(j) = b.at(i).at(j);
		for(j = 0; j < i; ++j){
			mu.at(i).at(j) = dot(b.at(i), GSOb.at(j)) / dot(GSOb.at(j), GSOb.at(j));
			for(k = 0; k < m; ++k) GSOb.at(i).at(k) -= mu.at(i).at(j) * GSOb.at(j).at(k);
		}
        B.at(i) = dot(GSOb.at(i), GSOb.at(i));
	}
	return std::forward_as_tuple(B, mu);
}


/* 部分サイズ基底簡約 */
void SizeReduce(std::vector<std::vector<int>>& b, std::vector<std::vector<double>>& mu, const int i, const int j){
	int q;
	const int m = b.at(0).size();
	if(mu.at(i).at(j) > 0.5 || mu.at(i).at(j) < -0.5){
		q = round(mu.at(i).at(j));
		for(int k = 0; k < m; ++k) b.at(i).at(k) -= q * b.at(j).at(k);
		for(int k = 0; k <= j; ++k) mu.at(i).at(k) -= mu.at(j).at(k) * q;
	}
}


/* LLL基底簡約 */
void LLLReduce(std::vector<std::vector<int>>& b, const float d = 0.99){
	const int n = b.size(), m = b.at(0).size(); int j, i, h;
	double t, nu, BB, C;
	std::vector<std::vector<double>> mu;
	std::vector<double> B; std::tie(B, mu) = Gram_Schmidt_squared(b);
	int tmp;
	for(int k = 1; k < n;){
		h = k - 1;

		for(j = h; j > -1; --j) SizeReduce(b, mu, k, j);

		//Checks if the lattice basis matrix b satisfies Lovasz condition.
		if(k > 0 && B.at(k) < (d - mu.at(k).at(h) * mu.at(k).at(h)) * B.at(h)){
			for(i = 0; i < m; ++i){tmp = b.at(h).at(i); b.at(h).at(i) = b.at(k).at(i); b.at(k).at(i) = tmp;}

			nu = mu.at(k).at(h); BB = B.at(k) + nu * nu * B.at(h); C = 1.0 / BB;
            mu.at(k).at(h) = nu * B.at(h) * C; B[k] *= B.at(h) * C; B.at(h) = BB;

            for(i = 0; i <= k - 2; ++i){
                t = mu.at(h).at(i); mu.at(h).at(i) = mu.at(k).at(i); mu.at(k).at(i) = t;
            }
            for(i = k + 1; i < n; ++i){
                t = mu.at(i).at(k); mu.at(i).at(k) = mu.at(i).at(h) - nu * t;
                mu.at(i).at(h) = t + mu.at(k).at(h) * mu.at(i).at(k);
            }

			--k;
		}else ++k;
	}
}
```

## Eigenライブラリ
```cpp: LLL簡約（Eigen）
#include <iostream>
#include <eigen3/Eigen/Dense>

/* Gram-Schmidtの直交化 */
void GSO(const Eigen::MatrixXi b, Eigen::VectorXd& B, Eigen::MatrixXd& mu, const int n, const int m){
    int j;
    Eigen::MatrixXd GSOb(n, m);

    for(int i = 0; i < n; ++i){
        mu.coeffRef(i, i) = 1.0;
        GSOb.row(i) = b.row(i).cast<double>();
        for(j = 0; j < i; ++j){
            mu.coeffRef(i, j) = b.row(i).cast<double>().dot(GSOb.row(j)) / GSOb.row(j).dot(GSOb.row(j));
            GSOb.row(i) -= mu.coeff(i, j) * GSOb.row(j);
        }
        B.coeffRef(i) = GSOb.row(i).dot(GSOb.row(i));
    }
}

/* サイズ基底簡約 */
void SizeReduce(Eigen::MatrixXi& b, Eigen::MatrixXd& mu, const int i, const int j, const int m){
    if(mu.coeff(i, j) > 0.5 || mu.coeff(i, j) < -0.5){
        const int q = round(mu.coeff(i, j));
        b.row(i) -= q * b.row(j);
        mu.row(i).head(j + 1) -= (double)q * mu.row(j).head(j + 1);
    }
}

/* LLL簡約 */
void LLLReduce(Eigen::MatrixXi& b, const long double d, const int n, const int m){
    double nu, BB, C, t;
    Eigen::VectorXd B(n), logB(n);
    Eigen::MatrixXd mu(n, n);
    GSO(b, B, mu, n, m);
    
    for(int k = 1, j, i, k1; k < n;){
        k1 = k - 1;
        for(j = k1; j > -1; --j) SizeReduce(b, mu, k, j, m);

        if(k > 0 && B.coeff(k) < (d - mu.coeff(k, k1) * mu.coeff(k, k1)) * B.coeff(k1)){
            b.row(k).swap(b.row(k1));
            
            nu = mu.coeff(k, k1); BB = B.coeff(k) + nu * nu * B.coeff(k1); C = 1.0 / BB;
            mu.coeffRef(k, k1) = nu * B.coeff(k1) * C;
            B.coeffRef(k) *= B.coeff(k1) * C; B.coeffRef(k1) = BB;

            mu.row(k1).head(k - 1).swap(mu.row(k).head(k - 1));
            for(i = k + 1; i < n; ++i){
                t = mu.coeff(i, k); mu.coeffRef(i, k) = mu.coeff(i, k1) - nu * t;
                mu.coeffRef(i, k1) = t + mu.coeff(k, k1) * mu.coeff(i, k);
            }
            
            k = k1;
        }else ++k;
    }
}
```

## NTLライブラリ

NTLライブラリには初めからLLL簡約する関数が用意されているので、それを利用します。また、```b```には基底行列が入っているものとします。

```cpp: LLL簡約（NTL）
#include <NTL/ZZ.h>
#include <NTL/mat_ZZ.h>
#include <NTL/LLL.h>

NTL::ZZ _;
NTL::LLL(_, b, 99, 100);
```

# Pythonでの実装
```python: LLL.py
import numpy as np

def Gram_Schmidt(b):
	m: int = b.shape[1]; n: int = b.shape[0]
	GSOb = np.zeros((n, m)); mu = np.eye(n)
	for i in range(n):
		GSOb[i] = b[i][:]
		for j in range(i):
			mu[i, j] = (b[i] @ GSOb[j]) / (GSOb[j] @ GSOb[j])
			GSOb[i] -= mu[i, j] * GSOb[j][:]
	return GSOb, mu

 
# 部分サイズ基底簡約
def SizeReduce(b, mu, i, j):
	n = b.shape[0]; m = b.shape[1]
	if abs(mu[i, j]) > 0.5:
		q = round(mu[i, j])
		b[i] -= q * b[j][:]
		mu[i][: j + 1] -= q * mu[j][: j + 1].copy()
	return b, mu


# GSO情報の更新（LLL）
def GSOupdate_LLL(mu, B, k):
	nu = mu[k, k - 1]
	b = B[k] + nu ** 2 * B[k - 1]; b_inv = 1. / b
	mu[k, k - 1] = nu * B[k - 1] * b_inv
	B[k] *= B[k - 1] * b_inv
	B[k - 1] = b

	tmp = mu[k - 1][: k - 1].copy()
	mu[k - 1][: k -1 ] = mu[k][: k - 1].copy()
	mu[k][: k - 1] = tmp.copy()
	#for j in range(k - 1):
	#	mu[k - 1, j], mu[k, j] = mu[k, j], mu[k - 1, j]
	for i in range(k + 1, B.size):
		t = mu[i, k]
		mu[i, k] = mu[i, k - 1] - nu * t
		mu[i, k - 1] = t + mu[k, k - 1] * mu[i, k]
	return mu, B


# LLL基底簡約
def LLL_reduce(b, d):
	n = b.shape[0]
	GSOb, mu = Gram_Schmidt(b)
	B = np.zeros(n)
	B = np.sum(GSOb * GSOb, axis = 1)
	k = 1
	while k < n:
		#partial size reduction
		for j in range(k)[::-1]: b, mu = SizeReduce(b, mu, k, j)
		#Satisfies Lovasz condition or not
		if B[k] >= (d - mu[k, k - 1] * mu[k, k - 1]) * B[k - 1]:
			k += 1
		else:
			b[k], b[k - 1] = b[k - 1], b[k].copy()
			mu, B = GSOupdate_LLL(mu, B, k)
			k = max(k - 1, 1)
	return b
```

# SageMathでの実装

SageにはLLL簡約を行う関数が最初から用意されているのでそれを利用します。また、```b```には基底行列が入っているものとします。
```python: LLL簡約（SageMath）
b = b.LLL(delta = 0.99)
```

# Risa/Asirでの実装

```: LLL.rr
/* 内積 */
def dot(X, Y){
	N = size(X);
	S = 0;
	for(I = 0; I < N[0]; I += 1){S += X[I] * Y[I];}
	return S;
}

def max(A, B){
    if(A > B) return A;
    else return B;
}

/* グラムシュミットの直交化法 */
def gram_schmidt_squared(B){
	M = size(B); N = M[0]; M = M[1];
	BB = newvect(N);
	MU = newmat(N, N);
	GSOB = newmat(N, M);
	for(I = 0; I < N; I += 1){
		MU[I][I] = 1.0;
		for(J = 0; J < M; J += 1){GSOB[I][J] = B[I][J];}
		for(J = 0; J < I; J += 1){
			MU[I][J] = dot(B[I], GSOB[J]) / dot(GSOB[J], GSOB[J]);
			for(K = 0; K < M; K += 1){GSOB[I][K] -= MU[I][J] * GSOB[J][K];}
		}
		BB[I] = dot(GSOB[I], GSOB[I]);
	}
	return [BB, MU];
}


/* 部分サイズ基底簡約 */
def size_reduce(B, MU, I, J){
	M = size(B); N = M[0]; M = M[1];
	if(MU[I][J] > 0.5 || MU[I][J] < -0.5){
		Q = rint(MU[I][J]);
		for(K = 0; K < M; K += 1){B[I][K] -= Q * B[J][K];}
		for(K = 0; K <= J; K += 1){MU[I][K] -= MU[J][K] * Q;}
	}
	return [B, MU];
}


/* LLL基底簡約 */
def lll(B, D){
	M = size(B); N = M[0]; M = M[1];
	MU = gram_schmidt_squared(B); BB = MU[0]; MU = MU[1];
	for(K = 1; K < N;){
		for(J = K - 1; J > -1; J -= 1){
			MU = size_reduce(B, MU, K, J);
			B = MU[0]; MU = MU[1];
		}

		if(BB[K] >= (D - MU[K][K - 1] * MU[K][K - 1]) * BB[K - 1]){K += 1;}
		else{
			for(I = 0; I < M; I += 1){
				TMP = B[K - 1][I]; B[K - 1][I] = B[K][I]; B[K][I] = TMP;
			}
			NU = MU[K][K - 1]; C = BB[K] + NU * NU * BB[K - 1]; C_INV = 1 / C; 
			MU[K][K - 1] = NU * BB[K - 1] * C_INV;
			BB[K] *= BB[K - 1] * C_INV; BB[K - 1] = C;
			for(J = 0; J < K - 1; J += 1){
				T = MU[K - 1][J]; MU[K - 1][J] = MU[K][J]; MU[K][J] = T;
			}
			for(I = K + 1; I < N; I += 1){
				T = MU[I][K]; MU[I][K] = MU[I][K - 1] - NU * T;
				MU[I][K - 1] = T + MU[K][K - 1] * MU[I][K];
			}

			K = max(K - 1, 1);
		}
	}
	return B;
}
end$
```

# それぞれの実行時間の比較
:::note info
ここでは、ランダムに生成した完全階数の40次元格子基底をLLL基底簡約できるまでにかかった時間を比較します。
:::

## 格子基底の生成方法

SageMathを用いて生成する：

```py: GenBasis.sage
A = random_matrix(ZZ, 40, 40)
```

実際に使う格子基底
```math
A=\begin{bmatrix}
-1 & 3 & 0 & -1 & -1 & 0 & -6 & 0 & -28 & -2 & 0 & -1 & 0 & -1 & 0 & -15 & 0 & -6 & 0 & -1 & -1 & -1 & -1 & -28 & -99 & -2 & 2 & -1 & -1 & -1 & 1 & 11 & 1 & -3 & 0 & 0 & 2 & 1 & -2 & 0 \\
1 & -1 & 13 & 1 & -1 & -1 & 3 & -10 & 1 & 1 & -3 & 1 & 1 & 2 & 2 & -1 & -1 & -1 & 1 & 0 & 2 & -2 & -1 & 1 & 1 & 3 & 1 & -1 & 2 & -2 & -9 & -3 & 1 & -3 & 0 & -2 & -1 & -4 & -1 & 1 \\
-3 & 0 & -3 & 1 & 1 & -3 & -3 & 6 & 2 & -1 & 1 & 0 & 0 & -4 & -1 & -1 & -4 & 1 & 1 & 3 & 0 & 2 & -1 & -1 & -9 & -1 & -1 & -1 & 1 & 0 & 0 & 1 & 0 & -1 & 1 & -1 & 0 & 2 & 13 & 0 \\
-1 & -1 & 1 & 1 & 0 & 0 & -2 & 0 & -9 & 2 & 1 & -1 & -1 & 1 & 2 & -2 & 1 & 1 & -1 & -20 & 0 & 1 & 0 & -1 & -1 & -1 & -1 & 296 & -3 & -1 & 1 & 2 & -11 & 2 & 8 & 14 & -1 & 0 & -3 & 1 \\
-2 & 0 & -1 & 1 & -1 & -1 & 1 & -1 & 1 & 6 & -1 & 0 & -1 & 1 & -11 & 1 & 1 & -1 & 1 & 1 & -3 & -1 & -1 & -1 & 4 & 5 & -4 & -1 & 3 & 0 & 1 & -5 & -1 & -2 & 1 & -1 & 1 & -3 & 0 & -5 \\
0 & 0 & 6 & 2 & -1 & 1 & 1 & 0 & 1 & 12 & 1 & -6 & 3 & 1 & 0 & -1 & 0 & 0 & 0 & 2 & 0 & 3 & -1 & 0 & 0 & -1 & 1 & 1 & 0 & 12 & -1 & 0 & 119 & -2 & -11 & 29 & 0 & 1 & 2 & -2 \\
0 & 1 & 14 & -2 & 1 & 2 & 1 & -3 & -1 & 0 & 1 & 1 & -3 & 0 & 2 & -2 & -1 & -1 & -2 & 0 & 1 & 15 & 14 & -4 & 1 & 0 & -1 & 1 & 1 & 1 & 4 & -1 & 17 & 1 & -2 & 1 & 3 & 0 & 3 & -1 \\
1 & -3 & -1 & 0 & -1 & 1 & -3 & 13 & -1 & -1 & 3 & -2 & 1 & 1 & 1 & 1 & -3 & 2 & -6 & -13 & 1 & 0 & 1 & -2 & -1 & 2 & 0 & -1 & -3 & 0 & -1 & 11 & -3 & 11 & -4 & 1 & 1 & -2 & -4 & -3 \\
0 & -7 & -1 & 0 & 1 & 0 & -37 & 0 & 0 & 1 & -16 & 0 & 1 & -2 & 0 & -7 & 5 & 2 & 2 & -1 & 1 & 2 & -12 & 0 & 0 & 4 & -1 & -1 & 1 & 17 & 13 & 0 & 2 & 2 & -1 & 0 & 0 & -1 & 5 & -8 \\
1 & 1 & -1 & 1 & 3 & 0 & -1 & 2 & -1 & 2 & -1 & -7 & 2 & -1 & 1 & -4 & 3 & 1 & 2 & -1 & 0 & 4 & 1 & -11 & 6 & -144 & 1 & -3 & 0 & -35 & 2 & 3 & 1 & -1 & 1 & 4 & 12 & 3 & 0 & 2 \\
-1 & -1 & 1 & -1 & 9 & 0 & 1 & 1 & 0 & 2 & 6 & -2 & -2 & 1 & -1 & 142 & 0 & 1 & 0 & 3 & 2 & -1 & 1 & 0 & -1 & -2 & -12 & -13 & -1 & 3 & 5 & 2 & -4 & -1 & -2 & -2 & -3 & 0 & 0 & 1 \\
0 & 17 & 0 & 2 & -16 & -3 & 0 & 0 & 1 & 26 & 1 & -3 & 3 & -2 & -4 & 0 & -2 & 0 & 1 & 1 & -1 & -6 & 1 & -5 & 1 & 1 & 1 & 0 & -1 & -3 & 18 & 4 & 2 & -3 & -107 & -5 & -6 & 0 & -3 & 1 \\
2 & -2 & 3 & 1 & 3 & -7 & -1 & 1 & 1 & 0 & 7 & -5 & 0 & -1 & 0 & 68 & 2 & 1 & -1 & 4 & 0 & 7 & 0 & 3 & 11 & 1 & -1 & -2 & 0 & 0 & 26 & -1 & 5 & -4 & -1 & 4 & -1 & -3 & 1 & -1 \\
1 & 1 & -5 & 1 & -1 & -1 & 0 & 1 & -40 & 1 & -1 & -1 & 1 & 1 & 2 & 1 & 0 & -30 & 1 & -1 & 0 & 2 & 1 & 0 & -2 & -1 & 0 & -6 & 1 & 2 & 1 & 2 & 1 & -1 & 0 & 0 & 1 & -2 & 1 & 3 \\
-7 & 5 & -2 & 2 & 3 & 0 & -1 & 1 & 1 & -218 & 2 & -1 & -1 & 2 & -2 & -1 & 0 & -1 & 0 & 5 & 1 & -1 & -1 & 1 & 0 & 1 & 1 & 1 & 1 & 1 & 0 & -3 & 3 & -1 & 1 & -1 & 0 & 1 & -2 & 1 \\
0 & 0 & -3 & -4 & 1 & -1 & -12 & -3 & 9 & 2 & 0 & -1 & -15 & -1 & -1 & 1 & -1 & -3 & 0 & -1 & -1 & 0 & 1 & 0 & 1 & 0 & -9 & 0 & 0 & 6 & 0 & 1 & -1 & -1 & -10 & -1 & 23 & 8 & 0 & 5 \\
2 & 1 & -4 & 0 & 0 & 1 & 1 & -1 & 2 & 6 & 1 & 2 & 4 & 0 & 2 & -1 & 22 & -1 & 0 & 1 & -3 & 2 & 1 & 0 & 1 & 0 & -3 & -2 & 0 & 4 & 0 & -2 & 1 & 3 & 1 & 1 & 0 & 1 & 1 & 0 \\
-12 & -1 & 0 & 1 & 0 & 1 & 1 & -1 & 1 & -27 & -1 & 2 & 1 & 1 & 1 & 5 & -1 & 2 & 1 & 0 & 2 & -8 & 0 & 0 & 11 & 95 & 2 & 1 & -1 & 1 & -1 & -2 & -3 & 0 & -2 & -2 & 2 & -1 & -3 & -2 \\
9 & -1 & -1 & 1 & -1 & 4 & -1 & -1 & 1 & 0 & -3 & 1 & 1 & 6 & 1 & -1 & -1 & 1 & 1 & 0 & 0 & 0 & 0 & -1 & 3 & 2 & 0 & -1 & -1 & -1 & 3 & 1 & 0 & -1 & 0 & -1 & -1 & -3 & 1 & -1 \\
-5 & -4 & -7 & -2 & 1 & -1 & 3 & 2 & -3 & -1 & 1 & 0 & 14 & 1 & 0 & 0 & 0 & 1 & -1 & 3 & -1 & -1 & 2 & -1 & -4 & 1 & -4 & 1 & 0 & 0 & 142 & -1 & 2 & -1 & 0 & -10 & 3 & -4 & 1 & -2 \\
0 & 1 & 0 & 1 & -2 & -1 & 0 & 0 & 1 & 0 & -1 & 1 & 2 & 0 & 0 & 3 & -1 & 1 & -1 & 0 & -1 & 0 & 156 & 5 & 0 & -3 & -3 & 7 & 1 & 184 & 4 & 1 & -33 & 1 & -5 & 0 & 8 & 0 & 1 & 0 \\
-1 & -1 & 1 & -1 & 2 & 0 & 0 & 1 & 1 & 5 & -47 & 1 & 1 & -3 & -15 & 0 & 3 & 1 & 2 & 3 & 1 & -2 & -9 & -15 & 49 & -3 & 0 & -16 & -1 & 0 & -1 & 1 & 0 & -1 & 3 & 1 & -1 & 6 & -29 & -1 \\
0 & -7 & -1 & -1 & -3 & 2 & -2 & 20 & 0 & -1 & 170 & 30 & 0 & -1 & 1 & 1 & -1 & 3 & 1 & 2 & -1 & 2 & 4 & 1 & -2 & 3 & -1 & 0 & 1 & 0 & -7 & 1 & -1 & 323 & -1 & -1 & 0 & 0 & 8 & -22 \\
-1 & -18 & 1 & -2 & -4 & -1 & 1 & 0 & 1 & -2 & 1 & -1 & 7 & -1 & 1 & -1 & 3 & -1 & -3 & 0 & 3 & 2 & 3 & -1 & -1 & 1 & 0 & -1 & 1 & 3 & -1 & 0 & -12 & -1 & 4 & 0 & -3 & 1 & -1 & 0 \\
1 & 1 & -17 & 0 & -1 & 3 & 4 & -1 & -4 & -1 & 2 & 1 & -2 & 2 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 0 & 2 & 3 & 5 & 1 & 2 & -1 & 2 & 1 & -4 & 0 & 0 & 4 & 1 & 24 & 0 & -1 & 2 \\
0 & -2 & 8 & 2 & -1 & 4 & 1 & -1 & 10 & 1 & -1 & 2 & 3 & 0 & 7 & -1 & 8 & -1 & -6 & 5 & -1 & 1 & -123 & -1 & 1 & 3 & -1 & -2 & 0 & 1 & 0 & 0 & 1 & 1 & 0 & 1 & 1 & 1 & 1 & -1 \\
3 & 1 & 7 & -1 & -3 & -1 & 3 & -1 & 3 & -1 & -65154 & -1 & 1 & 2 & -2 & 3 & 1 & -23 & 10 & -2 & 0 & 0 & 0 & 3 & -1 & 2 & 0 & -1 & -1 & 8 & -5 & 1 & 0 & -2 & 0 & 0 & 0 & 1 & -1 & 3 \\
0 & -1 & -4 & 2 & 4 & 19 & 3 & 0 & -1 & 1 & -2 & -2 & 2 & -1 & 0 & -1 & 0 & -1 & -1 & 0 & 1 & -1 & -1 & -2 & 0 & 8 & 0 & 0 & 0 & -8 & -1 & -1 & 1 & 2 & 1 & 3 & -1 & 9 & -2 & 1 \\
1 & -1 & 1 & 1 & 25 & 2 & -2 & 0 & 4 & 1 & 1 & 0 & -1 & 0 & 1 & -6 & 0 & 1 & -1 & -1 & -3 & 0 & 1 & -2 & 59 & 0 & 1 & -2 & 1 & 1 & -1 & -1 & 0 & 2 & 1 & -1 & 1 & 0 & 1 & -3 \\
1 & 4 & 1 & 4 & 1 & -2 & -58 & 0 & -9 & -1 & -1 & -1 & -1 & 0 & -20 & 0 & -1 & 0 & -2 & 0 & 1 & 4 & 2 & 5 & 19 & -1 & -3 & 1 & 140 & 1 & -2 & -1 & 0 & -1 & 3 & 2 & -1 & 0 & -2 & -1 \\
18 & -1 & 1 & -3 & -1 & -1 & 0 & 2 & 14 & -2 & -1 & 1 & 0 & -1 & 2 & -1 & -1 & 0 & -4 & 1 & 1 & -7 & -2 & -1 & -4 & 2 & 3 & -1 & -2 & -1 & -6 & -1 & -1 & 0 & -39 & 3 & 0 & 0 & -2 & 40 \\
-5 & 1 & 0 & -1 & 2 & 0 & 3 & 0 & 1 & 6 & 1 & 0 & 2 & 0 & -3 & -4 & -1270 & 1 & -1 & 1 & -2 & 2 & -2 & -1 & -4 & 1 & -1 & -2 & 1 & -2 & 2 & 0 & -3 & 2 & 0 & 1 & -1 & 0 & -1 & 1 \\
-49 & 0 & 0 & -1 & 1 & 1 & -3 & -9 & 0 & 0 & 0 & -2 & 1 & 1 & -47 & 0 & 0 & 0 & 4 & 6 & -2 & -2 & 0 & 2 & 1 & 2 & 1 & -1 & 1 & 4 & 7 & 0 & 1 & -1 & 0 & -1 & -1 & 1 & 0 & -1 \\
-2 & -1 & -2 & 1 & 17 & 1 & 10 & 6 & 0 & -1 & -1 & -2 & -3 & 1 & 1 & 1 & 5 & 18 & 1 & 0 & 49 & 7 & 0 & 2 & 0 & 0 & -1 & 0 & 0 & 0 & 2 & -2 & 1 & 0 & 0 & -1 & 0 & 0 & -2 & 3 \\
-1 & -1 & 166 & 0 & -46 & 0 & -2 & 0 & 1 & 1 & 15 & -2 & 1 & 12 & 1 & 4 & 0 & -2 & 1 & -1 & 0 & -1 & 0 & 0 & -25 & -6 & 0 & -1 & -4 & 0 & 1 & 2 & 78 & -15 & 2 & 0 & 0 & -1 & 17 & 17 \\
-13 & 4 & -16 & -1 & 0 & 1 & 0 & 1222 & 1 & -1 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0 & -1 & 0 & -5 & 3 & 1 & 1 & -1 & -1 & -1 & -1 & 3 & 0 & 0 & -4 & 2 & -1 \\
1 & 2 & -2 & 0 & 0 & -4 & 1 & -1 & 9 & -1 & -1 & -13 & 24 & -1 & -1 & -8 & 1 & -1 & -29 & 0 & 1 & 0 & -1 & 1 & 1 & 1 & 1 & 4 & 9 & -1 & -1 & 0 & -18 & -1 & -1 & 1 & 6 & -1 & -4 & 0 \\
6 & -13 & 0 & 0 & 0 & 15 & 0 & -3 & 0 & -2 & 0 & -2 & 1 & -1 & -2 & -1 & 0 & 1 & -2 & 0 & 1 & -1 & 1 & -1 & -9 & 1 & 6 & -5 & -9 & 0 & 1 & -2 & -1 & -2 & 9 & 1 & 1 & 3 & 0 & 1 \\
1 & 3 & -2 & -1 & -3 & 1 & -2 & -1 & -1 & -2 & 460 & 1 & 2 & -8 & 0 & 2 & 0 & 4 & 0 & 1 & -1 & 0 & 0 & -1 & -1 & 4 & 23 & -1 & -1 & -1 & -1 & -4 & 2 & 15 & -2 & 3 & 0 & -1 & 2 & -1 \\
0 & 0 & -1 & 1 & 69 & 4 & 0 & -2 & 0 & 8 & 1 & 4 & 1 & 2 & -2 & 0 & -2 & -1 & 8 & -4 & 1 & -11 & 0 & -1 & -1 & 0 & -1 & 0 & 1 & -1 & -11 & 0 & 505 & 3 & 8 & -1 & -1 & 1 & -3 & 1
\end{bmatrix}\in\mathbb{Z}^{40\times 40}
```

## c言語での結果

```shell
$ clang LLL.c -o LLL.exe -lm
$ time ./LLL.exe

9, -1, -1, 1, -1, 4, -1, -1, 1, 0, -3, 1, 1, 6, 1, -1, -1, 1, 1, 0, 0, 0, 0, -1, 3, 2, 0, -1, -1, -1, 3, 1, 0, -1, 0, -1, -1, -3, 1, -1,
-2, 0, -1, 1, -1, -1, 1, -1, 1, 6, -1, 0, -1, 1, -11, 1, 1, -1, 1, 1, -3, -1, -1, -1, 4, 5, -4, -1, 3, 0, 1, -5, -1, -2, 1, -1, 1, -3, 0, -5,
-3, 0, -3, 1, 1, -3, -3, 6, 2, -1, 1, 0, 0, -4, -1, -1, -4, 1, 1, 3, 0, 2, -1, -1, -9, -1, -1, -1, 1, 0, 0, 1, 0, -1, 1, -1, 0, 2, 13, 0,
...(略)...
6, -15, 15, -1, 43, 62, -20, 50, 98, 31, 30, 144, 35, -89, -11, -5, -3, -104, -3, 133, 44, 16, 13, 72, -19, 0, -25, 33, -1, -43, 12, 71, -7, -1, -6, 24, 31, -138, -66, -10,
norm = 14.035669

real    0m0.003s
user    0m0.003s
sys     0m0.000s
```

## C++での結果
```shell
$ clang++ LLL.cpp -o LLL.exe
$ time ./LLL.exe

9, -1, -1, 1, -1, 4, -1, -1, 1, 0, -3, 1, 1, 6, 1, -1, -1, 1, 1, 0, 0, 0, 0, -1, 3, 2, 0, -1, -1, -1, 3, 1, 0, -1, 0, -1, -1, -3, 1, -1,
-2, 0, -1, 1, -1, -1, 1, -1, 1, 6, -1, 0, -1, 1, -11, 1, 1, -1, 1, 1, -3, -1, -1, -1, 4, 5, -4, -1, 3, 0, 1, -5, -1, -2, 1, -1, 1, -3, 0, -5,
-3, 0, -3, 1, 1, -3, -3, 6, 2, -1, 1, 0, 0, -4, -1, -1, -4, 1, 1, 3, 0, 2, -1, -1, -9, -1, -1, -1, 1, 0, 0, 1, 0, -1, 1, -1, 0, 2, 13, 0,
...(略)...
6, -15, 15, -1, 43, 62, -20, 50, 98, 31, 30, 144, 35, -89, -11, -5, -3, -104, -3, 133, 44, 16, 13, 72, -19, 0, -25, 33, -1, -43, 12, 71, -7, -1, -6, 24, 31, -138, -66, -10,
norm = 14.035669

real    0m0.016s
user    0m0.006s
sys     0m0.007s
```

## C++での結果（Eigenライブラリ使用）
```shell
$ clang++ EigenLLL.cpp -o EigenLLL.exe
$ time ./EigenLLL.exe

   9   -1   -1    1   -1    4   -1   -1    1    0   -3    1    1    6    1   -1   -1    1    1    0    0    0    0   -1    3    2    0   -1   -1   -1    3    1    0   -1    0   -1   -1   -3    1   -1
  -2    0   -1    1   -1   -1    1   -1    1    6   -1    0   -1    1  -11    1    1   -1    1    1   -3   -1   -1   -1    4    5   -4   -1    3    0    1   -5   -1   -2    1   -1    1   -3    0   -5
  -3    0   -3    1    1   -3   -3    6    2   -1    1    0    0   -4   -1   -1   -4    1    1    3    0    2   -1   -1   -9   -1   -1   -1    1    0    0    1    0   -1    1   -1    0    2   13    0
...(略)...
   6  -15   15   -1   43   62  -20   50   98   31   30  144   35  -89  -11   -5   -3 -104   -3  133   44   16   13   72  -19    0  -25   33   -1  -43   12   71   -7   -1   -6   24   31 -138  -66  -10
norm = 14.035669

real    0m0.018s
user    0m0.016s
sys     0m0.000s
```

## C++での結果（NTLライブラリ使用）
```shell
$ clang++ NTLLLL.cpp -o NTLLLL.exe -lntl
$ time ./NTLLLL.exe

[[9 -1 -1 1 -1 4 -1 -1 1 0 -3 1 1 6 1 -1 -1 1 1 0 0 0 0 -1 3 2 0 -1 -1 -1 3 1 0 -1 0 -1 -1 -3 1 -1]
[-2 0 -1 1 -1 -1 1 -1 1 6 -1 0 -1 1 -11 1 1 -1 1 1 -3 -1 -1 -1 4 5 -4 -1 3 0 1 -5 -1 -2 1 -1 1 -3 0 -5]        
[-3 0 -3 1 1 -3 -3 6 2 -1 1 0 0 -4 -1 -1 -4 1 1 3 0 2 -1 -1 -9 -1 -1 -1 1 0 0 1 0 -1 1 -1 0 2 13 0]
...(略)...
[6 -15 15 -1 43 62 -20 50 98 31 30 144 35 -89 -11 -5 -3 -104 -3 133 44 16 13 72 -19 0 -25 33 -1 -43 12 71 -7 -1 -6 24 31 -138 -66 -10]
]
norm = 14.03566885

real    0m0.006s
user    0m0.005s
sys     0m0.000s
```

## Pythonでの結果

```shell
$ time python3 LLL.py

[[   9   -1   -1 ...   -3    1   -1]
 [  -2    0   -1 ...   -3    0   -5]
 [  -3    0   -3 ...    2   13    0]
 ...
 [  16   32    6 ...   54   13  -34]
 [  23   18    9 ...   13  -15   35]
 [   6  -15   15 ... -138  -66  -10]]
norm = 14.035668847618199

real    0m0.274s
user    0m0.525s
sys     0m0.909s
```

## SageMathでの結果
```shell
$ time sage LLL.sage
[   9   -1   -1    1   -1    4   -1   -1    1    0   -3    1    1    6    1   -1   -1    1    1    0    0    0    0   -1    3    2    0   -1   -1   -1    3    1    0   -1    0   -1   -1   -3    1   -1]
[  -2    0   -1    1   -1   -1    1   -1    1    6   -1    0   -1    1  -11    1    1   -1    1    1   -3   -1   -1   -1    4    5   -4   -1    3    0    1   -5   -1   -2    1   -1    1   -3    0   -5]
[  -3    0   -3    1    1   -3   -3    6    2   -1    1    0    0   -4   -1   -1   -4    1    1    3    0    2   -1   -1   -9   -1   -1   -1    1    0    0    1    0   -1    1   -1    0    2   13    0]
...(略)...
[   6  -15   15   -1   43   62  -20   50   98   31   30  144   35  -89  -11   -5   -3 -104   -3  133   44   16   13   72  -19    0  -25   33   -1  -43   12   71   -7   -1   -6   24   31 -138  -66  -10]
norm = 14.0356688476182

real    0m0.897s
user    0m0.781s
sys     0m0.071s
```

## Risa/Asirでの結果
```shell
$ asir
[2077] load("LLL.rr");
[2084] main();
[ 9 -1 -1 1 -1 4 -1 -1 1 0 -3 1 1 6 1 -1 -1 1 1 0 0 0 0 -1 3 2 0 -1 -1 -1 3 1 0 -1 0 -1 -1 -3 1 -1 ]
[ -2 0 -1 1 -1 -1 1 -1 1 6 -1 0 -1 1 -11 1 1 -1 1 1 -3 -1 -1 -1 4 5 -4 -1 3 0 1 -5 -1 -2 1 -1 1 -3 0 -5 ]      
[ -3 0 -3 1 1 -3 -3 6 2 -1 1 0 0 -4 -1 -1 -4 1 1 3 0 2 -1 -1 -9 -1 -1 -1 1 0 0 1 0 -1 1 -1 0 2 13 0 ]
...(略)...
[ 6 -15 15 -1 43 62 -20 50 98 31 30 144 35 -89 -11 -5 -3 -104 -3 133 44 16 13 72 -19 0 -25 33 -1 -43 12 71 -7 -1 -6 24 31 -138 -66 -10 ]
14.0357
0.6324sec(0.6324sec)
0
```

# それぞれの結果の比較

| 言語 | 実行時間(秒) |
|:-:|:-:|
| C言語 | 0.003 |
| C++(外部ライブラリ不使用)  | 0.016  |
| C++(外部ライブラリ不使用、O3)  | 0.004  |
| C++(Eigenライブラリ使用)  | 0.018  |
| C++(Eigenライブラリ使用、O3)  | 0.004  |
|C++(NTLライブラリ使用)|0.006|
|C++(NTLライブラリ使用、O3)|0.006|
|Python|0.274|
|SageMath|0.897|
|Risa/Asir|0.6324|

やっぱりCが一番早くて、Risa/AsirやSgaeMathがとても遅かった。また、参考までに```-O3```オプションを付けた場合も表の中に入れた。

# 最後に
ここで使用した各言語のソースコードはGitHub( https://github.com/satoshin-des/LLL_reduces )の方に公開してありますので、もし気になる方はそちらをのぞいてみてください。