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 /* computes c = lcm(a,b) using Gcd(a,b,&c) */
22 void Lcm3(Value a
, Value b
, Value
*c
)
26 if (value_zero_p(a
)) {
30 if (value_zero_p(b
)) {
35 value_multiply(tmp
, a
, b
);
36 value_absolute(tmp
, tmp
);
38 value_division(*c
, tmp
, *c
);
43 * Return the lcm of 'i' and 'j'
45 Value
*Lcm (Value i
, Value j
)
49 tmp
= (Value
*) malloc (sizeof(Value
));
56 * Return an identity matrix of size 'size'.
58 Matrix
*Identity(unsigned size
)
63 A
= Matrix_Alloc(size
, size
);
64 for (i
= 0; i
< size
; i
++)
65 value_set_si(A
->p
[i
][i
], 1);
70 * Exchange the rows 'Row1' and 'Row2' of the matrix 'M'
72 void ExchangeRows(Matrix
*M
, int Row1
, int Row2
) {
78 for (i
= 0; i
< (int)M
->NbColumns
; i
++) {
79 value_assign(temp
,M
->p
[Row1
][i
]);
80 value_assign(M
->p
[Row1
][i
],M
->p
[Row2
][i
]);
81 value_assign(M
->p
[Row2
][i
],temp
);
88 * Exchange the columns 'Column1' and 'Column2' of the matrix 'M'
90 void ExchangeColumns(Matrix
*M
, int Column1
, int Column2
) {
94 for (i
= 0; i
< M
->NbRows
; i
++)
95 value_swap(M
->p
[i
][Column1
],M
->p
[i
][Column2
]);
98 } /* ExchangeColumns */
101 * Return the Transpose of a matrix 'A'
103 Matrix
*Transpose (Matrix
*A
) {
108 transpose
= Matrix_Alloc (A
->NbColumns
, A
->NbRows
);
109 for (i
= 0; i
< (int)A
->NbRows
; i
++)
110 for (j
= 0; j
< (int)A
->NbColumns
; j
++)
111 value_assign(transpose
->p
[j
][i
],A
->p
[i
][j
]);
116 * Return a copy of the contents of a matrix 'Src'.
118 Matrix
*Matrix_Copy(Matrix
const *Src
)
123 Dst
= Matrix_Alloc(Src
->NbRows
, Src
->NbColumns
);
125 for (i
= 0; i
< Src
->NbRows
; i
++)
126 for (j
= 0; j
< Src
->NbColumns
; j
++)
127 value_assign(Dst
->p
[i
][j
],Src
->p
[i
][j
]);
132 * Test if the input matrix is integral or not.
134 Bool
isIntegral (Matrix
*A
) {
139 value_init(divisor
); value_init(tmp
);
140 value_assign(divisor
,A
->p
[A
->NbRows
-1][A
->NbColumns
-1]);
142 for (i
= 0; i
< A
->NbRows
; i
++)
143 for (j
= 0; j
< A
->NbColumns
; j
++) {
144 value_modulus(tmp
,A
->p
[i
][j
],divisor
);
145 if (value_notzero_p(tmp
)) {
146 value_clear(divisor
); value_clear(tmp
);
150 value_clear(divisor
); value_clear(tmp
);
155 * Check if the matrix 'A' is in Hermite normal form or not.
157 Bool
isinHnf(Matrix
*A
) {
164 temp
= Homogenise(A
,True
) ;
165 for (i
= 0; i
< temp
->NbRows
; i
++) {
166 value_assign(rem
,temp
->p
[i
][i
]);
167 for (j
= 0; j
< i
; j
++)
168 if (value_ge(temp
->p
[i
][j
],rem
)) {
173 for (j
= i
+1; j
< temp
->NbColumns
; j
++)
174 if (value_notzero_p(temp
->p
[i
][j
])) {
186 * Remove the row 'Rownumber' and place it at the end of the matrix 'X'
188 void PutRowLast (Matrix
*X
, int Rownumber
) {
193 if (Rownumber
== X
->NbRows
-1)
197 for (j
= 0; j
< X
->NbColumns
; j
++) {
198 value_assign(temp
,X
->p
[Rownumber
][j
]);
199 for (i
= Rownumber
; i
< X
->NbRows
-1; i
++)
200 value_assign(X
->p
[i
][j
],X
->p
[i
+1][j
]);
201 value_assign(X
->p
[i
][j
],temp
);
208 * Remove the row 'Rownumber' and place it at the begining of the matrix 'X'
210 void PutRowFirst(Matrix
*X
, int Rownumber
) {
216 for (j
= 0; j
< X
->NbColumns
; j
++) {
217 value_assign(temp
,X
->p
[Rownumber
][j
]);
218 for (i
= Rownumber
; i
> 0; i
--)
219 value_assign(X
->p
[i
][j
],X
->p
[i
-1][j
]);
220 value_assign(X
->p
[i
][j
],temp
);
227 * Remove the column 'Columnnumber' and place it at the begining of the matrix
230 void PutColumnFirst (Matrix
*X
, int Columnnumber
) {
236 for (i
= 0; i
< X
->NbRows
; i
++) {
237 value_assign(temp
,X
->p
[i
][Columnnumber
]);
238 for (j
= Columnnumber
; j
> 0; j
--)
239 value_assign(X
->p
[i
][j
],X
->p
[i
][j
-1]);
240 value_assign(X
->p
[i
][0],temp
);
244 } /* PutColumnFirst */
247 * Remove the column 'Columnnumber' and place it at the end of the matrix 'X'
249 void PutColumnLast (Matrix
*X
, int Columnnumber
) {
255 for (i
= 0; i
< X
->NbRows
; i
++) {
256 value_assign(temp
,X
->p
[i
][Columnnumber
]);
257 for (j
= Columnnumber
; j
< X
->NbColumns
-1; j
++)
258 value_assign(X
->p
[i
][j
],X
->p
[i
][j
+1]);
259 value_assign(X
->p
[i
][X
->NbColumns
-1],temp
);
263 } /* PutColumnLast */
266 * Add a row of zeros at the end of the matrix 'M' and return the new matrix
268 Matrix
*AddANullRow (Matrix
*M
) {
273 Result
= Matrix_Alloc(M
->NbRows
+1,M
->NbColumns
);
274 for (i
= 0;i
< M
->NbRows
; i
++)
275 for (j
= 0; j
< M
->NbColumns
; j
++)
276 value_assign(Result
->p
[i
][j
],M
->p
[i
][j
]);
277 for (j
= 0; j
< M
->NbColumns
; j
++)
278 value_set_si(Result
->p
[i
][j
],0);
283 * Add a column of zeros at the end of the matrix 'M' and return the new
286 Matrix
*AddANullColumn(Matrix
*M
) {
291 Result
= Matrix_Alloc(M
->NbRows
, M
->NbColumns
+1);
292 for (i
= 0;i
< M
->NbRows
; i
++)
293 for (j
= 0; j
< M
->NbColumns
; j
++)
294 value_assign(Result
->p
[i
][j
],M
->p
[i
][j
]);
295 for (i
= 0; i
< M
->NbRows
; i
++)
296 value_set_si(Result
->p
[i
][M
->NbColumns
],0);
298 } /* AddANullColumn */
301 * Remove a row 'Rownumber' from matrix 'M' and return the new matrix.
303 Matrix
*RemoveRow(Matrix
*M
, int Rownumber
) {
308 Result
= Matrix_Alloc(M
->NbRows
-1, M
->NbColumns
);
310 for (i
= 0; i
< Rownumber
; i
++)
311 Vector_Copy(M
->p
[i
], Result
->p
[i
], M
->NbColumns
);
312 for ( ; i
< Result
->NbRows
; i
++)
313 Vector_Copy(M
->p
[i
+1], Result
->p
[i
], M
->NbColumns
);
319 * Remove NumColumns columns starting at column number 'FirstColumnnumber' from matrix 'M'
320 * and return the new matrix.
322 Matrix
*RemoveNColumns (Matrix
*M
, int FirstColumnnumber
, int NumColumns
) {
327 Result
= Matrix_Alloc (M
->NbRows
, M
->NbColumns
-NumColumns
);
329 for (i
= 0; i
< Result
->NbRows
; i
++) {
330 Vector_Copy(M
->p
[i
], Result
->p
[i
], FirstColumnnumber
);
331 Vector_Copy(M
->p
[i
]+FirstColumnnumber
+NumColumns
, Result
->p
[i
]+FirstColumnnumber
,
332 M
->NbColumns
-NumColumns
-FirstColumnnumber
);
338 * Remove a column 'Columnnumber' from matrix 'M' and return the new matrix.
340 Matrix
*RemoveColumn (Matrix
*M
, int Columnnumber
) {
345 Result
= Matrix_Alloc (M
->NbRows
, M
->NbColumns
-1);
347 for (i
= 0; i
< Result
->NbRows
; i
++) {
348 Vector_Copy(M
->p
[i
], Result
->p
[i
], Columnnumber
);
349 Vector_Copy(M
->p
[i
]+Columnnumber
+1, Result
->p
[i
]+Columnnumber
,
350 M
->NbColumns
-1-Columnnumber
);
356 * Given a Matrix M of dimesnion n * l and rank l1, find a unimodular matrix
357 * 'Result' such that the Vector Space spanned by M is the subset of the vector
358 * Space spanned by the first l1 Rows of Result. The function returns the rank
359 * l1 and the Matrix Result.
361 int findHermiteBasis(Matrix
*M
, Matrix
**Result
) {
364 Matrix
*C
, *curMat
, *temp1
, *temp2
;
367 int dim
, curDim
, curVect
,rank
;
369 if (M
->NbRows
== 0) {
370 Result
[0] = Identity (M
->NbColumns
);
374 if (M
->NbRows
<= M
->NbColumns
) {
377 for (i
= 0; i
< H
->NbRows
; i
++)
378 if (value_zero_p(H
->p
[i
][i
]))
381 if (i
== H
->NbRows
) {
382 Result
[0] = Transpose(U
);
391 /* Eliminating the Zero Rows */
394 for (i
= 0; i
< C
->NbRows
; i
++) {
395 for (j
= 0; j
< C
->NbColumns
; j
++)
396 if (value_notzero_p(C
->p
[i
][j
]))
398 if (j
== C
->NbColumns
) {
400 temp
= RemoveRow(C
, i
);
402 C
= Matrix_Copy(temp
);
408 /* Eliminating the Redundant Rows */
414 curMat
= Matrix_Alloc(1,C
->NbColumns
);
415 for (i
= 0; i
< C
->NbColumns
; i
++)
416 value_assign(curMat
->p
[0][i
],C
->p
[0][i
]);
418 while((curVect
< C
->NbRows
) && (curDim
< dim
)) {
420 temp
= AddANullRow(curMat
);
421 for (i
= 0; i
< C
->NbColumns
; i
++)
422 value_assign(temp
->p
[temp
->NbRows
-1][i
],C
->p
[curVect
][i
]);
424 temp1
= AddANullRow(temp
);
425 temp2
= AddANullColumn(temp1
);
426 rank
= SolveDiophantine(temp2
, &U
, &V
);
427 if (rank
== temp
->NbRows
) {
429 curMat
= Matrix_Copy(temp
);
441 Hermite(curMat
, &H
, &U
);
442 rank
= curMat
->NbRows
;
445 Result
[0] = Transpose (U
);
449 } /* findHermiteBasis */