source/kernel/SolveDio.c

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