| 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 | } |
|---|
