source/kernel/NormalForms.c

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