/* ========================== C MeatAxe =============================
   ymatrix.c -  Matrix operations

   (C) Copyright 1993 Michael Ringe, Lehrstuhl D fuer Mathematik,
   RWTH Aachen, Germany  <mringe@tiffy.math.rwth-aachen.de>
   This program is free software; see the file COPYING for details.
   ================================================================== */

/* $Id: ymatrix.c,v 2.4 1994/02/13 18:26:56 mringe Exp $
 *
 * $Log: ymatrix.c,v $
 * Revision 2.4  1994/02/13  18:26:56  mringe
 * Neu: os.c, os.h.
 *
 * Revision 2.3  1994/02/12  04:10:13  mringe
 * UMFANGREICHE AENDERUNGEN AN VIELEN DATENTYPEN.
 *
 * Revision 2.2  1993/12/07  14:34:24  mringe
 * Teste bei realloc() auf size==0.
 *
 * Revision 2.1  1993/10/20  18:17:07  mringe
 * MeatAxe-2.0, Phase II.
 *
 * Revision 2.0  1993/10/14  18:54:18  mringe
 * MeatAxe-2.0, Phase I
 *
 * Revision 1.17  1993/10/11  19:05:28  mringe
 * Neue Library-Struktur.
 *
 * Revision 1.16  1993/10/06  04:41:05  mringe
 * utils Library eliminiert.
 *
 * Revision 1.15  1993/10/05  18:57:28  mringe
 * Neue Fehlermeldungen.
 *
 * Revision 1.14  1993/10/05  11:54:02  mringe
 * Benutze zreadvec()/zwritevec()
 *
 * Revision 1.13  1993/10/02  16:23:02  mringe
 * matread() und matwrite() in matload() bzw. matsave() umbenannt.
 *
 * Revision 1.12  1993/10/02  16:09:23  mringe
 * prmat() in matprint umbenannt.
 *
 * Revision 1.11  1993/08/31  12:57:03  mringe
 * ERR_DIV0 ersetzt ERR_SINGULAR
 *
 * Revision 1.10  1993/08/12  16:14:23  mringe
 * Include <string.h>.
 *
 * Revision 1.9  1993/08/06  14:01:59  mringe
 * Neuer File-header.
 *
 */


#include <string.h>
#include <stdlib.h>

#include "meataxe.h"





#define FOREACHROW(m,i,x)\
    for ((i)=1,(x)=(m)->d;(i)<=(m)->nor;++(i),zadvance(&(x),(long)1))



/* ------------------------------------------------------------------
   Global data
   ------------------------------------------------------------------ */

static long *piv = NULL;		/* Pivot table */
static long maxnoc = -1;		/* Table size */


/* ------------------------------------------------------------------
   Function prototypes
   ------------------------------------------------------------------ */

static void pivalloc _PL((long noc));
static void matpwr_ _PL((long n, PTR inp, PTR out, PTR tmp2));
static long ggt _PL((long a, long b));
static long kgv _PL((long a, long b));



/* ------------------------------------------------------------------
   matdesc() - Describe a matrix
   ------------------------------------------------------------------ */
    
char *matdesc(m)
matrix_t *m;

{
    static char buf[50];
    sprintf(buf,"%ld x %ld matrix over GF(%ld)",
	m->nor, m->noc, m->fl);
    return buf;
}

/* ------------------------------------------------------------------
   matprint() - Print a matrix to stdout
   ------------------------------------------------------------------ */

void matprint(name,m)
char *name;
matrix_t *m;

{
    PTR x;
    long i, k;

    zsetlen(m->fl,m->noc);
    x = m->d;
    if (name != NULL) printf("%s=\n",name);
    for (i = 1; i <= m->nor; ++i)
    {
	for (k = 1; k <= m->noc; ++k)
	    printf("%d",zextract(x,k));
	printf("\n");
	zadvance(&x,(long)1);
    }
}



/* ------------------------------------------------------------------
   pivalloc() - Allocate pivot table. All functions in this module
	use one global pivot table (piv). pivalloc() checks first
	if the currently allocated table is large enough.
   ------------------------------------------------------------------ */

