PIPS
Functions
determinant.c File Reference
#include <stdlib.h>
#include <stdio.h>
#include "linear_assert.h"
#include "boolean.h"
#include "arithmetique.h"
#include "matrix.h"
+ Include dependency graph for determinant.c:

Go to the source code of this file.

Functions

void matrix_determinant (Pmatrix a, result)
 package matrix More...
 
void matrix_sub_determinant (Pmatrix a, int i, int j, result)
 void matrix_sub_determinant(Pmatrix a,int i, int j, int result[]) calculate sub determinant of a matrix More...
 

Function Documentation

◆ matrix_determinant()

void matrix_determinant ( Pmatrix  a,
result   
)

package matrix

int matrix_determinant(Pmatrix a, int * result):
calculate determinant of an (nxn) rational matrix a

The result consists of two integers, a numerator result[1] and a denominator result[0].

Algorithm: integer matrix triangularisation by Hermite normal form (because the routine is available)

Complexity: O(n**3) – up to gcd computations...

Authors: Francois Irigoin, September 1989, Yi-qing Yang, January 1990 output

cette routine est FAUSSE car la factorisation de Hermite peut s'appuyer sur des matrice unimodulaires de determinant 1 ou -1. Il faudrait donc ecrire un code de calcul de determinant de matrice unimodulaire...

Ou tout reprendre avec une triangularisation rapide et n'utilisant que des matrice de determinant +1. Voir mon code des US. Cette triangularisation donnera en derivation une routine d'inversion rapide pour matrice unimodulaire

matrix_hermite expects an integer matrix

product of diagonal elements

product of determinants of three matrixs P, H, Q

free useless matrices

reduce result

Definition at line 54 of file determinant.c.

57 {
58 
59  Value determinant_gcd;
60  Value determinant_p;
61  Value determinant_q;
62  int n = MATRIX_NB_LINES(a);
63  Value denominator = MATRIX_DENOMINATOR(a);
64  int i;
65  /* cette routine est FAUSSE car la factorisation de Hermite peut
66  s'appuyer sur des matrice unimodulaires de determinant 1 ou -1.
67  Il faudrait donc ecrire un code de calcul de determinant de
68  matrice unimodulaire...
69 
70  Ou tout reprendre avec une triangularisation rapide et n'utilisant
71  que des matrice de determinant +1. Voir mon code des US.
72  Cette triangularisation donnera en derivation une routine
73  d'inversion rapide pour matrice unimodulaire */
74 
75  assert(n > 0 && denominator != 0);
76 
77  switch(n) {
78 
79  case 1:
80  result[1] = MATRIX_ELEM(a,1,1);
81  result[0] = denominator;
82  break;
83 
84  case 2:
85  {
86  Value v1 = value_mult(MATRIX_ELEM(a,1,1),MATRIX_ELEM(a,2,2)),
87  v2 = value_mult(MATRIX_ELEM(a,1,2),MATRIX_ELEM(a,2,1));
88 
89  result[1] = value_minus(v1,v2);
90  result[0] = value_mult(denominator,denominator);
91  break;
92  }
93  default: {
94  Pmatrix p = matrix_new(n,n);
95  Pmatrix h = matrix_new(n,n);
96  Pmatrix q = matrix_new(n,n);
97 
98  /* matrix_hermite expects an integer matrix */
99  MATRIX_DENOMINATOR(a) = VALUE_ONE;
100  matrix_hermite(a,p,h,q,&determinant_p,&determinant_q);
101  MATRIX_DENOMINATOR(a) = denominator;
102 
103  /* product of diagonal elements */
104  result[1] = MATRIX_ELEM(h,1,1);
105  result[0] = denominator;
106  for(i=2; i <= n && result[1]!=0; i++) {
107  value_product(result[1],MATRIX_ELEM(h,i,i));
108  value_product(result[0],denominator);
109  determinant_gcd = pgcd_slow(result[0],result[1]);
110  if(value_notone_p(determinant_gcd)) {
111  value_division(result[0],determinant_gcd);
112  value_division(result[1],determinant_gcd);
113  }
114  }
115 
116  /* product of determinants of three matrixs P, H, Q */
117  value_product(result[1],determinant_p);
118  value_product(result[1],determinant_q);
119 
120  /* free useless matrices */
121  matrix_free(p);
122  matrix_free(h);
123  matrix_free(q);
124  }
125  }
126  /* reduce result */
127  determinant_gcd = pgcd_slow(result[0],result[1]);
128  if(value_notone_p(determinant_gcd)) {
129  value_division(result[0],determinant_gcd);
130  value_division(result[1],determinant_gcd);
131  }
132 }
value_minus
#define value_minus(v1, v2)
Definition: arithmetique-local.h:346
value_notone_p
#define value_notone_p(val)
Definition: arithmetique-local.h:393
Value
int Value
Definition: arithmetique-local.h:296
value_division
#define value_division(ref, val)
Definition: arithmetique-local.h:365
value_product
#define value_product(v, w)
Definition: arithmetique-local.h:447
VALUE_ONE
#define VALUE_ONE
Definition: arithmetique-local.h:303
value_mult
#define value_mult(v, w)
whether the default is protected or not this define makes no sense any more...
Definition: arithmetique-local.h:444
pgcd_slow
Value pgcd_slow(Value, Value)
pgcd.c
Definition: pgcd.c:44
MATRIX_NB_LINES
#define MATRIX_NB_LINES(matrix)
Definition: matrix-local.h:87
MATRIX_DENOMINATOR
#define MATRIX_DENOMINATOR(matrix)
int MATRIX_DENONIMATOR(matrix): acces au denominateur global d'une matrice matrix
Definition: matrix-local.h:86
matrix_free
#define matrix_free(m)
Allocation et desallocation d'une matrice.
Definition: matrix-local.h:73
MATRIX_ELEM
#define MATRIX_ELEM(matrix, i, j)
Macros d'acces aux elements d'une matrice.
Definition: matrix-local.h:80
matrix_new
Pmatrix matrix_new(int m, int n)
package matrix
Definition: alloc.c:41
matrix_hermite
void matrix_hermite(Pmatrix MAT, Pmatrix P, Pmatrix H, Pmatrix Q, Value *det_p, Value *det_q)
package matrix
Definition: hermite.c:78
assert
#define assert(ex)
Definition: newgen_assert.h:41
Pmatrix
package matrice
Definition: matrix-local.h:63

