1 | /*--------------------------------------------------------------------------- |
---|
2 | / Integral Kernel using LLL Algorithm | |
---|
3 | / | |
---|
4 | / Implemented by | |
---|
5 | / Alfredo Sánchez-Navarro and Alberto Vigneron-Tenorio | |
---|
6 | / alfredo.sanchez@uca.es alberto.vigneron@uca.es, | |
---|
7 | / E.U.E.Empresariales de Jerez. | |
---|
8 | / Universidad de Cádiz | |
---|
9 | / Porvera 54. | |
---|
10 | / Jerez de la Frontera (Cádiz, Spain) | |
---|
11 | / | |
---|
12 | / | |
---|
13 | / This application computes the LLL-reduced basis of the Z-kernel of a | |
---|
14 | / matrix using the algorithm 2.7.2 appeared in: | |
---|
15 | / Cohen, H. A course in computational algebraic number theory. | |
---|
16 | / GTM 138, Springer 1996. | |
---|
17 | / | |
---|
18 | / The Input is a file with the structure: | |
---|
19 | / Number of coordinates, Number of vectors, | |
---|
20 | / b11 b12 b13 ... b1C | |
---|
21 | / b21 b22 b23 ... b2C | |
---|
22 | / ... | |
---|
23 | / bF1 bF2 bF3 ... bFC | |
---|
24 | / | |
---|
25 | / All these inputs are integer numbers. | |
---|
26 | / | |
---|
27 | / The Output is a LLL-reduced basis of the kernel of the matrix | |
---|
28 | / b11 b21 b31 ... bF1 | |
---|
29 | / b12 b22 b32 ... bF2 | |
---|
30 | / ... | |
---|
31 | / b1C b2C bC3 ... bFC | |
---|
32 | / | |
---|
33 | / | |
---|
34 | / All these outputs are integer numbers. | |
---|
35 | ---------------------------------------------------------------------------*/ |
---|
36 | |
---|
37 | |
---|
38 | /*--------------------------------------------------------------------------- |
---|
39 | | EXAMPLE: | |
---|
40 | | | |
---|
41 | | Let prf be the file: | |
---|
42 | | 4,3 | |
---|
43 | | 1 -3 2 | |
---|
44 | | 1 0 -1 | |
---|
45 | | 2 -4 2 | |
---|
46 | | 4 -7 3 | |
---|
47 | | then: | |
---|
48 | | | |
---|
49 | | ./ikernel <prf | |
---|
50 | | Output: | |
---|
51 | | -4 -2 3 0 | |
---|
52 | | -9 -5 5 1 | |
---|
53 | | | |
---|
54 | | The output is a LLL-reduced basis of the Z-kernel of ZxZxZxZ generated by | |
---|
55 | | (-4,-2,3,0), (-9,-5,5,1) | |
---|
56 | | | |
---|
57 | ---------------------------------------------------------------------------*/ |
---|
58 | |
---|
59 | |
---|
60 | #include <stdio.h> |
---|
61 | #include <stdlib.h> |
---|
62 | #include <string.h> |
---|
63 | #include <ctype.h> |
---|
64 | #include <new> |
---|
65 | #include <gmp.h> |
---|
66 | |
---|
67 | #define llln_max(a,b) (a > b ? a : b) |
---|
68 | |
---|
69 | /* Global variables */ |
---|
70 | mpz_t *llln_d; |
---|
71 | int *llln_f; |
---|
72 | mpz_t **llln_lambda; |
---|
73 | mpz_t **llln_b; |
---|
74 | mpz_t **llln_H; |
---|
75 | |
---|
76 | mpz_t **llln_proc_entrada (int *F, int *C); |
---|
77 | void llln_swapk (int k, int kmax, int C, int F); |
---|
78 | void llln_redi (int k, int l, int C, int F); |
---|
79 | void llln_test (int *k, int kma, int C, int n); |
---|
80 | void llln_producto (mpz_t *a, mpz_t *f, int n, mpz_t *prod); |
---|
81 | //int llln_max (int a, int b); |
---|
82 | |
---|
83 | void main (void) |
---|
84 | { |
---|
85 | int i, j, k, kmax, r; |
---|
86 | int F,C; |
---|
87 | mpz_t t; |
---|
88 | mpz_t prod; |
---|
89 | mpz_t u; |
---|
90 | mpz_t auxi1; |
---|
91 | mpz_t auxi2; |
---|
92 | |
---|
93 | mpz_init (prod); |
---|
94 | mpz_init (u); |
---|
95 | mpz_init (auxi1); |
---|
96 | mpz_init (auxi2); |
---|
97 | mpz_init (t); |
---|
98 | |
---|
99 | k = 2; |
---|
100 | kmax = 1; |
---|
101 | |
---|
102 | llln_b=llln_proc_entrada(&F,&C); |
---|
103 | |
---|
104 | llln_H=new (mpz_t *)[F+1]; |
---|
105 | for (i=1; i<=F; i++) { |
---|
106 | llln_H[i]=new mpz_t[F+1]; |
---|
107 | for (j=1; j<F+1; j++) { |
---|
108 | mpz_init (llln_H[i][j]); |
---|
109 | if (i==j) |
---|
110 | mpz_set_ui (llln_H[i][j],1); |
---|
111 | else |
---|
112 | mpz_set_ui (llln_H[i][j],0); |
---|
113 | } |
---|
114 | } |
---|
115 | |
---|
116 | |
---|
117 | llln_f=new int[F+1]; |
---|
118 | |
---|
119 | llln_d=new mpz_t[F+1]; |
---|
120 | llln_lambda=new (mpz_t *)[F+1]; |
---|
121 | mpz_init (llln_d[0]); |
---|
122 | for (i=1; i<=F; i++) { |
---|
123 | mpz_init (llln_d[i]); |
---|
124 | llln_lambda[i]=new mpz_t[F+1]; |
---|
125 | for (j=1; j<=F; j++) |
---|
126 | mpz_init(llln_lambda[i][j]); |
---|
127 | } |
---|
128 | |
---|
129 | mpz_set_ui (llln_d[0],1); |
---|
130 | llln_producto (llln_b[1],llln_b[1], C, &prod); |
---|
131 | mpz_set (t, prod); |
---|
132 | |
---|
133 | if (mpz_cmp_ui (t,0) != 0 ) { |
---|
134 | mpz_set (llln_d[1], t); |
---|
135 | llln_f[1]=1; |
---|
136 | } |
---|
137 | else { |
---|
138 | mpz_set_ui (llln_d[1],1); |
---|
139 | llln_f[1]=0; |
---|
140 | } |
---|
141 | |
---|
142 | while (k <= F) { |
---|
143 | if (k > kmax) { |
---|
144 | kmax = k; |
---|
145 | for (j=1;j<=k;j++) { |
---|
146 | if (llln_f[j]==0 & j<k) |
---|
147 | mpz_set_ui (llln_lambda[k][j], 0); |
---|
148 | else { |
---|
149 | llln_producto (llln_b[k],llln_b[j], C, &prod); |
---|
150 | mpz_set (u, prod); |
---|
151 | |
---|
152 | for (i=1; i<=j-1; i++) { |
---|
153 | if (llln_f[i] != 0) { |
---|
154 | mpz_mul (auxi1, llln_d[i], u); |
---|
155 | mpz_mul (auxi2, llln_lambda[k][i], llln_lambda[j][i]); |
---|
156 | mpz_sub (u, auxi1, auxi2); |
---|
157 | mpz_tdiv_q ( u, u, llln_d[i-1]); |
---|
158 | } |
---|
159 | } |
---|
160 | |
---|
161 | if (j<k) { |
---|
162 | mpz_set (llln_lambda[k][j], u); |
---|
163 | } |
---|
164 | else { |
---|
165 | if (mpz_sgn (u) == 0) { |
---|
166 | mpz_set (llln_d[k], llln_d[k-1]); |
---|
167 | llln_f[k]=0; |
---|
168 | } |
---|
169 | else { |
---|
170 | mpz_set (llln_d[k], u); |
---|
171 | llln_f[k]=1; |
---|
172 | } |
---|
173 | } |
---|
174 | } |
---|
175 | } |
---|
176 | } |
---|
177 | llln_test (&k, kmax, C, F); |
---|
178 | for (i=k-2; i>0; i--) { |
---|
179 | if (llln_f[i] != 0) |
---|
180 | llln_redi (k,i,C,F); |
---|
181 | } |
---|
182 | k = k+1; |
---|
183 | } |
---|
184 | |
---|
185 | |
---|
186 | r=0; |
---|
187 | i=1; |
---|
188 | while (llln_f[i] == 0 & i<=F) |
---|
189 | r = i++; |
---|
190 | |
---|
191 | for (i=1; i<=r; i++) { |
---|
192 | printf("\n"); |
---|
193 | for (j=1; j<=F; j++) { |
---|
194 | fprintf(stdout, " "); |
---|
195 | mpz_out_str (stdout, 0, llln_H[i][j]); |
---|
196 | fprintf(stdout, " "); |
---|
197 | } |
---|
198 | printf ("\n"); |
---|
199 | } |
---|
200 | mpz_clear (auxi1); |
---|
201 | mpz_clear (auxi2); |
---|
202 | mpz_clear (u); |
---|
203 | mpz_clear (prod); |
---|
204 | mpz_clear (t); |
---|
205 | |
---|
206 | for (i = 1; i <= F; i++) |
---|
207 | { |
---|
208 | mpz_clear (llln_d[i]); |
---|
209 | for (j = 1; j <= F; j++) |
---|
210 | mpz_clear (llln_lambda[i][j]); |
---|
211 | mpz_clear (llln_H[i][j]); |
---|
212 | delete[]llln_lambda[i]; |
---|
213 | delete[]llln_H[i]; |
---|
214 | } |
---|
215 | delete[]llln_lambda; |
---|
216 | delete[]llln_H; |
---|
217 | delete[]llln_d; |
---|
218 | delete[]llln_f; |
---|
219 | for (i = 1; i < F + 1; i++) |
---|
220 | { |
---|
221 | for (j = 1; j < C + 1; j++) |
---|
222 | mpz_clear (llln_b[i][j]); |
---|
223 | delete[]llln_b[i]; |
---|
224 | } |
---|
225 | delete[]llln_b; |
---|
226 | } |
---|
227 | |
---|
228 | /*--------------------------------------------------------------------------- |
---|
229 | | proc_entrada {}: | |
---|
230 | | Read the input file. | |
---|
231 | |___________________________________________________________________________|*/ |
---|
232 | |
---|
233 | mpz_t **llln_proc_entrada (int *F, int *C) /* Inputs */ |
---|
234 | { |
---|
235 | int numF, numC; |
---|
236 | int i, j, c; |
---|
237 | mpz_t **base; |
---|
238 | mpz_t valor; |
---|
239 | char *aux, caux; |
---|
240 | int tamano=100; |
---|
241 | |
---|
242 | aux=new char[tamano]; |
---|
243 | fscanf (stdin, "%i, %i\n", &numF, &numC); |
---|
244 | base = new (mpz_t *)[numF+1]; |
---|
245 | for (i = 1; i <= numF; i++) { |
---|
246 | base[i] = new mpz_t[numC+1]; |
---|
247 | for (j = 1; j <= numC; j++) { |
---|
248 | do { |
---|
249 | caux=getc (stdin); |
---|
250 | } |
---|
251 | while (!feof(stdin) && !(caux=='-' || isdigit(caux))); |
---|
252 | mpz_init (base[i][j]); |
---|
253 | c=0; |
---|
254 | if (!feof(stdin)) { |
---|
255 | do { |
---|
256 | aux[c++]=caux; |
---|
257 | if (c==tamano) { |
---|
258 | char *nuevo=new char[3*tamano/2]; |
---|
259 | memcpy (nuevo, aux, tamano); |
---|
260 | delete[] aux; |
---|
261 | aux=nuevo; |
---|
262 | tamano=3*tamano/2; |
---|
263 | } |
---|
264 | caux=getc(stdin); |
---|
265 | } |
---|
266 | while (!feof(stdin) && caux!=' ' && caux!='\n' && caux!='\t'); |
---|
267 | aux[c]='\0'; |
---|
268 | mpz_set_str (base[i][j], aux, 0); |
---|
269 | } |
---|
270 | else { |
---|
271 | printf("Empty file.\n"); |
---|
272 | mpz_set_ui (base[i][j],0); |
---|
273 | } |
---|
274 | } |
---|
275 | } |
---|
276 | *F = numF; |
---|
277 | *C = numC; |
---|
278 | delete[] aux; |
---|
279 | return (base); |
---|
280 | } |
---|
281 | |
---|
282 | /*--------------------------------------------------------------------------- |
---|
283 | | llln_producto {}: | |
---|
284 | | Multiply the integer vectors a and f with n coordinates. | |
---|
285 | |___________________________________________________________________________|*/ |
---|
286 | |
---|
287 | void llln_producto (mpz_t * a, mpz_t * f, int n, mpz_t * prod) |
---|
288 | { |
---|
289 | int i; |
---|
290 | mpz_t auxi; |
---|
291 | |
---|
292 | mpz_init (auxi); |
---|
293 | mpz_set_ui (*prod, 0); |
---|
294 | for (i = 1; i < n + 1; i++) { |
---|
295 | mpz_mul (auxi, a[i], f[i]); |
---|
296 | mpz_add (*prod, *prod, auxi); |
---|
297 | } |
---|
298 | mpz_clear (auxi); |
---|
299 | } |
---|
300 | |
---|
301 | |
---|
302 | |
---|
303 | void llln_test (int *k, int kma, int C, int n) |
---|
304 | { |
---|
305 | int i; |
---|
306 | |
---|
307 | if (llln_f[*k-1] != 0) |
---|
308 | llln_redi (*k, *k -1, C, n); |
---|
309 | |
---|
310 | if (llln_f[*k-1] != 0 & llln_f[*k] == 0) { |
---|
311 | llln_swapk (*k, kma, C, n); |
---|
312 | *k = llln_max (2, *k-1); |
---|
313 | llln_test (k, kma, C, n); |
---|
314 | } |
---|
315 | } |
---|
316 | |
---|
317 | |
---|
318 | /*--------------------------------------------------------------------------- |
---|
319 | | lll_redi {}: | |
---|
320 | | Make some changes in global variables. | |
---|
321 | |___________________________________________________________________________|*/ |
---|
322 | |
---|
323 | |
---|
324 | void llln_redi (int k, int l, int C, int F) |
---|
325 | { |
---|
326 | int i; |
---|
327 | mpz_t q; |
---|
328 | mpz_t auxi1; |
---|
329 | |
---|
330 | mpz_init (q); |
---|
331 | mpz_init (auxi1); |
---|
332 | |
---|
333 | mpz_mul_ui (auxi1, llln_lambda[k][l], 2); |
---|
334 | mpz_abs (auxi1, auxi1); |
---|
335 | |
---|
336 | if (mpz_cmp (auxi1, llln_d[l]) > 0) { |
---|
337 | mpz_mul_ui (auxi1, llln_lambda[k][l], 2); |
---|
338 | mpz_add (auxi1, auxi1, llln_d[l]); |
---|
339 | mpz_mul_ui (q, llln_d[l], 2); |
---|
340 | mpz_fdiv_q (q, auxi1, q); |
---|
341 | |
---|
342 | for (i=1; i<F+1; i++) { |
---|
343 | mpz_mul (auxi1, q, llln_H[l][i]); |
---|
344 | mpz_sub (llln_H[k][i], llln_H[k][i], auxi1); |
---|
345 | } |
---|
346 | for (i=1; i<C+1; i++) { |
---|
347 | mpz_mul (auxi1, q, llln_b[l][i]); |
---|
348 | mpz_sub (llln_b[k][i], llln_b[k][i], auxi1); |
---|
349 | } |
---|
350 | mpz_mul (auxi1, q, llln_d[l]); |
---|
351 | |
---|
352 | mpz_sub (llln_lambda[k][l], llln_lambda[k][l], auxi1); |
---|
353 | |
---|
354 | for (i=1; i<=l-1; i++){ |
---|
355 | mpz_mul (auxi1, q, llln_lambda[l][i]); |
---|
356 | mpz_sub (llln_lambda[k][i], llln_lambda[k][i], auxi1); |
---|
357 | } |
---|
358 | } |
---|
359 | |
---|
360 | mpz_clear (q); |
---|
361 | mpz_clear (auxi1); |
---|
362 | } |
---|
363 | |
---|
364 | |
---|
365 | /*--------------------------------------------------------------------------- |
---|
366 | | lll_swapi {}: | |
---|
367 | | Exchange some global variables | |
---|
368 | |___________________________________________________________________________|*/ |
---|
369 | |
---|
370 | |
---|
371 | void llln_swapk (int k, int kmax, int C, int F) |
---|
372 | { |
---|
373 | int i,j; |
---|
374 | mpz_t lamb; |
---|
375 | mpz_t auxi1; |
---|
376 | mpz_t t; |
---|
377 | mpz_t swap; |
---|
378 | |
---|
379 | mpz_init (lamb); |
---|
380 | mpz_init (auxi1); |
---|
381 | mpz_init (t); |
---|
382 | mpz_init (swap); |
---|
383 | |
---|
384 | for (j=1; j<F+1 ;j++) { |
---|
385 | mpz_set (swap, llln_H[k][j]); |
---|
386 | mpz_set (llln_H[k][j], llln_H[k-1][j]); |
---|
387 | mpz_set (llln_H[k-1][j], swap); |
---|
388 | } |
---|
389 | |
---|
390 | |
---|
391 | for (j=1; j<C+1 ;j++) { |
---|
392 | mpz_set (swap, llln_b[k][j]); |
---|
393 | mpz_set (llln_b[k][j], llln_b[k-1][j]); |
---|
394 | mpz_set (llln_b[k-1][j], swap); |
---|
395 | } |
---|
396 | |
---|
397 | if (k>2) { |
---|
398 | for (j=1;j<=k-2;j++) { |
---|
399 | mpz_set (swap, llln_lambda[k][j]); |
---|
400 | mpz_set (llln_lambda[k][j], llln_lambda[k-1][j]); |
---|
401 | mpz_set (llln_lambda[k-1][j], swap); |
---|
402 | } |
---|
403 | } |
---|
404 | mpz_set (lamb, llln_lambda[k][k-1]); |
---|
405 | if (mpz_sgn (lamb) == 0) { |
---|
406 | mpz_set (llln_d[k-1], llln_d[k-2]); |
---|
407 | llln_f[k-1]= 0; |
---|
408 | llln_f[k]= 1; |
---|
409 | mpz_set_ui (llln_lambda[k][k-1], 0); |
---|
410 | |
---|
411 | for (j=k+1; j <= kmax ; j++) { |
---|
412 | mpz_set (llln_lambda[j][k], llln_lambda[j][k-1]); |
---|
413 | mpz_set_ui (llln_lambda[j][k-1],0); |
---|
414 | } |
---|
415 | } |
---|
416 | else |
---|
417 | { |
---|
418 | for (j=k+1; j <= kmax ; j++)s |
---|
419 | { |
---|
420 | mpz_mul (auxi1, lamb, llln_lambda[j][k-1]); |
---|
421 | mpz_tdiv_q (llln_lambda[j][k-1], auxi1, llln_d[k-1]); |
---|
422 | } |
---|
423 | mpz_set (t, llln_d[k]); |
---|
424 | mpz_mul (auxi1, lamb, lamb); |
---|
425 | mpz_tdiv_q (llln_d[k-1], auxi1, llln_d[k-1]); |
---|
426 | mpz_set (llln_d[k], llln_d[k-1]); |
---|
427 | |
---|
428 | for (j=k+1; j <= kmax -1 ; j++) |
---|
429 | { |
---|
430 | for (i=j+1; i<= kmax; i++) |
---|
431 | { |
---|
432 | mpz_mul (auxi1, llln_lambda[i][j], llln_d[k-1]); |
---|
433 | mpz_tdiv_q (llln_lambda[i][j], auxi1, t); |
---|
434 | } |
---|
435 | } |
---|
436 | |
---|
437 | for (j=k+1; j<= kmax; j++) |
---|
438 | { |
---|
439 | mpz_mul (auxi1, llln_d[j], llln_d[k-1]); |
---|
440 | mpz_tdiv_q (llln_d[j], auxi1, t); |
---|
441 | } |
---|
442 | } |
---|
443 | |
---|
444 | mpz_clear (auxi1); |
---|
445 | mpz_clear (lamb); |
---|
446 | mpz_clear (t); |
---|
447 | mpz_clear (swap); |
---|
448 | } |
---|