static void pivalloc(noc)
long noc;

{	if (noc <= maxnoc)
		return;
	maxnoc = noc;
	if (piv != NULL)
		free(piv);
	piv = (long *)malloc((size_t)(noc+1)*sizeof(long));
	if (piv == NULL)
	{
	    fprintf(stderr,"pivalloc(): out of memory\n");
	    exit(1);
	}
}


/* ------------------------------------------------------------------
   matalloc() - Allocate a matrix
   matfree() - Free a matrix
   ------------------------------------------------------------------ */

matrix_t *matalloc(fl, nor, noc)
long fl, nor, noc;

{
    matrix_t *m;

    if (fl == 0) fl = zfl;	/* Default value */
    zsetlen(fl,noc);

    m = (matrix_t *) malloc(sizeof(matrix_t));
    if (m == NULL)
	MTXFAIL(ERR_NOMEM,NULL);
    m->id = T_MATRIX;
    m->d = zalloc(nor);
    if (m->d == NULL)
    {
	free(m);
	MTXFAIL(ERR_NOMEM,NULL);
    }
    m->fl = fl;
    m->nor = nor;
    m->noc = noc;
    return m;
}


void matfree(m)
matrix_t *m;

{
    if (m->d != NULL)
	free(m->d);
    free(m);
}




/* ------------------------------------------------------------------
   matid() - Identity matrix
   ------------------------------------------------------------------ */

matrix_t *matid(fl,nor)
long fl,nor;

{
    matrix_t *m;
    PTR x;
    long i;

    m = matalloc(fl,nor,nor);
    if (m == NULL) return NULL;
    FOREACHROW(m,i,x)
    {
	zmulrow(x,F_ZERO);
	zinsert(x,i,F_ONE);
    }
    return m;
}




/* ------------------------------------------------------------------
   matdup() - Duplicate a matrix
   matmove() - Move a matrix
   ------------------------------------------------------------------ */

matrix_t *matdup(src)
matrix_t *src;

{
    matrix_t *m;

    m = matalloc(src->fl,src->nor,src->noc);
    if (m == NULL) return NULL;
    memcpy(m->d,src->d,zsize(src->nor));
    return m;
}


int matmove(dest,src)
matrix_t *dest, *src;

{
    if (dest->fl != src->fl || dest->nor != src->nor ||
	dest->noc != src->noc)
    {
	MTXFAIL(ERR_INCOMPAT,-1);
    }
    zsetlen(src->fl,src->noc);
    memcpy(dest->d,src->d,zsize(src->nor));
    return 0;
}



/* ------------------------------------------------------------------
   matextract() - Extract rows out of a matrix.
   ------------------------------------------------------------------ */

matrix_t *matextract(src,first,n)
matrix_t *src;
long first;		/* First row */
long n;			/* Number of rows */

{
    matrix_t *m;
    PTR p;

    zsetlen(src->fl,src->noc);
    if (first < 1 || first+n-1 > src->nor || n < 0)
    {
	MTXFAIL(ERR_RANGE,NULL);
    }
    m = matalloc(src->fl,n,src->noc);
    if (m == NULL) return NULL;
    p = src->d;
    zadvance(&p,first-1);
    memcpy(m->d,p,zsize(n));
    return m;
}


/* ------------------------------------------------------------------
   matread(), matwrite() - Read and write matrices
   ------------------------------------------------------------------ */

matrix_t *matread(f)
FILE *f;

{
    matrix_t *m;
    long hdr[3];

    if (zreadlong(f,hdr,3) != 3) MTXFAIL(ERR_FILEREAD,NULL);
    if (hdr[0] < 2) MTXFAIL(ERR_BADTYPE,NULL);
    m = matalloc(hdr[0],hdr[1],hdr[2]);
    if (m == NULL) return NULL;
    if (zreadvec(f,m->d,(size_t)m->nor) != m->nor)
    {
	matfree(m);
	return NULL;
    }
    return m;
}

int matwrite(f, mat)
FILE *f;
matrix_t *mat;

