/* ========================== C MeatAxe =============================
   zcp.c -  Characteristic polynomial of a matrix

   (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: zcp.c,v 2.4 1994/03/14 12:59:28 mringe Exp $
 *
 * $Log: zcp.c,v $
 * Revision 2.4  1994/03/14  12:59:28  mringe
 * Option -f.
 *
 * Revision 2.3  1994/02/13  18:26:56  mringe
 * Neu: os.c, os.h.
 *
 * Revision 2.2  1993/10/26  10:47:35  mringe
 * Compiler Warnings.
 *
 * 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.15  1993/10/11  19:05:28  mringe
 * Neue Library-Struktur.
 *
 * Revision 1.14  1993/10/06  04:41:05  mringe
 * utils Library eliminiert.
 *
 * Revision 1.13  1993/08/10  14:51:42  mringe
 * Include string.h
 *
 * Revision 1.12  1993/08/06  14:01:59  mringe
 * Neuer File-header.
 *
 * Revision 1.11  1993/08/05  15:48:54  mringe
 * Neues message.c
 *
 * Revision 1.10  1993/07/23  11:40:20  mringe
 * GAP-Output: Benutze MeatAxe.CharPol.
 *
 * Revision 1.9  1993/07/17  09:14:30  mringe
 * Bugs korrigiert.
 *
 * Revision 1.8  1993/07/17  08:59:05  mringe
 * Benutze message library. Optionen: -G, -Q, -V.
 *
 * Revision 1.7  1993/07/13  20:30:59  mringe
 * Neue File i/o library.
 *
 * Revision 1.6  1993/05/21  07:55:49  mringe
 * Newlines auf dem Hilfstext entfernt.
 *
 * Revision 1.5  1993/02/17  11:16:12  mringe
 * Include-Files...
 *
 * Revision 1.4  1993/02/10  19:40:54  mringe
 * Libraries angelegt (YYY und ZZZ).
 *
 * Revision 1.3  1993/01/07  15:02:51  mringe
 * GAP-Output
 *
 */

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




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


static void init _PL((int argc, char *argv[]));
static void write_init _PL((void));
static void write_end _PL((void));
static void write_one _PL((poly_t *pol));
void spinup _PL((void));
long spinup_cyclic _PL((PTR base, long first));
void writeout _PL((long n));
int main _PL((int argc, char *argv[]));




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

long fl, nor, noc;		/* Field and size of the generators */
long *piv;			/* Pivot table */
char *ispiv;			/* Pivot flags */
PTR mat,			/* The matrix */
    A,				/* Work space (for spin-up) */
    B;				/* Work space II (coefficients) */
int opt_G = 0;			/* GAP output */
int opt_f = 0;			/* Factorization */
long dim = 0;			/* Dimension reached so far */
char *fname;			/* Matrix file name */
factor_t *factors = NULL;	/* Factorization */

static char *helptext[] = {
    "SYNTAX",
    "    zcp [-GQVf] <File>",
    "",
    "OPTIONS",
    "    -G   GAP output",
    "    -Q   Quiet, no messages",
    "    -V   Verbose, more messages",
    "    -f   Factorize the polynomial",
    "",
    "FILES",
    "    <File>    i  A square matrix",
    NULL
    };

static proginfo_t pinfo =
   { "zcp", "Characteristic Polynomial",
     "$Revision: 2.4 $", helptext };




/* ------------------------------------------------------------------
   init()
   ------------------------------------------------------------------ */

static void init(argc, argv)
int argc;
char *argv[];

{	
    int i;
    FILE *f;


    /* Parse command line
       ------------------ */
    mtxinit();
    initargs(argc, argv, &pinfo);
    while ((i = zgetopt("Gf")) != OPT_END)
    {
	switch (i)
	{
	    case 'G': opt_G = 1; msg_level = -100; break;
	    case 'f': opt_f = 1; break;
	}
    }
    if (opt_ind != argc - 1) errexit(ERR_NARGS,"zcp");
    fname = argv[opt_ind];

    /* Read the matrix
       --------------- */
    f = zreadhdr(fname,&fl,&nor,&noc);
    if (f == NULL) errexit(-1,fname);
    if (fl <= 1) errexit(ERR_NOTMATRIX,fname);
    if (nor != noc) errexit(ERR_NOTSQUARE,fname);
    zsetlen(fl,noc);
    mat = zalloc(nor);
    zreadvec(f,mat,(size_t)nor);
    fclose(f);
    MESSAGE(1,("%s: %ldx%ld-matrix over GF(%ld)\n",fname,nor,noc,fl));

    /* Allocate memory
       --------------- */
    A = zalloc(nor+1);
    B = zalloc(nor+1);
    piv = (long *) malloc((size_t)(nor+2)*sizeof(long));
    ispiv = (char *) malloc((size_t)(nor+2));
    memset(ispiv,0,(size_t)(nor+2));
}


/* ------------------------------------------------------------------
   write_init() - Called once before factors are written
   write_factor() - Write one factor
   write_end() - Called once after all factors have been written
   ------------------------------------------------------------------ */

static int first;

static void write_init()

