source: git/IntegerProgramming/LLL_lattice_to_send.cc @ 543931

spielwiese
Last change on this file since 543931 was 543931, checked in by Hans Schönemann <hannes@…>, 18 years ago
*hannes: sparc-cc git-svn-id: file:///usr/local/Singular/svn/trunk@9229 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 11.9 KB
Line 
1/*---------------------------------------------------------------------------
2/                         Integral LLL Algorithm                            |
3/                                                                           |
4/                           Implemented by                                  |
5/          Alberto Vigneron-Tenorio and Alfredo Sánchez-Navarro             |
6/            alberto.vigneron@uca.es, alfredo.sanchez@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 a Z-lattice using      |
14/ the algorithm 2.6.7 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 vectors, Number of coordinates                                  |
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 Z-lattice generated by           |
28/ (b11,b12,...,b1C),...,(bF1,...,bFC)                                       |
29/                                                                           |
30/ All these outputs are integer numbers.                                    |
31---------------------------------------------------------------------------*/
32
33
34/*---------------------------------------------------------------------------
35| EXAMPLE:                                                                  |
36|                                                                           |
37| Let prf be the file:                                                      |
38| 2,3                                                                       |
39| 1 -2 4                                                                    |
40| 2 1 -7                                                                    |
41| then:                                                                     |
42|                                                                           |
43| ./LLL_lattice <prf                                                        |
44| Output:                                                                   |
45|  1  -2  4                                                                 |
46|  3  -1  -3                                                                |
47|                                                                           |
48| The output is a LLL-reduced basis of the lattice of ZxZxZ generated by    |
49| (1,-2,4),(2,1,-7)                                                         |
50|                                                                           |
51---------------------------------------------------------------------------*/
52
53
54#include <stdio.h>
55#include <stdlib.h>
56#include <string.h>
57#include <ctype.h>
58#include <new>
59#include <cstring>
60#include <si_gmp.h>
61
62#define LLL_MAX(a,b) (a > b ? a : b)
63
64/* Global variables */
65mpz_t * lll_d;
66mpz_t ** lll_lambda;
67mpz_t ** lll_b;
68
69mpz_t **proc_entrada (int *F, int *C);
70void lll_swapi (int k, int kmax, int C);
71void lll_redi (int k, int l, int C);
72void lll_test (int *k, int kma, int C);
73void lll_product (mpz_t * a, mpz_t * f, int n, mpz_t * prod);
74
75int main (void)
76{
77  int i, j, k, kmax;
78  int F, C;
79  mpz_t prod;
80  mpz_t u;
81  mpz_t auxi1;
82  mpz_t auxi2;
83  mpz_init (prod);
84  mpz_init (u);
85  mpz_init (auxi1);
86  mpz_init (auxi2);
87
88  k = 2;    /* Initialize */
89  kmax = 1; /* Initialize */
90
91  lll_b = proc_entrada (&F, &C);
92
93  lll_d = new mpz_t[F + 1];
94  lll_lambda = new mpz_t *[F + 1];
95  mpz_init (lll_d[0]);
96  for (i = 1; i <= F; i++) {
97    mpz_init (lll_d[i]);
98    lll_lambda[i] = new mpz_t[F + 1];
99    for (j = 1; j <= F; j++)
100      mpz_init (lll_lambda[i][j]);
101  }
102  mpz_set_ui (lll_d[0], 1); /* Initialize */
103  lll_product (lll_b[1], lll_b[1], C, &prod);
104  mpz_set (lll_d[1], prod); /* Initialize */
105
106
107/* Incremental Gram-Schmidt */
108  while (k <= F)
109  {
110    if (k > kmax)
111    {
112      kmax = k;
113      for (j = 1; j < k + 1; j++)
114      {
115        lll_product (lll_b[k], lll_b[j], C, &prod);
116        mpz_set (u, prod);
117        for (i = 1; i < j - 1 + 1; i++)
118        {
119          mpz_mul (auxi1, lll_d[i], u);
120          mpz_mul (auxi2, lll_lambda[k][i], lll_lambda[j][i]);
121          mpz_sub (u, auxi1, auxi2);
122          mpz_tdiv_q (u, u, lll_d[i - 1]);
123        }
124        if (j < k)
125        {
126          mpz_set (lll_lambda[k][j], u);
127        }
128        else
129        {
130          mpz_set (lll_d[k], u);
131          if (mpz_cmp_ui (lll_d[k], 0) == 0)
132          { /* The input is not a basis */
133            fprintf (stderr, "\n INPUT NOT BASIS !!!! \n");
134            exit (1);
135          }
136        }
137      }
138    }
139    lll_test (&k, kmax, C);
140    for (i = k - 2; i > 0; i--)
141      lll_redi (k, i, C);
142    k = k + 1;
143  }
144  fprintf (stdout, "list res=");
145  for (i = 1; i < F + 1; i++)
146  {
147    printf ("\nlist(");
148    for (j = 1; j < C + 1; j++)
149    {
150      fprintf (stdout, " ");
151      mpz_out_str (stdout, 0, lll_b[i][j]);
152      if (j!=C) fprintf (stdout, ",");
153    }
154    if (i!=F) printf("),\n");
155    else      printf(");\n");
156  }
157  mpz_clear (auxi1);
158  mpz_clear (auxi2);
159  mpz_clear (u);
160  mpz_clear (prod);
161  mpz_clear (lll_d[0]);
162  for (i = 1; i <= F; i++)
163  {
164    mpz_clear (lll_d[i]);
165    for (j = 1; j <= F; j++)
166      mpz_clear (lll_lambda[i][j]);
167    delete[]lll_lambda[i];
168  }
169  delete[]lll_lambda;
170  delete[]lll_d;
171  for (i = 1; i < F + 1; i++)
172  {
173    for (j = 1; j < C + 1; j++)
174       mpz_clear (lll_b[i][j]);
175    delete[]lll_b[i];
176  }
177  delete[]lll_b;
178  return 0;
179}
180
181/*---------------------------------------------------------------------------
182| proc_entrada {}:                                                          |
183| Read the input file.                                                      |
184|___________________________________________________________________________|*/
185
186mpz_t **proc_entrada (int *F, int *C) /* Inputs */
187{
188  int numF, numC;
189  int i, j, c;
190  mpz_t **base;
191  mpz_t valor;
192  char *aux, caux;
193  int tamano=100;
194
195  aux=new char[tamano];
196  fscanf (stdin, "%i, %i\n", &numF, &numC);
197  base = new mpz_t *[numF+1];
198  for (i = 1; i <= numF; i++) {
199    base[i] = new mpz_t[numC+1];
200    for (j = 1; j <= numC; j++)
201    {
202      do { caux=getc (stdin); }
203      while (!feof(stdin) && !(caux=='-' || isdigit(caux)));
204      mpz_init (base[i][j]);
205      c=0;
206      if (!feof(stdin))
207      {
208        do
209        {
210          aux[c++]=caux;
211          if (c==tamano)
212          {
213            char *nuevo=new char[3*tamano/2];
214            memcpy (nuevo, aux, tamano);
215            delete[] aux;
216            aux=nuevo;
217            tamano=3*tamano/2;
218          }
219          caux=getc(stdin);
220        }
221        while (!feof(stdin) && caux!=' ' && caux!='\n' && caux!='\t');
222        aux[c]='\0';
223        mpz_set_str (base[i][j], aux, 0);
224      }
225      else
226      {
227        fprintf(stderr,"Empty file.\n");
228        mpz_set_ui (base[i][j],0);
229      }
230    }
231  }
232  *F = numF;
233  *C = numC;
234  delete[] aux;
235  return (base);
236}
237
238/*---------------------------------------------------------------------------
239| lll_product {}:                                                           |
240| Multiply the integer vectors a and f with n coordinates.                  |
241|___________________________________________________________________________|*/
242
243void  lll_product (mpz_t * a, mpz_t * f, int n, mpz_t * prod)
244{
245  int i;
246  mpz_t auxi;
247
248  mpz_init (auxi);
249  mpz_set_ui (*prod, 0);
250  for (i = 1; i < n + 1; i++)
251  {
252    mpz_mul (auxi, a[i], f[i]);
253    mpz_add (*prod, *prod, auxi);
254  }
255  mpz_clear (auxi);
256}
257
258void lll_test (int *k, int kma, int C) /* Test LLL condition */
259{
260  int i;
261  mpz_t auxi1;
262  mpz_t auxi2;
263  mpz_t auxi3;
264  mpz_t auxi4;
265
266  mpz_init (auxi1);
267  mpz_init (auxi2);
268  mpz_init (auxi3);
269
270  lll_redi (*k, *k - 1, C);
271  mpz_mul (auxi1, lll_d[*k], lll_d[*k - 2]);
272  mpz_mul_ui (auxi1, auxi1, 4);
273  mpz_mul (auxi2, lll_d[*k - 1], lll_d[*k - 1]);
274  mpz_mul_ui (auxi2, auxi2, 3);
275  mpz_mul (auxi3, lll_lambda[*k][*k - 1], lll_lambda[*k][*k - 1]);
276  mpz_mul_ui (auxi3, auxi3, 4);
277  mpz_sub (auxi3, auxi2, auxi3);
278  if (mpz_cmp (auxi1, auxi3) < 0)
279  {
280    lll_swapi (*k, kma, C);
281    *k = LLL_MAX (2, *k - 1);
282    lll_test (&*k, kma, C);
283  }
284
285  mpz_clear (auxi1);
286  mpz_clear (auxi2);
287  mpz_clear (auxi3);
288}
289
290/*---------------------------------------------------------------------------
291| lll_redi {}:                                                              |
292| Make some changes in global variables.                                    |
293|___________________________________________________________________________|*/
294
295
296void lll_redi (int k, int l, int C)
297{
298  int i;
299  mpz_t q;
300  mpz_t auxi1;
301
302  mpz_init (q);
303  mpz_init (auxi1);
304
305  mpz_mul_ui (auxi1, lll_lambda[k][l], 2);
306  mpz_abs (auxi1, auxi1);
307
308    if (mpz_cmp (auxi1, lll_d[l]) > 0)
309    {
310      mpz_mul_ui (auxi1, lll_lambda[k][l], 2);
311      mpz_add (auxi1, auxi1, lll_d[l]);
312      mpz_mul_ui (q, lll_d[l], 2);
313      mpz_fdiv_q (q, auxi1, q);
314      for (i = 1; i < C + 1; i++)
315      {
316        mpz_mul (auxi1, q, lll_b[l][i]);
317        mpz_sub (lll_b[k][i], lll_b[k][i], auxi1);
318      }
319      mpz_mul (auxi1, q, lll_d[l]);
320      mpz_sub (lll_lambda[k][l], lll_lambda[k][l], auxi1);
321      for (i = 1; i <= l - 1; i++)
322      {
323        mpz_mul (auxi1, q, lll_lambda[l][i]);
324        mpz_sub (lll_lambda[k][i], lll_lambda[k][i], auxi1);
325      }
326    }
327    mpz_clear (q);
328    mpz_clear (auxi1);
329}
330
331/*---------------------------------------------------------------------------
332| lll_swapi {}:                                                             |
333| Exchange some global variables                                            |
334|___________________________________________________________________________|*/
335
336void lll_swapi (int k, int kmax, int C)
337{
338  int j;
339  mpz_t lamb;
340  mpz_t B;
341  mpz_t auxi1;
342  mpz_t auxi2;
343  mpz_t t;
344  mpz_t swap;
345
346  mpz_init (lamb);
347  mpz_init (B);
348  mpz_init (auxi1);
349  mpz_init (auxi2);
350  mpz_init (t);
351  mpz_init (swap);
352
353  for (j = 1; j < C + 1; j++)
354  {
355    mpz_set (swap, lll_b[k][j]);
356    mpz_set (lll_b[k][j], lll_b[k - 1][j]);
357    mpz_set (lll_b[k - 1][j], swap);
358  }
359  if (k > 2)
360  {
361    for (j = 1; j <= k - 2; j++)
362    {
363      mpz_set (swap, lll_lambda[k][j]);
364      mpz_set (lll_lambda[k][j], lll_lambda[k - 1][j]);
365      mpz_set (lll_lambda[k - 1][j], swap);
366    }
367  }
368
369  mpz_set (lamb, lll_lambda[k][k - 1]);
370  mpz_mul (auxi1, lll_d[k - 2], lll_d[k]);
371  mpz_mul (auxi2, lamb, lamb);
372  mpz_add (B, auxi1, auxi2);
373  mpz_tdiv_q (B, B, lll_d[k - 1]);
374  for (j = k + 1; j <= kmax; j++)
375  {
376    mpz_set (t, lll_lambda[j][k]);
377    mpz_mul (auxi1, lll_d[k], lll_lambda[j][k - 1]);
378    mpz_mul (auxi2, lamb, t);
379    mpz_sub (auxi1, auxi1, auxi2);
380    mpz_tdiv_q (lll_lambda[j][k], auxi1, lll_d[k - 1]);
381    mpz_mul (auxi2, lamb, lll_lambda[j][k]);
382    mpz_mul (auxi1, B, t);
383    mpz_add (auxi1, auxi1, auxi2);
384    mpz_tdiv_q (lll_lambda[j][k - 1], auxi1, lll_d[k]);
385  }
386  mpz_set (lll_d[k - 1], B);
387
388  mpz_clear (auxi1);
389  mpz_clear (auxi2);
390  mpz_clear (lamb);
391  mpz_clear (B);
392  mpz_clear (t);
393  mpz_clear (swap);
394}
Note: See TracBrowser for help on using the repository browser.