/* ========================== C MeatAxe =============================
   zmo.c - Make orbits under permutations

   (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: zmo.c,v 2.10 1994/03/12 13:03:53 mringe Exp $
 *
 * $Log: zmo.c,v $
 * Revision 2.10  1994/03/12  13:03:53  mringe
 * Verschiedene Bugs beseitigt.
 *
 * Revision 2.9  1994/03/09  14:03:42  mringe
 * Warnung, wenn bei -g > 1 permutation in einem File steht.
 *
 * Revision 2.8  1994/03/09  13:53:30  mringe
 * Option -g
 *
 * Revision 2.7  1994/03/09  12:14:33  mringe
 * Fehler im Sizes-File behoben.
 *
 * Revision 2.6  1994/02/15  10:28:33  mringe
 * MSDOS_BCC entfernt.
 *
 * Revision 2.5  1994/02/13  18:26:56  mringe
 * Neu: os.c, os.h.
 *
 * Revision 2.4  1993/10/26  10:47:35  mringe
 * Compiler Warnings.
 *
 * Revision 2.3  1993/10/21  21:57:35  mringe
 * Permutationen.
 *
 * Revision 2.2  1993/10/21  08:57:15  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.8  1993/10/11  19:05:28  mringe
 * Neue Library-Struktur.
 *
 * Revision 1.7  1993/10/06  04:41:05  mringe
 * utils Library eliminiert.
 *
 * Revision 1.6  1993/08/06  14:01:59  mringe
 * Neuer File-header.
 *
 * Revision 1.5  1993/08/06  12:53:45  mringe
 * help() eingebaut.
 *
 * Revision 1.4  1993/07/23  13:46:27  mringe
 * OS-Symbole neu (SYS_xxx)
 *
 * Revision 1.3  1993/07/13  20:30:59  mringe
 * Neue File i/o library.
 *
 * Revision 1.2  1992/06/27  17:20:11  mringe
 * Sagt jetzt, wenn es mehr als 200 verschiedene
 * Bahnl"angen gibt.
 *
 * Revision 1.1  1992/05/26  17:39:25  mringe
 * Initial revision
 *
 */


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




#define MAXPERM 10		/* Max. number of permutations */

#define MAXLEV 100000
#define MAXNORB 20000


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

static void err _PL((int c));



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

long *perm[MAXPERM+1];			/* Permutations */
long sstt[201], ors[MAXNORB+1];
long fl;
char *permname, *orbname, *sizname;

static char *helptext[] = {
"SYNTAX",
"    zmo [-QV] [-g <#Perms>] <Perms> <Orbits> <OrbSz>",
"",
"OPTIONS",
"    -Q   Quiet, no messages",
"    -V   Verbose output",
"    -g   Set number of permutations. Read from <Perm>.1, <Perm>.2 ...",
"",
"FILES",
"    <Perms>    i  A set of Permutations",
"    <Orbits>   o  Orbits under the permutations (Permutation format)",
"    <OrbSz>    o  Orbit sizes (Permutation format)",
"",
NULL};

static proginfo_t pinfo =
   { "zmo", "Make Orbits", "$Revision: 2.10 $", helptext };




/* ------------------------------------------------------------------
   err() - Print error message and exit
   ------------------------------------------------------------------ */

static void err(c)
int c;

{
    fprintf(stderr,"ZMO ERROR - ");
    switch (c)
    {
	case 'o':
	    fprintf(stderr,"MORE THAN %ld ORBITS\n",(long)MAXNORB);
	    break;
	}
	exit(EXIT_ERR);
}


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

