-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathLLL.c
158 lines (136 loc) · 9.44 KB
/
LLL.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double dot_dbl_dbl(double *x, double *y, const int n);
double dot_int_dbl(int *x, double *y, const int n);
double dot_int_int(int *x, int *y, const int n);
void GSO(int **b, double *B, double **mu, const int n, const int m);
void SizeReduce(int **b, double **mu, const int i, const int j, const int m);
void LLLReduce(int **b, const double d, const int n, const int m);
void PrintMat(int **b, int n, int m){
for(int i = 0, j; i < n; ++i){
for(j = 0; j < m; ++j) printf("%d, ", b[i][j]);
puts("\b\b");
}
}
int main(){
int **b, t[][40] = {
{-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}
};
b = (int **)malloc(40 * sizeof(int *));
for(int i = 0, j; i < 40; ++i){
b[i] = (int *)malloc(40 * sizeof(int));
for(j = 0; j < 40; ++j) b[i][j] = t[i][j];
}
LLLReduce(b, 0.99, 40, 40);
PrintMat(b, 40, 40);
printf("norm = %lf\n", sqrt(dot_int_int(b[0], b[0], 40)));
return 0;
}
/* 内積 */
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;
}
double dot_int_int(int *x, int *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;
}
}