/* ========================== C MeatAxe =============================
   znullsp.c - Calculate null-space

   (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: znullsp.c,v 1.2 1993/10/20 18:17:07 mringe Exp $
 *
 * $Log: znullsp.c,v $
 * Revision 1.2  1993/10/20  18:17:07  mringe
 * MeatAxe-2.0, Phase II.
 *
 * Revision 1.1  1993/10/19  09:30:59  mringe
 * Initial revision
 *
 *
 */



#include <stdlib.h>
#include "meataxe.h"



long znullsp(matrix, nor, piv, nsp)
PTR matrix, nsp;
long nor, *piv;

{
    PTR x, y, op;
    int i;
    long noc = znoc;
    long dim;

    /* Make identity matrix in nsp
       --------------------------- */
    zsetlen(zfl,nor);
    op = zalloc((long)1);
    x = nsp;
    for (i = 1; i <= nor; ++i)
    {
	piv[i] = 0;
	zmulrow(x,F_ZERO);
	zinsert(x,i,F_ONE);
	zadvance(&x,(long)1);
    }

    /* Gaussian elimination
       -------------------- */
    x = matrix;
    y = nsp;
    for (i = 1; i <= nor; ++i)
    {
	PTR xx = matrix, yy = nsp;
	long k, p;
	FEL f;

	for (k = 1; k < i; ++k)
	{
	    if ((p = piv[k]) != 0 && (f = zextract(x,p)) != F_ZERO)
	    {
		f = zsub(F_ZERO,f);
		zsetlen(zfl,noc);
		zaddmulrow(x,xx,f);
		zsetlen(zfl,nor);
		zaddmulrow(y,yy,f);
	    }
	    zsetlen(zfl,noc);
	    zadvance(&xx,(long)1);;
	    zsetlen(zfl,nor);
	    zadvance(&yy,(long)1);
	}

	piv[i] = p = zfindpiv(x,&f);
	if (p != 0 && f != F_ONE)
	{
	    f = zdiv(F_ONE,f);
	    zsetlen(zfl,noc);
	    zmulrow(x,f);
	    zsetlen(zfl,nor);
	    zmulrow(y,f);
	}
	zsetlen(zfl,noc);
	zadvance(&x,(long)1);
	zsetlen(zfl,nor);
	zadvance(&y,(long)1);
    }

    /* Step 2: Make the null space in echelon form
       ------------------------------------------- */
    zsetlen(zfl,nor);
    dim = 0;
    y = nsp;
    for (i = 1; i <= nor; ++i)
    {
	if (piv[i] == 0)
	{
	    FEL f;
	    PTR yy = y;
	    long k, p;

	    ++dim;
	    p = zfindpiv(y,&f);
	    if (f != F_ONE) zmulrow(y,zdiv(F_ONE,f));
	    for (k = i; k <= nor; ++k)
	    {
		if (k > i && piv[k] == 0)
		{
		    FEL f1 = zextract(yy,p);
		    if (f1 != F_ZERO)
			zaddmulrow(yy,y,zsub(F_ZERO,f1));
		}
		zadvance(&yy,(long)1);
	    }
	}
	zadvance(&y,(long)1);
    }
    zsetlen(zfl,nor);
    piv[0] = nor - dim;
    return dim;
}