{
    long hdr[3];

    hdr[0] = mat->fl;
    hdr[1] = mat->nor;
    hdr[2] = mat->noc;
    if (zwritelong(f,hdr,3) != 3) MTXFAIL(ERR_FILEWRITE,-1);
    zsetlen(mat->fl,mat->noc);
    if (zwritevec(f,mat->d,(size_t)mat->nor) != mat->nor)
	return -1;
    return 0;
}



/* ------------------------------------------------------------------
   matload(), matsave() - Read and write matrices
   ------------------------------------------------------------------ */

matrix_t *matload(fn)
char *fn;

{
    FILE *f;
    matrix_t *m;

    if ((f = os_fopen(fn,FM_READ)) == NULL)
    {
	perror(fn);
	errexit(ERR_FILEREAD,fn);
    }
    m = matread(f);
    fclose(f);
    if (m == NULL)
    {
	mtxerror(fn);
	errexit(ERR_FILEREAD,fn);
    }
    return m;
}


int matsave(mat,fn)
matrix_t *mat;
char *fn;

{
    FILE *f;
    int i;

    if ((f = os_fopen(fn,FM_CREATE)) == NULL)
    {
	perror(fn);
	errexit(ERR_FILEWRITE,fn);
    }
    i = matwrite(f,mat);
    fclose(f);
    if (i != 0)
    {
	mtxerror(fn);
	errexit(ERR_FILEWRITE,fn);
    }
    return i;
}



/* ------------------------------------------------------------------
   matadd(), matmul() - Basic matrix arithmetic
   ------------------------------------------------------------------ */

matrix_t *matadd(dest, src)		/* dest += src */
matrix_t *dest, *src;

{
    PTR s, d;
    long i;

    if (src->fl != dest->fl || src->noc != dest->noc ||
	src->nor != dest->nor)
	MTXFAIL(ERR_INCOMPAT,NULL);
    zsetlen(src->fl,src->noc);
    d = dest->d;
    s = src->d;
    for (i = src->noc; i != 0; --i)
    {
	zaddrow(d,s);
	zadvance(&d,(long)1);
	zadvance(&s,(long)1);
    }
    return dest;
}



matrix_t *matmul(dest, src)		/* dest *= src */
matrix_t *dest, *src;

{
    PTR x, tmp, result;
    long i;

    if (src->fl != dest->fl || src->nor != dest->noc)
    {
	MTXFAIL(ERR_INCOMPAT,NULL);
    }
    zsetlen(src->fl,src->noc);
    result = tmp = zalloc(dest->nor);
    if (result == NULL)
    {
	MTXFAIL(ERR_NOMEM,NULL);
    }
    x = dest->d;
    for (i = dest->nor; i != 0; --i)
    {
	zsetlen(zfl,src->noc);
	zmaprow(x,src->d,src->nor,tmp);
	zadvance(&tmp,(long)1);
	zsetlen(zfl,dest->noc);
	zadvance(&x,(long)1);
    }
    free(dest->d);
    dest->d = result;
    dest->noc = src->noc;
    return dest;
}


/* ------------------------------------------------------------------
   mattr() - Transpose a matrix
   ------------------------------------------------------------------ */

matrix_t *mattr(src)
matrix_t *src;

{
    PTR s, d;
    long i, j;
    matrix_t *dest;

    dest = matalloc(src->fl,src->noc,src->nor);
    if (dest == NULL) return NULL;
    d = dest->d;
    for (i = 1; i <= src->noc; ++i)
    {
	zsetlen(zfl,dest->noc);
	zmulrow(d,F_ZERO);
	zsetlen(zfl,src->noc);
	FOREACHROW(src,j,s)
	    zinsert(d,j,zextract(s,i));
	zsetlen(zfl,dest->noc);
	zadvance(&d,(long)1);
    }
    return dest;
}


/* ------------------------------------------------------------------
   matinv() - Matrix inversion
   ------------------------------------------------------------------ */

matrix_t *matinv(src)
matrix_t *src;

