/*
** Ihermite.c Birk Huber, 2/99
** -- compute hermite normal form of an integer matrix
**
**
** TiGERS, Toric Groebner Basis Enumeration by Reverse Search
** copyright (c) 1999 Birk Huber
**
*/
#include<stdio.h>
#include<stdlib.h>
#include "utils.h"
/*
**
** Imatrix_hermite:
** Input: S an n by m Imatrix.
** Output: True
**
** Sidefects: U points to a unitary nxn matrix
** S gets replaced by its hermite normal form (U*S)
**
** Method:
** U = Inxn
** c = 1
** while (c<=n) do
** find mc>=c such that abs(S(mc,c)) has minimum non-zero value
** if (mc!=c) then switch rows mc and c (in S, and U)
** if (S(c,c)<0) then multiply row c by -1
** if (S(c,c)!=0) then
** for i in c+1:n add S(i,c)/S(c,c) times row c to row i;
** end(if)
** if S(i,c)==0 for all i in c+1 ... n set c=c+1
** end(while)
**
*/
#define Mref(U,i,j) U[(i)-1][(j)-1]
#define Mrows(U) m
#define Mcols(U) n
int ihermite(int **S,int **U,int m, int n){
int c=1,mc,i,j,done,sign;
int t=0,mv=0, mn,crk=0;
if (m>n) mn=n; else mn=m;
/* Initialize U to nxn identity */
for(i=1;i<=n;i++){
for(j=1;j<=n;j++){
if (i==j) Mref(U,i,j)=1;
else Mref(U,i,j)=0;
}
}
while(c<=mn){
/* find minimum entry in col c */
mv=Mref(S,c,c);
if (mv<0) mv*=-1;
mc=c;
for(i=c+1;i<=n;i++){
t=Mref(S,i,c);
if (t<0) t*=-1;
if(mv==0 || (mv>t && t!=0)){
mv=t;
mc=i;
}
}
/* if nescesary pivot to put min in row c and multiply by+-1
to ensure diagonal entry is positive */
if (mc!=c||Mref(S,mc,c)<0){
if (Mref(S,mc,c)<0) sign=-1;
else sign=+1;
for(j=c;j<=m;j++){
t=Mref(S,mc,j);
Mref(S,mc,j)=Mref(S,c,j);
Mref(S,c,j)=sign*t;
}
for(j=1;j<=n;j++){
t=Mref(U,mc,j);
Mref(U,mc,j)=Mref(U,c,j);
Mref(U,c,j)=sign*t;
}
}
/* if collumb is not zero do a reduction step */
done=TRUE;
if (Mref(S,c,c)!=0){
crk=c;
for(i=c+1;i<=n;i++){
t=Mref(S,i,c)/Mref(S,c,c);
for(j=c;j<=m;j++){
Mref(S,i,j)-=t*Mref(S,c,j);
}
for(j=1;j<=n;j++){
Mref(U,i,j)-=t*Mref(U,c,j);
}
if (Mref(S,i,c)!=0) done=FALSE;
}
}
/* if all entrees of col c bellow row c are zero go tonext col */
if (done==TRUE) c++;
}
return crk;
}