References assert, MATRIX_DENOMINATOR, MATRIX_ELEM, matrix_free, matrix_hermite(), MATRIX_NB_LINES, matrix_new(), pgcd_slow(), value_division, value_minus, value_mult, value_notone_p, VALUE_ONE, and value_product.

Referenced by matrix_sub_determinant().

+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ matrix_sub_determinant()

void matrix_sub_determinant ( Pmatrix  a,
int  i,
int  j,
result   
)

void matrix_sub_determinant(Pmatrix a,int i, int j, int result[]) calculate sub determinant of a matrix

The sub determinant indicated by (i,j) is the one obtained by deleting row i and column j from matrix of dimension n. The result passed back consists of two integers, a numerator result[1] and a denominator result[0].

Precondition: n >= i > 0 && n >= j > 0

Algorithm: put 0 in row i and in column j, and 1 in a[i,j], and then call matrix_determinant, restore a

Complexity: O(n**3) output

save column j

save row i, but for the (i,j) elements

restore row i

restore column j, with proper (i,j) element

Parameters
iThe sub determinant indicated by (i,j) is the one obtained by deleting row i and column j from matrix of dimension n. The result passed back consists of two integers, a numerator result[1] and a denominator result[0].

Precondition: n >= i > 0 && n >= j > 0

Algorithm: put 0 in row i and in column j, and 1 in a[i,j], and then call matrix_determinant, restore a

Complexity: O(n**3) input matrix

Definition at line 149 of file determinant.c.

153 {
154  int n = MATRIX_NB_LINES(a);
155  int r;
156  int c;
157  Pmatrix save_row = matrix_new(1,n);
158  Pmatrix save_column = matrix_new(1,n);
159 
160  /* save column j */
161  for(r=1; r <= n; r++) {
162  MATRIX_ELEM(save_column,1,r) = MATRIX_ELEM(a,r,j);
163  MATRIX_ELEM(a,r,j) = 0;
164  }
165  /* save row i, but for the (i,j) elements */
166  for(c=1; c <= n; c++) {
167  MATRIX_ELEM(save_row,1,c) = MATRIX_ELEM(a,i,c);
168  MATRIX_ELEM(a,i,c) = VALUE_ZERO;
169  }
170  MATRIX_ELEM(a,i,j) = VALUE_ONE;
171 
172  matrix_determinant(a,result);
173 
174  /* restore row i */
175  for(c=1; c <= n; c++) {
176  MATRIX_ELEM(a,i,c) = MATRIX_ELEM(save_row,1,c);
177  }
178  /* restore column j, with proper (i,j) element */
179  for(r=1; r <= n; r++) {
180  MATRIX_ELEM(a,r,j) = MATRIX_ELEM(save_column,1,r);
181  }
182 
183  matrix_free(save_row);
184  matrix_free(save_column);
185 }
VALUE_ZERO
#define VALUE_ZERO
Definition: arithmetique-local.h:302
matrix_determinant
void matrix_determinant(Pmatrix a, result)
package matrix
Definition: determinant.c:54

References matrix_determinant(), MATRIX_ELEM, matrix_free, MATRIX_NB_LINES, matrix_new(), VALUE_ONE, and VALUE_ZERO.

Referenced by matrix_triangular_inversion().

+ Here is the call graph for this function:
+ Here is the caller graph for this function: