source/kernel/SolveDio.c

   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 static void RearrangeMatforSolveDio(Matrix *M);
  22
  23 /*
  24  *  Solve Diophantine Equations :
  25  *        This function takes as input a system of equations in the form
  26  *          Ax + C = 0 and finds the solution for it, if it exists
  27  *
  28  *        Input : The matrix form the system of the equations Ax + C = 0
  29  *                 ( a pointer to a Matrix. )
  30  *                 A pointer to the pointer, where the matrix U
  31  *                  corresponding to the free variables of the equation
  32  *                  is stored.
  33  *                 A pointer to the pointer of a vector is a solution to T.
  34  *
  35  *
  36  *        Output : The above matrix U and the vector T.
  37  *
  38  *        Algorithm :
  39  *                    Given an integral matrix A, we can split it such that
  40  *                    A = HU, where H is in HNF (lowr triangular)
  41  *                      and U is unimodular.
  42  *                    So Ax = c -> HUx = c -> Ht = c ( where Ux = t).
  43  *                       Solving for Ht = c is easy.
  44  *                       Using 't' we find x = U(inverse) * t.
  45  *
  46  *        Steps :
  47  *                   1) For the above algorithm to work correctly to
  48  *                      need the condition that the first 'rank' rows are
  49  *                      the rows which contribute to the rank of the matrix.
  50  *                      So first we copy Input into a matrix 'A' and
  51  *                      rearrange the rows of A (if required) such that
  52  *                      the first rank rows contribute to the rank.
  53  *                   2) Extract A and C from the matrix 'A'. A = n * l matrix.
  54  *                   3) Find the Hermite normal form of the matrix A.
  55  *                       ( the matrices the lower tri. H and the unimod U).
  56  *                   4) Using H, find the values of T one by one.
  57  *                      Here we use a sort of Gaussian elimination to find
  58  *                      the solution. You have a lower triangular matrix
  59  *                      and a vector,
  60  *                      [ [a11, 0], [a21, a22, 0] ...,[arank1...a rankrank 0]]
  61  *                       and the solution vector [t1.. tn] and the vector
  62  *                      [ c1, c2 .. cl], now as we are traversing down the
  63  *                      rows one by one, we will have all the information
  64  *                      needed to calculate the next 't'.
  65  *
  66  *                      That is to say, when you want to calculate t2,
  67  *                      you would have already calculated the value of t1
  68  *                      and similarly if you are calculating t3, you will
  69  *                      need t1 and t2 which will be available by that time.
  70  *                      So, we apply a sort of Gaussian Elimination inorder
  71  *                      to find the vector T.
  72  *
  73  *                   5) After finding t_rank, the remaining (l-rank) t's are
  74  *                      made equal to zero, and we verify, if these values
  75  *                      agree with the remaining (n-rank) rows of A.
  76  *
  77  *                   6) If a solution exists, find the values of X using
  78  *                        U (inverse) * T.
  79  */
  80
  81 int SolveDiophantine(Matrix *M, Matrix **U, Vector **X) {
  82
  83   int i, j, k1, k2, min, rank;
  84   Matrix *A, *temp, *hermi, *unimod,  *unimodinv ;
  85   Value *C; /* temp storage for the vector C */
  86   Value *T; /* storage for the vector t */
  87   Value sum, tmp;
  88
  89 #ifdef DOMDEBUG
  90   FILE *fp;
  91   fp = fopen("_debug", "a");
  92   fprintf(fp,"\nEntered SOLVEDIOPHANTINE\n");
  93   fclose(fp);
  94 #endif
  95
  96   value_init(sum); value_init(tmp);
  97
  98   /* Ensuring that the first rank row of A contribute to the rank*/
  99   A = Matrix_Copy(M);
 100   RearrangeMatforSolveDio(A);
 101   temp = Matrix_Alloc(A->NbRows-1, A->NbColumns-1);
 102
 103   /* Copying A into temp, ignoring the Homogeneous part */
 104   for (i = 0; i < A->NbRows -1; i++)
 105     for (j = 0; j < A->NbColumns-1; j++)
 106       value_assign(temp->p[i][j],A->p[i][j]);
 107
 108   /* Copying C into a temp, ignoring the Homogeneous part */
 109   C = (Value *) malloc (sizeof(Value) * (A->NbRows-1));
 110   k1 = A->NbRows-1;
 111
 112   for (i = 0; i < k1; i++) {
 113     value_init(C[i]);
 114     value_oppose(C[i],A->p[i][A->NbColumns-1]);
 115   }
 116   Matrix_Free (A);
 117
 118   /* Finding the HNF of temp */
 119   Hermite(temp, &hermi, &unimod);
 120
 121   /* Testing for existence of a Solution */
 122
 123   min=(hermi->NbRows <= hermi->NbColumns ) ? hermi->NbRows : hermi->NbColumns ;
 124   rank = 0;
 125   for (i = 0; i < min ; i++) {
 126     if (value_notzero_p(hermi->p[i][i]))
 127       rank ++;
 128     else
 129       break ;
 130   }
 131
 132   /* Solving the Equation using Gaussian Elimination*/
 133
 134   T = (Value *) malloc(sizeof(Value) * temp->NbColumns);
 135   k2 = temp->NbColumns;
 136   for(i=0;i< k2; i++)
 137     value_init(T[i]);
 138
 139   for (i = 0; i < rank ; i++) {
 140     value_set_si(sum,0);
 141     for (j = 0; j < i; j++) {
 142       value_addmul(sum, T[j], hermi->p[i][j]);
 143     }
 144     value_subtract(tmp,C[i],sum);
 145     value_modulus(tmp,tmp,hermi->p[i][i]);
 146     if (value_notzero_p(tmp)) { /* no solution to the equation */
 147       *U = Matrix_Alloc(0,0);
 148       *X = Vector_Alloc (0);
 149       Matrix_Free (unimod);
 150       Matrix_Free (hermi);
 151       Matrix_Free (temp);
 152       value_clear(sum); value_clear(tmp);
 153       for (i = 0; i < k1; i++)
 154         value_clear(C[i]);
 155       for (i = 0; i < k2; i++)
 156         value_clear(T[i]);
 157       free(C);
 158       free(T);
 159       return (-1);
 160     };
 161     value_subtract(tmp,C[i],sum);
 162     value_division(T[i],tmp,hermi->p[i][i]);
 163   }
 164
 165   /** Case when rank < Number of Columns; **/
 166
 167   for (i = rank; i < hermi->NbColumns; i++)
 168     value_set_si(T[i],0);
 169
 170   /** Solved the equtions **/
 171   /** When rank < hermi->NbRows; Verifying whether the solution agrees
 172       with the remaining n-rank rows as well. **/
 173
 174   for (i = rank; i < hermi->NbRows; i++) {
 175     value_set_si(sum,0);
 176     for (j = 0; j < hermi->NbColumns; j++) {
 177       value_addmul(sum, T[j], hermi->p[i][j]);
 178     }
 179     if (value_ne(sum,C[i])) {
 180       *U = Matrix_Alloc(0,0);
 181       *X = Vector_Alloc (0);
 182       Matrix_Free (unimod);
 183       Matrix_Free (hermi);
 184       Matrix_Free (temp);
 185       value_clear(sum); value_clear(tmp);
 186       for (i = 0; i < k1; i++)
 187         value_clear(C[i]);
 188       for (i = 0; i < k2; i++)
 189         value_clear(T[i]);
 190       free(C);
 191       free(T);
 192       return (-1);
 193     }
 194   }
 195   unimodinv = Matrix_Alloc(unimod->NbRows, unimod->NbColumns);
 196   Matrix_Inverse(unimod, unimodinv);
 197   Matrix_Free(unimod);
 198   *X = Vector_Alloc(M->NbColumns-1);
 199
 200   if (rank == hermi->NbColumns)
 201     *U = Matrix_Alloc(0,0);
 202   else { /* Extracting the General solution form U(inverse) */
 203
 204     *U = Matrix_Alloc(hermi->NbColumns, hermi->NbColumns-rank);
 205     for (i = 0; i < U[0]->NbRows; i++)
 206       for (j = 0; j < U[0]->NbColumns; j++)
 207         value_assign(U[0]->p[i][j],unimodinv->p[i][j+rank]);
 208   }
 209
 210   for (i = 0; i < unimodinv->NbRows; i++) {
 211
 212     /* Calculating the vector X = Uinv * T */
 213     value_set_si(sum,0);
 214     for (j = 0; j < unimodinv->NbColumns; j++) {
 215       value_addmul(sum, unimodinv->p[i][j], T[j]);
 216     }
 217     value_assign(X[0]->p[i],sum);
 218   }
 219
 220   /*
 221     for (i = rank; i < A->NbColumns; i ++)
 222     X[0]->p[i] = 0;
 223   */
 224   Matrix_Free (unimodinv);
 225   Matrix_Free (hermi);
 226   Matrix_Free (temp);
 227   value_clear(sum); value_clear(tmp);
 228   for (i = 0; i < k1; i++)
 229     value_clear(C[i]);
 230   for (i = 0; i < k2; i++)
 231     value_clear(T[i]);
 232   free(C);
 233   free(T);
 234   return (rank);
 235 } /* SolveDiophantine */
 236
 237 /*
 238  * Rearrange :
 239  *            This function takes as input a matrix M (pointer to it)
 240  *            and it returns the tranformed matrix M, such that the first
 241  *            'rank' rows of the new matrix M are the ones which contribute
 242  *            to the rank of the matrix M.
 243  *
 244  *            1) For a start we try to put all the zero rows at the end.
 245  *            2) Then cur = 1st row of the remaining matrix.
 246  *            3) nextrow = 2ndrow of M.
 247  *            4) temp = cur + nextrow
 248  *            5) If (rank(temp) == temp->NbRows.) {cur = temp;nextrow ++}
 249  *            6) Else (Exchange the nextrow of M with the currentlastrow.
 250  *                     and currentlastrow --).
 251  *            7) Repeat steps 4,5,6 till it is no longer possible.
 252  *
 253  */
 254 static void RearrangeMatforSolveDio(Matrix *M) {
 255
 256   int i, j, curend, curRow, min, rank=1;
 257   Bool add = True;
 258   Matrix *A, *L, *H, *U;
 259
 260   /* Though I could have used the Lattice function
 261         Extract Linear Part, I chose not to use it so that
 262         this function can be independent of Lattice Operations */
 263
 264   L = Matrix_Alloc(M->NbRows-1,M->NbColumns-1);
 265   for (i = 0; i < L->NbRows; i++)
 266     for (j = 0; j < L->NbColumns; j++)
 267       value_assign(L->p[i][j],M->p[i][j]);
 268
 269   /* Putting the zero rows at the end */
 270   curend = L->NbRows-1;
 271   for (i = 0; i < curend; i++) {
 272     for (j = 0; j < L->NbColumns; j++)
 273       if (value_notzero_p(L->p[i][j]))
 274         break;
 275     if (j == L->NbColumns) {
 276       ExchangeRows(M,i,curend);
 277       curend --;
 278     }
 279   }
 280
 281   /* Trying to put the redundant rows at the end */
 282
 283   if (curend > 0) { /* there are some useful rows */
 284
 285     Matrix *temp;
 286     A = Matrix_Alloc(1, L->NbColumns);
 287
 288     for (i = 0; i <L->NbColumns; i++)
 289       value_assign(A->p[0][i],L->p[0][i]);
 290     curRow = 1;
 291     while (add == True ) {
 292       temp= AddANullRow(A);
 293       for (i = 0;i <A->NbColumns; i++)
 294         value_assign(temp->p[curRow][i],L->p[curRow][i]);
 295       Hermite(temp, &H, &U);
 296       for (i = 0; i < H->NbRows; i++)
 297         if (value_zero_p(H->p[i][i]))
 298           break;
 299       if (i != H->NbRows) {
 300         ExchangeRows(M, curRow, curend);
 301         curend --;
 302       }
 303       else {
 304         curRow ++;
 305         rank ++;
 306         Matrix_Free (A);
 307         A = Matrix_Copy (temp);
 308         Matrix_Free (temp);
 309       }
 310       Matrix_Free (H);
 311       Matrix_Free (U);
 312       min = (curend >= L->NbColumns) ? L->NbColumns : curend ;
 313       if (rank==min || curRow >= curend)
 314         break;
 315     }
 316     Matrix_Free (A);
 317   }
 318   Matrix_Free (L);
 319   return;
 320 } /* RearrangeMatforSolveDio */