-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathLLL.rr
128 lines (117 loc) · 7.6 KB
/
LLL.rr
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
/* 内積 */
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;
}
def main(){
B = matrix(40, 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]
]);
tstart();
B = lll(B, 0.99);
print(B);
print(deval(dot(B[0], B[0]) ^ (1 / 2)));
tstop();
}
end$