{
    PTR tmp = NULL;	/* Workspace */
    matrix_t *dest;

    if (src->nor != src->noc)
    {
	MTXFAIL(ERR_NOTSQUARE,NULL);
    }
    dest = matid(src->fl,src->nor);
    if (dest == NULL) return NULL;

    /* Copy matrix into workspace
       -------------------------- */
    tmp = zalloc(src->nor);
    memcpy(tmp,src->d,zsize(src->nor));

    /* Inversion
       --------- */
    if (zmatinv(tmp,dest->d) != 0) 
    {
	matfree(dest);
	return NULL;
    }
    return dest;
}


/* ------------------------------------------------------------------
   nullity() - Put a matrix into semi-echelon form (without removing
	zero rows) and return its nullity.
   ------------------------------------------------------------------ */

long nullity(mat)
matrix_t *mat;

{
    long j, k, l;
    long rank = 0;
    PTR xj, xk;
    FEL f, g;

    zsetlen(zfl,mat->noc);
    pivalloc(mat->nor);
    for (j=1, xj=mat->d; j <= mat->nor; ++j, zadvance(&xj,(long)1))
    {
	for (k=1, xk=mat->d; k < j; ++k, zadvance(&xk,(long)1))
	{
	    if ((l = piv[k]) == 0) continue;
	    g = zsub(F_ZERO,zdiv(zextract(xj,l),zextract(xk,l)));
	    zaddmulrow(xj,xk,g);
	}
	if ((piv[j] = zfindpiv(xj,&f)) != 0)
		++rank;
    }
    return (mat->nor - rank);
}



/* ------------------------------------------------------------------
   nullsp() - Put matrix into full echelon form (without removing
	zero rows) and return its null space.
   ------------------------------------------------------------------ */

matrix_t *nullsp(mat)
matrix_t *mat;

{
    long j, k, l;
    long rank = 0;
    PTR xj, xk, yj, yk;
    FEL f, g;
    matrix_t *nsp;
    size_t size;

    pivalloc(mat->nor);
    nsp = matalloc(mat->fl,mat->nor,mat->nor);
    for (j=1, xj=nsp->d; j <= nsp->nor; ++j, zadvance(&xj,(long)1))
    {
	zmulrow(xj,F_ZERO);
	zinsert(xj,j,F_ONE);
    }

    xj = mat->d;
    yj = nsp->d;
    for (j=1; j <= mat->nor; ++j)
    {
	zsetlen(zfl,mat->noc);
	if ((l = piv[j] = zfindpiv(xj,&f)) != 0)
	{
	    ++rank;
	    xk = mat->d;
	    yk = nsp->d;
	    for (k=1; k <= mat->nor; ++k)
	    {	if (k != j)
		{
		    g = zsub(F_ZERO,zdiv(zextract(xk,l),f));
		    zsetlen(zfl,mat->noc);
		    zaddmulrow(xk,xj,g);
		    zsetlen(zfl,nsp->noc);
		    zaddmulrow(yk,yj,g);
		}
		zsetlen(zfl,mat->noc);
		zadvance(&xk,(long)1);
		zsetlen(zfl,nsp->noc);
		zadvance(&yk,(long)1);
	    }
	}
	zsetlen(zfl,mat->noc);
	zadvance(&xj,(long)1);
	zsetlen(zfl,nsp->noc);
	zadvance(&yj,(long)1);
    }
    zsetlen(mat->fl,nsp->noc);
    yj = yk = nsp->d;
    for (j=1; j <= nsp->nor; ++j)
    {
	if (piv[j] == 0)
	{
	    if (yk != yj)
		zmoverow(yk,yj);
	    zadvance(&yk,(long)1);
	}
	zadvance(&yj,(long)1);
    }
    nsp->nor = (mat->nor - rank);

    size = zsize(nsp->nor);
    if (size == 0) size = 1;
    nsp->d = (PTR) realloc(nsp->d,size);
    return nsp;
}


/* ------------------------------------------------------------------
   echelon() - Put a matrix into echelon form
   ------------------------------------------------------------------ */

matrix_t *echelon(src)
matrix_t *src;