{
    int ngen = -1;
    long i, nor, nop, seed, nout, siz, dis, norb;
    long *m1, *stk, *mout;
    long f1, f3,f5, lev;
    FILE *f;

    /* Parse command line
       ------------------ */
    mtxinit();
    initargs(argc, argv, &pinfo);
    while ((i = zgetopt("g:")) != OPT_END)
    {
	switch (i)
	{
	    case 'g':
		ngen = (int) getint();
		if (ngen < 1 || ngen > MAXPERM || *opt_text_ptr != 0)
		    errexit(ERR_OPTION,"-g");
		break;
	}
    }
    if (opt_ind != argc-3) errexit(ERR_NARGS,"zmo");
    permname = argv[opt_ind];
    orbname = argv[opt_ind+1];
    sizname = argv[opt_ind+2];

    /* Read permutations
       ----------------- */
    if (ngen == -1)
    {
    	if ((f = zreadhdr(permname,&fl,&nor,&nop)) == NULL)
	    errexit(-1,permname);
    	if (fl != -1) errexit(ERR_NOTPERM,argv[opt_ind]);
    	if (nop > MAXPERM)
    	{
	    fprintf(stderr,"ZMO WARNING - ONLY %ld OF %ld READ\n",
	    	(long) MAXPERM,nop);
	    nop = MAXPERM;
    	}
    	MESSAGE(1,("%s: %ld Permutation(s) of degree %ld\n",permname,
	    nop,nor));
        m1 = (long *) malloc(sizeof(long) * nor * nop);
        if (m1 == NULL) errexit(ERR_NOMEM,"zmo");
        if (zreadlong(f,m1,nor*nop) != nor * nop)
	    errexit(ERR_FILEREAD,permname);
        fclose(f);
        for (i = 1; i <= nop; ++i)
	    perm[i] = m1 + (nor * (i-1) - 1);
    }
    else
    {
	int k;
	nop = ngen;
	for (k = 1; k <= nop; ++k)
	{
	    char buf[200];
	    long nor2, nop2;
	    sprintf(buf,"%s.%d",permname,k);
    	    if ((f = zreadhdr(buf,&fl,&nor2,&nop2)) == NULL)
	        errexit(-1,buf);
    	    MESSAGE(1,("%s: 1 Permutation of degree %ld\n",buf,nor));
	    if (nop2 != 1)
	        fprintf(stderr,"ZMO WARNING - ONLY 1 OF %ld READ\n",
	    		nop2);
    	    if (fl != -1) errexit(ERR_NOTPERM,argv[opt_ind]);
	    if (k == 1)
		nor = nor2;
	    else if (nor != nor2)
		errexit(ERR_INCOMPAT,buf);
	    perm[k] = malloc(sizeof(long) * (size_t) nor);
	    if (perm[k] == NULL) errexit(ERR_NOMEM,"zmo");
	    if (zreadlong(f,perm[k],nor) != nor)
		errexit(ERR_FILEREAD,buf);
	    --perm[k];
	    fclose(f);
	}
    }

    stk = (long *) malloc(sizeof(long) * MAXLEV);
    mout = (long *) malloc(sizeof(long) * 250);
    if (stk == NULL || mout == NULL) errexit(ERR_NOMEM,"zmo");

    if ((f = zwritehdr(orbname,fl,nor,(long)1)) == NULL)
	errexit(-1,orbname);
    norb = 0;
    nout = 0;
    for (seed = 1; seed <= nor; ++seed)
    {
	if (norb >= MAXNORB) err('o');
	f1 = perm[1][seed];
	if (f1 < 0) continue;
	perm[1][seed] = -f1;	/* Mark as done */
	siz = 0;
	stk[lev = 1] = seed;	/* Put first point on stack */

	while (lev > 0)		/* Stack not empty */
	{
	    f1 = stk[lev--];	/* Take from stack */
	    if (nout == 250)	/* If buffer full, write out */
	    {	
		if (zwritelong(f,mout,250) != 250)
		    errexit(ERR_FILEWRITE,orbname);
		nout = 0;
	    }
	    MESSAGE(2,("mout[%ld]=%ld\n",(long)nout,(long)f1));
	    mout[nout++] = f1;	/* Put into buffer */
	    ++siz;

	    /* Apply all permutations and collect new points
	       --------------------------------------------- */
	    for (i = 1; i <= nop; ++i)
	    {	
		f3 = perm[i][f1];
	        MESSAGE(3,("perm%d(%ld)=%ld\n",(int)nop,f1,f3));
		if (f3 < 0) f3 = -f3;
		f5 = perm[1][f3];
		if (f5 <= 0) continue;	/* Already done, ignore */
		perm[1][f3] = -f5;	/* Mark as done */
		if (lev >= MAXLEV - 1)
		    errexit(ERR_NOMEM,"zmo (stack check)");
		stk[++lev] = f3;	/* Put on stack */
	    }
	}
	ors[norb++] = siz;
    }

    /* Flush buffer and close file
       --------------------------- */
    if (zwritelong(f,mout,nout) != nout)
   	 errexit(ERR_FILEWRITE,orbname);
    fclose(f);

    /* Write the orbit sizes file
       -------------------------- */
    if ((f = zwritehdr(sizname,fl,norb,(long)1)) == NULL)
	errexit(-1,sizname);
    if (zwritelong(f,ors,norb) != norb)
   	 errexit(ERR_FILEWRITE,sizname);
    fclose(f);

    /* Print orbit sizes
       ----------------- */
    if (msg_level >= 0)
    {
	dis = 0;
	for (i = 0; dis < 200 && i < norb; ++i)
	{
	    int k;
	    siz = ors[i];

	    for (k = 0; k < dis && ors[k] != siz; ++k);
	    if (k < dis) ++sstt[k];
	    else if (dis++ < 200)
	    {
		sstt[k] = 1;
		ors[k] = siz;
	    }
	}
	if (dis > 200)
	    printf("More than 200 distinct orbit sizes\n");
	else if (dis > 15)
	    printf("%ld distinct orbit sizes\n", (long) dis);
	else
	    for (i = 0; i < dis; ++i)
		printf("%6ld ORBIT%c OF SIZE %6ld\n",
		    sstt[i],sstt[i] > 1 ? 'S':' ',ors[i]);
    }
    return (EXIT_OK);
}