{
    first = 1;
    if (opt_G)
	printf("MeatAxe.CharPol:=[\n");
    else
	printf("CHARACTERISTIC POLYNOMIAL:\n");
}


static void write_end()

{
    factor_t *f;

    if (opt_G)
	printf("];\n");
    else
    {
	if (opt_f)
	    for (f = factors; f->n != 0; ++f)
    	    {
		printf("( ");
		polprint(NULL,f->p);
		printf(" )^%ld\n",f->n);
    	    }
    }
}



static void write_one(pol)
poly_t *pol;

{ 
    if (opt_G)
    {
	int i;
	if (!first)
	    printf(",\n");
	printf("[");
	for (i = 0; i < pol->deg; ++i)
	    printf("%s,",zftogap(pol->buf[i]));
	printf("%s]",zftogap(F_ONE));
    }
    else if (!opt_f)
    {
	polprint(NULL,pol);
	printf("\n");
    }
    else
    {
	factor_t *f = factorization(pol);
    	if (factors == NULL)
	    factors = f;
    	else
    	{
	    int i, k, l, nf, nfact;

	    for (nf = 0; f[nf].n != 0; ++nf);
	    for (nfact = 0; factors[nfact].n != 0; ++nfact);
	    factors = realloc(factors,(size_t)(nf+nfact+1) *
		sizeof(factor_t));
	    for (i = 0; i < nf; ++i)
	    {
	    	for (k = 0; k < nfact; ++k)
		    if ((l = polcmp(factors[k].p,f[i].p)) >= 0) break;
	    	if (k == nfact || l > 0)
	    	{
		    int m;
		    for (m = nfact; m > k; --m)
		        factors[m] = factors[m-1];
		    factors[k] = f[i];
		    ++nfact;
	        }
	        else
	        {
		    factors[k].n += f[i].n;
		    polfree(f[i].p);
	        }
	    }
	    factors[nfact].p = NULL;
	    factors[nfact].n = 0;
            free(f);
    	}
    }
    first = 0;
    zsetlen(fl,noc);		/* Reset row size */
}

/* ------------------------------------------------------------------
   writeout() - Write out one factor of the polynomial
   ------------------------------------------------------------------ */

void writeout(n)
long n;

{
    int k;
    PTR x;
    static poly_t *pol = NULL;
    
    MESSAGE(2,("Cyclic subspace found, dim=%ld\n",n));
    if (pol == NULL) pol = polalloc(fl,nor);
    x = B;
    zadvance(&x,n);
    pol->deg = n;
    for (k = 1; k <= n; ++k)
	pol->buf[k-1] = zextract(x,k);
    pol->buf[n] = F_ONE;
    write_one(pol);
}


/* ------------------------------------------------------------------
   spinup()
   ------------------------------------------------------------------ */

void spinup()

{   PTR a;
    long i, n, dim;

    write_init();
    a = A;
    for (dim = 0; dim < nor; dim += n)
    {
    	for (i = 1; i <= nor && ispiv[i] != 0; ++i);
	MESSAGE(2,("Seed = %ld\n",i));
        zmulrow(a,F_ZERO);
	zinsert(a,i,F_ONE);
	n = spinup_cyclic(a,dim+1);
	zadvance(&a,n);
    }
    write_end();
}


/* ------------------------------------------------------------------
   spinup_cyclic() - Spin up one cyclic subaspace
   ------------------------------------------------------------------ */

long spinup_cyclic(base,first)
PTR base;
long first;

{   PTR a, b, x, y;
    long pv, i, k;
    FEL f;

    MESSAGE(2,("Next cyclic subspace, starting at %ld\n",first));
    a = base;
    b = B;
    zmulrow(b,F_ZERO);
    i = 0;
    while ((pv = zfindpiv(a,&f)) != 0)
    {	
    	/* Add new vector to basis
	   ----------------------- */
	if (MSG3 && i % 20 == 0)
	{
	    printf("dim=%ld\n",i);
	    fflush(stdout);
	}
	piv[first+i] = pv;
	ispiv[pv] = 1;
	++i;
	zinsert(b,i,F_ONE);

	/* Calculate the next vector
	   ------------------------- */
	x = a;
	zadvance(&a,(long)1);
	zmaprow(x,mat,nor,a);
	y = b;
	zadvance(&b,(long)1);
	zmulrow(b,F_ZERO);
	for (k = 1; k < nor; ++k)
	    zinsert(b,k+1,zextract(y,k));
    
	/* Clean with existing basis vectors
	   --------------------------------- */
	x = A;
	y = B;
	for (k = 1; k < first+i; ++k)
	{   f = zdiv(zextract(a,piv[k]),zextract(x,piv[k]));
	    zaddmulrow(a,x,zsub(F_ZERO,f));
	    if (k >= first)
	    {	zaddmulrow(b,y,zsub(F_ZERO,f));
		zadvance(&y,(long)1);
	    }
	    zadvance(&x,(long)1);
	}

    }
    writeout(i);
    return i;
}


/* ------------------------------------------------------------------
   main()
   ------------------------------------------------------------------ */

int main(argc, argv)
int argc;
char *argv[];

{	
	init(argc,argv);
	spinup();
	return EXIT_OK;
}