2 This file is part of PolyLib.
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.
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.
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/>.
19 #include <polylib/polylib.h>
21 static void RearrangeMatforSolveDio(Matrix
*M
);
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
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
33 * A pointer to the pointer of a vector is a solution to T.
36 * Output : The above matrix U and the vector T.
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.
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
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'.
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.
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.
77 * 6) If a solution exists, find the values of X using
81 int SolveDiophantine(Matrix
*M
, Matrix
**U
, Vector
**X
) {
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 */
91 fp
= fopen("_debug", "a");
92 fprintf(fp
,"\nEntered SOLVEDIOPHANTINE\n");
96 value_init(sum
); value_init(tmp
);
98 /* Ensuring that the first rank row of A contribute to the rank*/
100 RearrangeMatforSolveDio(A
);
101 temp
= Matrix_Alloc(A
->NbRows
-1, A
->NbColumns
-1);
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
]);
108 /* Copying C into a temp, ignoring the Homogeneous part */
109 C
= (Value
*) malloc (sizeof(Value
) * (A
->NbRows
-1));
112 for (i
= 0; i
< k1
; i
++) {
114 value_oppose(C
[i
],A
->p
[i
][A
->NbColumns
-1]);
118 /* Finding the HNF of temp */
119 Hermite(temp
, &hermi
, &unimod
);
121 /* Testing for existence of a Solution */
123 min
=(hermi
->NbRows
<= hermi
->NbColumns
) ? hermi
->NbRows
: hermi
->NbColumns
;
125 for (i
= 0; i
< min
; i
++) {
126 if (value_notzero_p(hermi
->p
[i
][i
]))
132 /* Solving the Equation using Gaussian Elimination*/
134 T
= (Value
*) malloc(sizeof(Value
) * temp
->NbColumns
);
135 k2
= temp
->NbColumns
;
139 for (i
= 0; i
< rank
; i
++) {
141 for (j
= 0; j
< i
; j
++) {
142 value_addmul(sum
, T
[j
], hermi
->p
[i
][j
]);
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 value_clear(sum
); value_clear(tmp
);
150 for (i
= 0; i
< k1
; i
++)
152 for (i
= 0; i
< k2
; i
++)
158 value_subtract(tmp
,C
[i
],sum
);
159 value_division(T
[i
],tmp
,hermi
->p
[i
][i
]);
162 /** Case when rank < Number of Columns; **/
164 for (i
= rank
; i
< hermi
->NbColumns
; i
++)
165 value_set_si(T
[i
],0);
167 /** Solved the equtions **/
168 /** When rank < hermi->NbRows; Verifying whether the solution agrees
169 with the remaining n-rank rows as well. **/
171 for (i
= rank
; i
< hermi
->NbRows
; i
++) {
173 for (j
= 0; j
< hermi
->NbColumns
; j
++) {
174 value_addmul(sum
, T
[j
], hermi
->p
[i
][j
]);
176 if (value_ne(sum
,C
[i
])) {
177 *U
= Matrix_Alloc(0,0);
178 *X
= Vector_Alloc (0);
179 value_clear(sum
); value_clear(tmp
);
180 for (i
= 0; i
< k1
; i
++)
182 for (i
= 0; i
< k2
; i
++)
189 unimodinv
= Matrix_Alloc(unimod
->NbRows
, unimod
->NbColumns
);
190 Matrix_Inverse(unimod
, unimodinv
);
192 *X
= Vector_Alloc(M
->NbColumns
-1);
194 if (rank
== hermi
->NbColumns
)
195 *U
= Matrix_Alloc(0,0);
196 else { /* Extracting the General solution form U(inverse) */
198 *U
= Matrix_Alloc(hermi
->NbColumns
, hermi
->NbColumns
-rank
);
199 for (i
= 0; i
< U
[0]->NbRows
; i
++)
200 for (j
= 0; j
< U
[0]->NbColumns
; j
++)
201 value_assign(U
[0]->p
[i
][j
],unimodinv
->p
[i
][j
+rank
]);
204 for (i
= 0; i
< unimodinv
->NbRows
; i
++) {
206 /* Calculating the vector X = Uinv * T */
208 for (j
= 0; j
< unimodinv
->NbColumns
; j
++) {
209 value_addmul(sum
, unimodinv
->p
[i
][j
], T
[j
]);
211 value_assign(X
[0]->p
[i
],sum
);
215 for (i = rank; i < A->NbColumns; i ++)
218 Matrix_Free (unimodinv
);
221 value_clear(sum
); value_clear(tmp
);
222 for (i
= 0; i
< k1
; i
++)
224 for (i
= 0; i
< k2
; i
++)
229 } /* SolveDiophantine */
233 * This function takes as input a matrix M (pointer to it)
234 * and it returns the tranformed matrix M, such that the first
235 * 'rank' rows of the new matrix M are the ones which contribute
236 * to the rank of the matrix M.
238 * 1) For a start we try to put all the zero rows at the end.
239 * 2) Then cur = 1st row of the remaining matrix.
240 * 3) nextrow = 2ndrow of M.
241 * 4) temp = cur + nextrow
242 * 5) If (rank(temp) == temp->NbRows.) {cur = temp;nextrow ++}
243 * 6) Else (Exchange the nextrow of M with the currentlastrow.
244 * and currentlastrow --).
245 * 7) Repeat steps 4,5,6 till it is no longer possible.
248 static void RearrangeMatforSolveDio(Matrix
*M
) {
250 int i
, j
, curend
, curRow
, min
, rank
=1;
252 Matrix
*A
, *L
, *H
, *U
;
254 /* Though I could have used the Lattice function
255 Extract Linear Part, I chose not to use it so that
256 this function can be independent of Lattice Operations */
258 L
= Matrix_Alloc(M
->NbRows
-1,M
->NbColumns
-1);
259 for (i
= 0; i
< L
->NbRows
; i
++)
260 for (j
= 0; j
< L
->NbColumns
; j
++)
261 value_assign(L
->p
[i
][j
],M
->p
[i
][j
]);
263 /* Putting the zero rows at the end */
264 curend
= L
->NbRows
-1;
265 for (i
= 0; i
< curend
; i
++) {
266 for (j
= 0; j
< L
->NbColumns
; j
++)
267 if (value_notzero_p(L
->p
[i
][j
]))
269 if (j
== L
->NbColumns
) {
270 ExchangeRows(M
,i
,curend
);
275 /* Trying to put the redundant rows at the end */
277 if (curend
> 0) { /* there are some useful rows */
280 A
= Matrix_Alloc(1, L
->NbColumns
);
282 for (i
= 0; i
<L
->NbColumns
; i
++)
283 value_assign(A
->p
[0][i
],L
->p
[0][i
]);
285 while (add
== True
) {
286 temp
= AddANullRow(A
);
287 for (i
= 0;i
<A
->NbColumns
; i
++)
288 value_assign(temp
->p
[curRow
][i
],L
->p
[curRow
][i
]);
289 Hermite(temp
, &H
, &U
);
290 for (i
= 0; i
< H
->NbRows
; i
++)
291 if (value_zero_p(H
->p
[i
][i
]))
293 if (i
!= H
->NbRows
) {
294 ExchangeRows(M
, curRow
, curend
);
301 A
= Matrix_Copy (temp
);
306 min
= (curend
>= L
->NbColumns
) ? L
->NbColumns
: curend
;
307 if (rank
==min
|| curRow
>= curend
)
314 } /* RearrangeMatforSolveDio */