{
    long j, k, l;
    long rank = 0;
    PTR xj, xk;
    FEL f, g;
    matrix_t *mat;
    size_t size;

    mat = matdup(src);
    if (mat == NULL) return NULL;
    pivalloc(mat->nor);
    zsetlen(mat->fl,mat->noc);
    xj = mat->d;
    for (j=1; j <= mat->nor; ++j)
    {
	if ((l = piv[j] = zfindpiv(xj,&f)) != 0)
	{
	    ++rank;
	    xk = mat->d;
	    for (k=1; k <= mat->nor; ++k)
	    {
		if (k != j)
		{	g = zsub(F_ZERO,
			  zdiv(zextract(xk,l),f));
			zaddmulrow(xk,xj,g);
		}
		zadvance(&xk,(long)1);
	    }
	}
	zadvance(&xj,(long)1);
    }
    xj = xk = mat->d;
    for (j=1; j <= mat->nor; ++j)
    {
	if (piv[j] != 0)
	{
	    if (xk != xj)
		zmoverow(xk,xj);
	    zadvance(&xk,(long)1);
	}
	zadvance(&xj,(long)1);
    }
    mat->nor = rank;
    size = zsize(rank);
    if (size == 0) size = 1;
    mat->d = (PTR) realloc(mat->d,size);
    return mat;
}


/* ------------------------------------------------------------------
   chkechelon() - Check if a matrix is in full echelon form
   ------------------------------------------------------------------ */

int chkechelon(mat)
matrix_t *mat;

{
    PTR x, y;
    long i, k, piv;
    FEL f;

    zsetlen(mat->fl,mat->noc);
    x = mat->d;
    for (i = 1; i <= mat->nor; ++i)
    {	piv = zfindpiv(x,&f);
    	if (f != F_ONE) return 1;
    	y = mat->d;
    	for (k = 1; k <= mat->nor; ++k, zadvance(&y,(long)1))
    	{
	    if (k == i) continue;
    	    if (zextract(y,piv) != F_ZERO) return 1;
    	}
    	zadvance(&x,(long)1);
    }
    return 0;
}


/* ------------------------------------------------------------------
   matpower() - Calculate the <n>-th power of <mat> using the
	binary method.
   ------------------------------------------------------------------ */

static void matpwr_(n,inp,out,tmp2)
long n;
PTR inp, out, tmp2;

{
    PTR x, y;
    long i;
    int first = 1;

    while (n > 0)
    {	if (n % 2 == 1)
	{
	    if (first)
	    {
		memcpy(out,inp,zsize(znoc));
		first = 0;
	    }
	    else
	    {
		x = out;
		for (i = 1; i <= znoc; ++i)
		{
		    zmaprow(x,inp,znoc,tmp2);
		    zmoverow(x,tmp2);
		    zadvance(&x,(long)1);
		}
	    }
	}
	x = inp;
	y = tmp2;
	for (i = 1; i <= znoc; ++i)
	{
	    zmaprow(x,inp,znoc,y);
	    zadvance(&x,(long)1);
	    zadvance(&y,(long)1);
	}
	memcpy(inp,tmp2,zsize(znoc));
	n /= 2;
    }
}


matrix_t *matpower(mat,n)
matrix_t *mat;
long n;

{
    matrix_t *result;
    PTR tmp, tmp2;
    if (mat->nor != mat->noc)
    {
	MTXFAIL(ERR_NOTSQUARE,NULL);
    }
    zsetlen(mat->fl,mat->noc);
    tmp = zalloc(znoc);
    memcpy(tmp,mat->d,zsize(znoc));
    tmp2 = zalloc(znoc);
    result = matalloc(mat->fl,mat->nor,mat->noc);
    if (result == NULL)
    {
	free(tmp);
	free(tmp2);
	return NULL;
    }
    matpwr_(n,tmp,result->d,tmp2);
    free(tmp);
    free(tmp2);
    return result;
}


/* ------------------------------------------------------------------
   matorder() - Order of a matrix. Returns -1 on failure.
   ------------------------------------------------------------------ */

