目次
1.C言語で実装
2.C++で実装
3.Pythonで実装
4.SageMathで実装
5.Risa/Asirで実装
6.実行してみた
7.比較してみた
8.GitHubにあげました
コメントでの指摘にもある通り、Cなどでヘッダに関数を直接書くのは危ないようです。そもそもこのなんちゃって実験においてわざわざヘッダに書く必要もないのでcファイル、cppファイルに直接関数を定義することにしました。
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++での実装
外部ライブラリ不使用版
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ライブラリ
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には基底行列が入っているものとします。
LLL簡約（NTL）
#include <NTL/ZZ.h>
#include <NTL/mat_ZZ.h>
#include <NTL/LLL.h>
NTL::ZZ _;
NTL::LLL(_, b, 99, 100);
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には基底行列が入っているものとします。
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$
それぞれの実行時間の比較
ここでは、ランダムに生成した完全階数の40次元格子基底をLLL基底簡約できるまでにかかった時間を比較します。
格子基底の生成方法
SageMathを用いて生成する：
GenBasis.sage
A = random_matrix(ZZ, 40, 40)
実際に使う格子基底
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言語での結果
$ 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++での結果
$ 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ライブラリ使用）
$ 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ライブラリ使用）
$ 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での結果
$ 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での結果
$ 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での結果
$ 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 )の方に公開してありますので、もし気になる方はそちらをのぞいてみてください。