1 /*
2 This file is part of PolyLib.
3
4 PolyLib is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
8
9 PolyLib is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with PolyLib. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18 #include <stdlib.h>
19 #include <polylib/polylib.h>
20
21 /* nota bene: on stocke les matrices par lignes */
22 /* nota bene: matrices are stored in row major order */
23
24 /*------------------------------------------------------------------------
25 change les signes de la ligne i de la matrice A
26 change the sign of row i of matrix A
27 ------------------------------------------------------------------------*/
28
moins_l(Value * a,int i,int n,int p)29 static void moins_l(Value *a,int i,int n,int p) {
30
31 int k;
32 Value *c;
33
34 c=a+(i-1)*p;
35
36 for(k=1; k<=p; k++) {
37 value_oppose(*c,*c);
38 c++;
39 }
40 return;
41 } /* moins_l */
42
43 /*------------------------------------------------------------------------
44 change les signes de la colonne i de la matrice A
45 change the sign of column i of matrix A
46 ------------------------------------------------------------------------*/
47
moins_c(Value * a,int i,int n,int p)48 static void moins_c(Value *a,int i,int n,int p) {
49
50 int k;
51 Value *c;
52
53 c=a+(i-1);
54
55 for(k=1;k<=n;k++) {
56 value_oppose(*c,*c);
57 c=c+p;
58 }
59 return;
60 } /* moins_c */
61
62 /*------------------------------------------------------------------------
63 echange les lignes i et j de la matrice A
64 exchange row i and j of matrix A
65 ------------------------------------------------------------------------*/
66
echange_l(Value * a,int i,int j,int n,int p)67 static void echange_l(Value *a,int i,int j,int n,int p) {
68
69 int k;
70 Value s;
71 Value *c1,*c2;
72
73 value_init(s);
74 c1=a+(i-1)*p;
75 c2=a+(j-1)*p;
76
77 for(k=1;k<=p;k++) {
78 value_assign(s,*c1);
79 value_assign(*c1,*c2);
80 value_assign(*c2,s);
81 c1++;
82 c2++;
83 }
84 value_clear(s);
85 return;
86 } /* echange_l */
87
88 /*------------------------------------------------------------------------
89 echange les colonnes i et j de la matrice A
90 exchange columns i and j of matrix A
91 ------------------------------------------------------------------------*/
92
echange_c(Value * a,int i,int j,int n,int p)93 static void echange_c(Value *a,int i,int j,int n,int p) {
94
95 int k;
96 Value s;
97 Value *c1,*c2;
98
99 value_init(s);
100 c1=a+(i-1);
101 c2=a+(j-1);
102
103 for(k=1;k<=n;k++) {
104 value_assign(s,*c1);
105 value_assign(*c1,*c2);
106 value_assign(*c2,s);
107 c1=c1+p;
108 c2=c2+p;
109 }
110 value_clear(s);
111 return;
112 } /* echange_c */
113
114 /*------------------------------------------------------------------------
115 A est une matrice de taille n*p
116 on ajoute a la ligne i, x fois la ligne j
117 A is a matrix of size n*p
118 we add x times the row j to the row i
119 ------------------------------------------------------------------------*/
120
ligne(Value * a,int i,int j,Value x,int n,int p)121 static void ligne(Value *a,int i,int j,Value x,int n,int p) {
122
123 int k;
124 Value *c1,*c2;
125
126 c1=a+(i-1)*p;
127 c2=a+(j-1)*p;
128
129 for(k=1;k<=p;k++) {
130
131 value_addmul(*c1, x, *c2);
132 c1++;
133 c2++;
134 }
135 return;
136 } /* ligne */
137
138 /*------------------------------------------------------------------------
139 A est une matrice de taille n*p
140 on ajoute a la colonne i, x fois la colonne j
141 A is a matrix of size n*p
142 we add x times the column j to the column i
143
144 ------------------------------------------------------------------------*/
145
colonne(Value * a,int i,int j,Value x,int n,int p)146 static void colonne(Value *a,int i,int j,Value x,int n,int p) {
147
148 int k;
149 Value *c1,*c2;
150
151 c1=a+(i-1);
152 c2=a+(j-1);
153
154 for(k=1;k<=n;k++) {
155 value_addmul(*c1, x, *c2);
156 c1=c1+p;
157 c2=c2+p;
158 }
159 return;
160 } /* colonne */
161
162 /*----------------------------------------------------------------------
163 trouve le numero de colonne du plus petit element non nul de
164 la ligne q, d'indice superieur a q, de la matrice A, mais
165 renvoie zero si tous les elements sont nuls sauf le qieme.
166
167 find the column number (greater than q) of the smallest non
168 null element of row q . it returns
169 0 if all the last elements of row q are null except the qth
170
171 ----------------------------------------------------------------------*/
172
petit_l(Value * a,int n,int p,int q)173 static int petit_l(Value *a,int n,int p,int q) {
174
175 int numero=0, k, tousnuls;
176 Value minus, comp;
177 Value *c;
178
179 value_init(minus); value_init(comp);
180 c=a+(q-1)*p+(p-1);
181 tousnuls=1;
182 for(k=p;k>q;k--) {
183 value_absolute(comp,*c);
184 if (value_notzero_p(comp)) {
185 if (tousnuls==1) {
186 value_assign(minus,comp);
187 numero=k;
188 tousnuls=0;
189 }
190 else if (value_ge(minus,comp)) {
191 value_assign(minus,comp);
192 numero=k;
193 }
194 }
195 c--;
196 }
197 if (tousnuls==1) {
198 value_clear(minus); value_clear(comp);
199 return(0);
200 }
201 else {
202 value_absolute(comp,*c);
203 if ((value_notzero_p(comp))&&(value_ge(minus,comp))) {
204 value_clear(minus); value_clear(comp);
205 return(q);
206 }
207 else {
208 value_clear(minus); value_clear(comp);
209 return(numero);
210 }
211 }
212 } /* petit_l */
213
214 /*----------------------------------------------------------------------
215 trouve le numero de ligne du plus petit element non nul de
216 la colonne q, d'indice superieur a q, de la matrice A, mais
217 renvoie zero si tous les elements sont nuls sauf le qieme.
218
219 find the row number (greater than q) of the smallest non
220 null element of column q . it returns
221 0 if all the last elements of column q are null except the qth
222
223 ----------------------------------------------------------------------*/
224
petit_c(Value * a,int n,int p,int q)225 static int petit_c(Value *a,int n,int p,int q) {
226
227 int numero=0, k, tousnuls;
228 Value minus, comp;
229 Value *c;
230
231 value_init(minus); value_init(comp);
232 c = a+q-1+p*(n-1);
233 tousnuls=1;
234 for(k=n;k>q;k--) {
235 value_absolute(comp,*c);
236 if (value_notzero_p(comp)) {
237 if (tousnuls==1) {
238 value_assign(minus,comp);
239 numero=k;
240 tousnuls=0;
241 }
242 else if (value_ge(minus,comp)) {
243 value_assign(minus,comp);
244 numero=k;
245 }
246 }
247 c=c-p;
248 }
249
250 if (tousnuls==1) {
251 value_clear(minus); value_clear(comp);
252 return(0);
253 }
254 else {
255 value_absolute(comp,*c);
256 if ((value_notzero_p(comp)) && (value_ge(minus,comp))) {
257 value_clear(minus); value_clear(comp);
258 return(q);
259 }
260 else {
261 value_clear(minus); value_clear(comp);
262 return(numero);
263 }
264 }
265 } /* petit_c */
266
267 /*-----------------------------------------------------------------------
268 A est initialisee a une matrice "identite" de taille n*p
269 A is initialized to an "identity" matrix of size n*p
270 -----------------------------------------------------------------------*/
271
identite(Value * a,int n,int p)272 static void identite(Value *a,int n,int p) {
273
274 int i,j;
275 Value *b;
276
277 b = a;
278 for(i=1;i<=n;i++) {
279 for(j=1;j<=p;j++) {
280 if (i==j)
281 value_set_si(*b,1);
282 else
283 value_set_si(*b,0);
284 b++;
285 }
286 }
287 return;
288 } /* identite */
289
290 /*-----------------------------------------------------------------------
291 transpose la sous-matrice de A dont les indices sont
292 superieurs ou egaux a q
293 transposition of the sub-matrix of A composed of the elements
294 whose indices are greater or equal to q
295 -----------------------------------------------------------------------*/
296
transpose(Value * a,int n,int q)297 static void transpose(Value *a,int n,int q) {
298
299 int i;
300 Value val;
301 Value *b,*c;
302
303 value_init(val);
304 if (q<n) {
305 b=a+(q-1)+(q-1)*n;
306 c=b;
307 for(i=q+1;i<=n;i++) {
308 b++;
309 c=c+n;
310 value_assign(val,*b);
311 value_assign(*b,*c);
312 value_assign(*c,val);
313 }
314 transpose(a,n,q+1);
315 }
316 value_clear(val);
317 return;
318 } /* transpose */
319
320 /*----------------------------------------------------------------------
321 trouve le numero de colonne du premier element non divisible
322 par val dans la sous-matrice d'indices > q mais renvoie 0 sinon
323
324 find the column number of the first non-divisible by val element
325 in the sub-matrix of a composed of the elements of indices >q.
326 Return 0 is there is no such element
327 ----------------------------------------------------------------------*/
328
encore(Value * a,int n,int p,int q,Value val)329 static int encore(Value *a,int n,int p,int q,Value val) {
330
331 int k,l;
332 Value comp, tmp;
333 Value *c;
334
335 if ((q>=n)||(q>=p)) return(0);
336
337 value_init(comp); value_init(tmp);
338 c=a+q*p+q;
339 k=q+1;
340 l=q+1;
341 while (k<=n) {
342 value_absolute(comp,*c);
343 if (value_zero_p(val)) {
344 if (value_notzero_p(comp)) {
345 value_clear(comp);
346 value_clear(tmp);
347 return(l);
348 }
349 }
350 else {
351 value_modulus(tmp,comp,val);
352 if (value_notzero_p(tmp)) {
353 value_clear(comp);
354 value_clear(tmp);
355 return(l);
356 }
357 }
358 if (l==p) {
359 k=k+1;
360 l=q+1;
361 c=c+q+1;
362 }
363 else {
364 l++;
365 c++;
366 }
367 }
368 value_clear(comp);
369 value_clear(tmp);
370 return(0);
371 } /* encore */
372
373 /*----------------------------------------------------------------------
374 transforme la matrice A (de taille n*p), a partir de la position
375 q*q en sa forme de smith, les matrices b et c subissant les memes
376 transformations respectivement sur les lignes et colonnes
377
378 transform the matrix A (size n*p), from position (q,q) into
379 its smith normal form. the same transformations are applied to
380 matrices b and c, respectively on rows and on columns
381 ----------------------------------------------------------------------*/
382
smith(Value * a,Value * b,Value * c,Value * b_inverse,Value * c_inverse,int n,int p,int q)383 static void smith(Value *a,Value *b,Value *c,Value *b_inverse,Value *c_inverse,int n,int p,int q) {
384
385 int i,k,fini;
386 Value x, pivot, tmp, x_inv;
387 Value *f;
388
389 value_init(pivot); value_init(tmp);
390 value_init(x); value_init(x_inv);
391
392 if ((q<=n)&&(q<=p)) {
393 fini = 1;
394
395 while (fini!=0) {
396 i=petit_c(a,n,p,q);
397 while (i!=0) {
398 if (i!=q) {
399 echange_l(a,i,q,n,p);
400 echange_l(b,i,q,n,n);
401 echange_c(b_inverse,i,q,n,n);
402 }
403 f=a+(q-1)+(q-1)*p;
404 value_assign(pivot,*f);
405 if (value_neg_p(pivot)) {
406 moins_l(a,q,n,p);
407 moins_l(b,q,n,n);
408 moins_c(b_inverse,q,n,n);
409 value_oppose(pivot,pivot);
410 }
411 for(k=q+1;k<=n;k++) {
412 f=f+p;
413 if (value_notzero_p(*f)) {
414 value_division(x,*f,pivot);
415 value_modulus(tmp,*f,pivot);
416 if (value_neg_p(tmp))
417 value_decrement(x,x);
418 value_oppose(x_inv,x);
419 ligne(a,k,q,x_inv,n,p);
420 ligne(b,k,q,x_inv,n,n);
421 colonne(b_inverse,q,k,x,n,n);
422 }
423 }
424 i=petit_c(a,n,p,q);
425 }
426 fini=0;
427 i=petit_l(a,n,p,q);
428 while (i!=0) {
429 if (i!=q) {
430 echange_c(a,i,q,n,p);
431 echange_c(c,i,q,p,p);
432 echange_l(c_inverse,i,q,p,p);
433 fini=1;
434 }
435 f=a+(q-1)+(q-1)*p;
436 value_assign(pivot,*f);
437 if (value_neg_p(pivot)) {
438 moins_c(a,q,n,p);
439 moins_c(c,q,p,p);
440 moins_l(c_inverse,q,p,p);
441 value_oppose(pivot,pivot);
442 }
443 for(k=q+1;k<=p;k++) {
444 f++;
445 if (value_notzero_p(*f)) {
446 value_division(x,*f,pivot);
447 value_modulus(tmp,*f,pivot);
448 if (value_neg_p(tmp))
449 value_decrement(x,x);
450 value_oppose(x_inv,x);
451 colonne(a,k,q,x_inv,n,p);
452 colonne(c,k,q,x_inv,p,p);
453 ligne(c_inverse,q,k,x,p,p);
454 }
455 }
456 i=petit_l(a,n,p,q);
457 }
458 }
459 value_assign(pivot,*(a+(q-1)+(q-1)*p));
460 if (value_neg_p(pivot)) {
461 moins_l(a,q,n,p);
462 moins_l(b,q,n,n);
463 moins_c(b_inverse,q,n,n);
464 value_oppose(pivot,pivot);
465 }
466
467 i=encore(a,n,p,q,pivot);
468 if (i!=0) {
469 value_set_si(x,1);
470 value_set_si(x_inv,-1);
471 colonne(a,q,i,x,n,p);
472 colonne(c,q,i,x,p,p);
473 ligne(c_inverse,i,q,x_inv,p,p);
474 smith(a,b,c,b_inverse,c_inverse,n,p,q);
475 }
476 else smith(a,b,c,b_inverse,c_inverse,n,p,q+1);
477 }
478 value_clear(pivot); value_clear(tmp);
479 value_clear(x); value_clear(x_inv);
480
481 return;
482 } /* smith */
483
484 /*----------------------------------------------------------------------
485 decompose la matrice A en le produit d'une matrice B
486 unimodulaire et d'une matrice triangulaire superieure
487 D est l'inverse de B
488 Decompose the matrix A in the product of a matrix B unimodular
489 and a matrix upper triangular, D is the inverse of B
490 ----------------------------------------------------------------------*/
491
hermite(Value * a,Value * b,Value * d,int n,int p,int q)492 static void hermite(Value *a,Value *b,Value *d,int n,int p,int q) {
493
494 int i,k;
495 Value x, pivot, x_inv, tmp;
496 Value *c1;
497
498 value_init(pivot); value_init(tmp);
499 value_init(x); value_init(x_inv);
500
501 if ((q<=p)&&(q<=n)) {
502 i=petit_c(a,n,p,q);
503 while (i!=0) {
504 if (i!=q) {
505 echange_l(a,i,q,n,p);
506 echange_c(b,i,q,n,n);
507 echange_l(d,i,q,n,n);
508 }
509
510 c1=a+(q-1)+(q-1)*p;
511 value_assign(pivot,*c1);
512 if (value_neg_p(pivot)) {
513 moins_l(a,q,n,p);
514 moins_l(d,q,n,n);
515 moins_c(b,q,n,n);
516 value_oppose(pivot,pivot);
517 }
518 for(k=q+1;k<=n;k++) {
519 c1=c1+p;
520 if (value_notzero_p(*c1)) {
521 value_division(x,*c1,pivot);
522 value_modulus(tmp,*c1,pivot);
523 if (value_neg_p(tmp))
524 value_decrement(x,x);
525 value_oppose(x_inv,x);
526 ligne(a,k,q,x_inv,n,p);
527 colonne(b,q,k,x,n,n);
528 ligne(d,k,q,x_inv,n,n);
529 }
530 }
531 i=petit_c(a,n,p,q);
532 }
533 c1=a+(q-1);
534 value_assign(pivot,*(c1+(q-1)*p));
535 if (value_neg_p(pivot)) {
536 moins_l(a,q,n,p);
537 moins_l(d,q,n,n);
538 moins_c(b,q,n,n);
539 value_oppose(pivot,pivot);
540 }
541 if (value_notzero_p(pivot)) {
542 for(k=1;k<q;k++) {
543 if (value_notzero_p(*c1)) {
544 value_division(x,*c1,pivot);
545 value_modulus(tmp,*c1,pivot);
546 if (value_neg_p(tmp))
547 value_decrement(x,x);
548 value_oppose(x_inv,x);
549 ligne(a,k,q,x_inv,n,p);
550 colonne(b,q,k,x,n,n);
551 ligne(d,k,q,x_inv,n,n);
552 }
553 c1=c1+p;
554 }
555 }
556 hermite(a,b,d,n,p,q+1);
557 }
558 value_clear(pivot); value_clear(tmp);
559 value_clear(x); value_clear(x_inv);
560 return;
561 } /* hermite */
562
563 /*----------------------------------------------------------------------
564 B est une g_base de dimension n dont les p premiers vecteurs
565 colonnes sont colineaires aux vecteurs colonnes de A
566 C est l'invers de B
567
568 B is an integral (g_base ?) of dimension n whose p first columns
569 vectors are proportionnal to the column vectors of A. C is the
570 inverse of B.
571 ----------------------------------------------------------------------*/
572
573 /** Convert PolmattoDarmat :
574 *** This function converts the matrix of a Polylib to a int * as necessary
575 *** for the functions in Darte's implementation.
576 ***
577 *** Input : A Polylib Matrix ( a pointer to it );
578 *** Output : An int * with the matrix copied into it
579 **/
ConvertPolMattoDarMat(Matrix * A)580 static Value *ConvertPolMattoDarMat(Matrix *A ) {
581
582 int i,j;
583 Value *result;
584
585 result = (Value *)malloc(sizeof(Value) * A->NbRows * A->NbColumns);
586 for(i=0;i<A->NbRows * A->NbColumns;i++)
587 value_init(result[i]);
588
589 for (i = 0; i < A->NbRows; i++)
590 for (j = 0 ; j < A->NbColumns; j++)
591 value_assign(result[i*A->NbColumns + j],A->p[i][j]);
592 return result;
593 } /* ConvertPolMattoDarMat */
594
595 /** Convert DarmattoPolmat
596 *** This function converts the matrix from Darte representation to a matrix
597 *** in PolyLib.
598 ***
599 *** Input : The matrix (a pointer to it), number of Rows, number of columns
600 *** Output : The matrix of the PolyLib
601 ***
602 **/
603
ConvertDarMattoPolMat(Value * A,int NbRows,int NbCols)604 static Matrix *ConvertDarMattoPolMat (Value *A, int NbRows, int NbCols) {
605
606 int i,j;
607 Matrix *result;
608
609 result = Matrix_Alloc(NbRows, NbCols);
610
611 for (i = 0; i < NbRows; i ++)
612 for (j = 0; j < NbCols; j ++)
613 value_assign(result->p[i][j],A[i*NbCols + j]);
614 return result;
615 } /* ConvertDarMattoPolMat */
616
617 /**
618 *** Smith : This function takes a Matrix A of dim n * l as its input
619 *** and returns the three matrices U, V and Product such that
620 *** A = U * Product * V, where U is an unimodular matrix of dimension
621 *** n * n and V is an unimodular matrix of dimension l * l.
622 *** Product is a diagonal matrix of dimension n * l.
623 ***
624 *** We use Alan Darte's implementation of Smith for computing
625 *** the Smith Normal Form
626 **/
Smith(Matrix * A,Matrix ** U,Matrix ** V,Matrix ** Product)627 void Smith(Matrix *A, Matrix **U, Matrix **V, Matrix **Product)
628 {
629 int i;
630 Matrix *u, *v;
631
632 u = Identity(A->NbRows);
633 v = Identity(A->NbColumns);
634
635 *U = Identity(A->NbRows);
636 *V = Identity(A->NbColumns);
637
638 *Product = Matrix_Copy(A);
639 smith((*Product)->p_Init, u->p_Init, v->p_Init, (*U)->p_Init, (*V)->p_Init,
640 A->NbRows, A->NbColumns, 1);
641
642 Matrix_Free(u);
643 Matrix_Free(v);
644 } /* Smith */
645
646 /** Hermite :
647 *** This function takes a Matrix as its input and finds its HNF
648 *** ( Left form )
649 ***
650 *** Input : A Matrix A (The Matrix A is not necessarily a Polylib matrix.
651 *** It is just a matrix as far as Hermite is concerned. It
652 *** does not even check if the matrix is a Polylib matrix or not)
653 *** Output : The Hnf matrix H and the Unimodular matrix U such that
654 *** A = H * U.
655 ***
656 *** We use Alan Darte's implementation of Hermite to compute the HNF.
657 *** Alan Darte's implementation computes the Upper Triangular HNF.
658 *** So We work on the fact that if A = H * U then
659 *** A (transpose) = U(transpose) * H (transpose)
660 *** There are a set of interface functions written in Interface.c
661 *** which convert a matrix from Polylib representationt to that of
662 *** Alan Darte's and vice versa.
663 ***
664 *** This Function Does the Following
665 *** Step 1 : Given the matrix A it finds its Transpose.
666 *** Step 2 : Finds the HNF (Right Form) using Alan Darte's Algorithm.
667 *** Step 3 : The H1 and U1 obtained in Step2 are both Transposed to get
668 *** the actual H and U such that A = HU.
669 **/
Hermite(Matrix * A,Matrix ** H,Matrix ** U)670 void Hermite (Matrix *A, Matrix **H, Matrix **U) {
671
672 int i;
673 Matrix *transpose, *tempH, *tempU ;
674 Value *darte_matA, *darte_identite, *darte_id_inv;
675
676 transpose = Transpose(A);
677 darte_matA = ConvertPolMattoDarMat(transpose);
678
679 darte_identite = (Value *)malloc(sizeof(Value) * A->NbColumns* A->NbColumns);
680 darte_id_inv = (Value *) malloc(sizeof(Value) * A->NbColumns*A->NbColumns);
681 for (i=0; i< A->NbColumns * A->NbColumns; i++)
682 value_init(darte_identite[i]);
683 for (i=0; i< A->NbColumns * A->NbColumns; i++)
684 value_init(darte_id_inv[i]);
685
686 identite(darte_identite, transpose->NbRows, transpose->NbRows);
687 identite(darte_id_inv, transpose->NbRows, transpose->NbRows);
688 hermite(darte_matA, darte_identite, darte_id_inv,A->NbColumns,A->NbRows, 1);
689 Matrix_Free (transpose);
690 transpose = ConvertDarMattoPolMat(darte_matA, A->NbColumns, A->NbRows);
691 tempU = ConvertDarMattoPolMat(darte_identite, A->NbColumns, A->NbColumns);
692
693 /* Finding H Transpose */
694 tempH = Transpose(transpose);
695 Matrix_Free(transpose);
696
697 /* Finding U Transpose */
698 transpose = Transpose(tempU);
699
700 H[0] = Matrix_Copy(tempH);
701 U[0] = Matrix_Copy(transpose);
702
703 Matrix_Free (tempH);
704 Matrix_Free (transpose);
705 Matrix_Free (tempU);
706
707 for (i=0; i< A->NbRows * A->NbColumns; i++)
708 value_clear(darte_matA[i]);
709 for (i=0; i< A->NbColumns * A->NbColumns; i++)
710 value_clear(darte_identite[i]);
711 for (i=0; i< A->NbColumns * A->NbColumns; i++)
712 value_clear(darte_id_inv[i]);
713 free (darte_matA);
714 free (darte_identite);
715 free (darte_id_inv);
716 return;
717 } /* Hermite */
718
719
720