long matorder(mat)
matrix_t *mat;

{	PTR m1, v1, v2, v3;
	PTR basis, bend, bptr;
	char *done;
	long dim;
	long j1, i;
	FEL f1;
	long tord;
	long ord;
	int flag;

    if (mat->nor != mat->noc)
    {
	MTXFAIL(ERR_NOTSQUARE,-1);
    }
    zsetlen(mat->fl,mat->noc);
    m1 = zalloc(mat->nor);
    memcpy(m1,mat->d,zsize(mat->nor));
    bend = basis = zalloc(mat->nor+1);
    pivalloc(mat->nor+1);	/* BUG */
    done = (char *) malloc((size_t)((mat->nor+1)));
    memset(done,0,(size_t)((mat->nor+1)));
    v1 = zalloc((long)1);
    v2 = zalloc((long)1);
    v3 = zalloc((long)1);
    tord = ord = 1;
    dim = 0;
    j1 = 1;
    while (dim < mat->nor && tord <= 1000 && ord <= 1000000)
    {
	/* Get next start vector
	   --------------------- */
	for (j1 = 1; j1 <= mat->nor && done[j1]; ++j1);
	if (j1 > mat->nor) break;	/* Done! */
	zmulrow(v1,F_ZERO);
	zinsert(v1,j1,F_ONE);

	/* Calculate order on cyclic subspace
	   ---------------------------------- */
	tord = 0;
	flag = 1;
	zmoverow(v3,v1);
	do
	{   zmoverow(v2,v3);
	    if (flag)
	    {   zmoverow(bend,v3);
	        bptr = basis;
	        for (i = 1; i <= dim; ++i)
	        {   f1 = zextract(bend,piv[i]);
		    if (f1 != 0)
			zaddmulrow(bend,bptr,zsub(F_ZERO,
			zdiv(f1,zextract(bptr,piv[i]))));
		    zadvance(&bptr,(long)1);
	        }

		if ((piv[dim+1] = zfindpiv(bend,&f1)) != 0)
		{   done[piv[++dim]] = 1;
		    zadvance(&bend,(long)1);
		}
		else
		    flag = 0;
	    }
	    zmaprow(v2,m1,mat->nor,v3);
	    ++tord;
	}
	while (tord <= 1000 && zcmprow(v3,v1) != 0);

	/* Order = kgv(orders on cyclic subspaces)
	   --------------------------------------- */
	ord = kgv(ord,tord);
    }
    free(done);
    free(v1);
    free(v2);
    free(v3);
    free(m1);
    free(basis);
    if (tord > 1000 || ord > 1000000)
	return -1;
    return ord;
}



/* ------------------------------------------------------------------
   ggt() - Greatest common divisor
   kgv() - Least common multiple
   ------------------------------------------------------------------ */

static long ggt(a,b)
long a,b;

{	if (a == 0) return b;
	if (b == 0) return a;
	while (1)
	{	if ((a %= b) == 0) return b;
		if ((b %= a) == 0) return a;
	}
}

static long kgv(a,b)
long a,b;

{	return (a / ggt(a,b)) * b;
}




/* ------------------------------------------------------------------
   matquot() - Projection on quotient
   ------------------------------------------------------------------ */

matrix_t *matquot(subsp, mat)
matrix_t *subsp, *mat;

{
    long sdim = subsp->nor;
    long qdim = subsp->noc - sdim;
    PTR q;
    matrix_t *result;

    if (subsp->fl != mat->fl || subsp->noc != mat->noc)
    {
	MTXFAIL(ERR_INCOMPAT,NULL);
    }
    pivalloc(subsp->nor);
    zsetlen(subsp->fl,subsp->noc);
    zmkpivot(subsp->d,subsp->nor,piv);
    zquotinit(subsp->d,subsp->nor,piv);
    q = zquot(mat->d,mat->nor);
    result = (matrix_t *) malloc(sizeof(matrix_t));
    result->id = T_MATRIX;
    result->fl = mat->fl;
    result->nor = mat->nor;
    result->noc = qdim;
    result->d = q;
    return result;
}