/* newgen.c  */

#include "def.h"
#include "macro.h"

static INT o_grenze;

static INT my_operate_perm_vector(ii,a,b,c) OP a,b,c; INT ii;
/****************************************************************
Input ii = INT, Index, der die Tiefe der Rekursion anzeigt, und damit
      zusaetzlich noch angibt, dass die anzuwendende Permutation die
      Stellen 0,...,ii-1 nicht veraendert.
      a = PERMUTATION
      b = VECTOR, auf dem a operiert, die Eintraege an den Stellen 
      0,...,ii-2 sind dabei nicht mehr angegeben
      natuerliche Operation einer Permutation aus S_n auf einem 
      n-elementigen Vector
Output c = VECTOR, Ergebnis der Operation von a auf b, a und b stimmen
       auf den Stellen 0,...,ii-1 ueberein, d.h. c wird erst ab der 
       Stelle ii angegeben.
Return OK
******************************************************************/
{
INT i;
INT erg = OK;
for (i=(INT)ii;i<S_V_LI(c);i++)
  M_I_I(S_V_II(b,i),S_V_I(c,S_P_II(a,i) -1L ));
return(OK);
}

static INT vergleich(a,b,j)
int j;
OP a,b;
{
int k;
for (k=j;k<S_V_LI(a);++k)
{
  if (S_V_II(a,k)<S_V_II(b,k)) return(4L);
  else if (S_V_II(a,k)>S_V_II(b,k))
  {
    if (k==j) return(1L);
    else return(2L);
  }
}
return(3L);
}

static INT my_sims_kette(a,in,b,anz,zibasis)
OP a,b,anz;
INT in,*zibasis;
/***********************************************************************
Input a = VECTOR von PERMUTATIONen, der ein Erzeugendensystem einer
      Permutationsgruppe enthaelt
      in = INT, Maechtigkeit des Definitionsbereichs
Output b = MATRIX, die ein System von strong generators der 
       Permutationsgruppe enthaelt, dabei sind alle Eintragungen 
       in b invertiert worden
       die Diagonalelemente sind stets die Identitaet
       anz = VECTOR, S_V_I(anz,i) enthaelt die Anzahl der Eintragungen
       in der i-ten Zeile von b, die verschieden von dem EMPTY Objekt
       sind. Diese Objekte liegen an den Stelle S_M_IJ(b,i,j) fuer
       i<=j<i+S_V_II(anz,i)
       zibasis = *INT, gibt die Laenge einer Basis der 
       Permutationsgruppe an, d.h. die hoechste Zeilennummer von b, 
       die nicht nur die Identitaet enthaelt.
Return OK oder ERROR;
***********************************************************************/
{
INT erg=OK;
INT i,j,l;
erg+=strong_generators(a,b);
erg+=m_il_v(in,anz);
for (i=0L;i<in;++i)
{
  l=0L;
  for (j=i;j<in ;++j)
  {
    if (not emptyp(S_M_IJ(b,i,j)))
    {
      if (j!=i+l) erg+=copy(S_M_IJ(b,i,j),S_M_IJ(b,i,i+l));
      l=l+1L;
    }
  }
  erg+=M_I_I(l,S_V_I(anz,i));
  for (j=i+l;j<S_M_LI(b);++j) if (!emptyp(S_M_IJ(b,i,j))) freeself(S_M_IJ(b,i,j));
}
i=in-1L;
while (S_V_II(anz,i)==1L && i>0L) i=i-1L;
*zibasis=i;
ENDR("my_sims_kette");
}

static INT next_vector_ogrenze(i,v,n)
OP v;
INT i,*n;
/***********************************************************************
Input i = INT, gibt an ob der zuletzt berechnete charakteristische 
      Vektor zu einem kanonischen Rep. gefuehrt hat (i=1L) oder 
      nicht (i=0L)
Input und
Output v = VECTOR, mit 0,1-Eintragungen, es wird lexikografisch der 
       naechstmoegliche Vektor (unter Verwendung der "Lerneffekte" 
       bestimmt, der eine maximal vorgegebene Anzahl von 1-en nicht 
       ueberschreitet)
       n = *INT (k+Anzahl der Einsen in w), wobei k die Dimension der 
       linearen Codes darstellt
Return 2L falls alle Kandidaten aufgelistet wurden
       1L sonst
***********************************************************************/
{
INT j,k=-1L,l,max=S_V_LI(v)-1L;
for (j=max;j>=0L;--j)
  if (S_V_II(v,j)==1L) 
  {
    k=j; j=0L;
  }
if (k==max)
{
  lab030697:
  l=k;
  for (j=k-1L;j>=0L;--j)
    if (S_V_II(v,j)==1L)
    {
      l=j; j=0L;
    }
  if (l==k) return(2L);
  /*******************************************
  alle Moeglichkeiten aufgelistet
  ********************************************/
  M_I_I(0L,S_V_I(v,l)); --(*n);
  ++l;
  if (S_V_II(v,l)==0L) 
  {
    M_I_I(1L,S_V_I(v,l));
    ++(*n);
  }
  for (j=++l;j<S_V_LI(v);++j)
    if (S_V_II(v,j)==1L) 
    {
      M_I_I(0L,S_V_I(v,j));
      --(*n);
    }
  return(1L);
}
++k;
if ((i==1L) && (*n<o_grenze))
{
  M_I_I(1L,S_V_I(v,k));
  ++(*n);
  return(1L);
}
else goto lab030697;
}

INT list_all_connected_graphs(g,og,vec)
OP g,og,vec;
{
OP b=callocobject();
OP anz=callocobject();
OP vv=callocobject();
OP w=callocobject(); /* neu */
INT in,ibasis,i,l,nn;
INT vert=2L; /* neu */
o_grenze=S_I_I(og);
in=S_P_LI(S_V_I(g,0L));
sims_kette(g,in,b,anz,&ibasis);
m_il_v(ibasis+2L,vv);
for (i=0L;i<S_V_LI(vv);++i) m_il_v(in,S_V_I(vv,i));
m_il_v(2L,vec);
m_il_nv(in,S_V_I(vec,0L));
m_il_nv(in,S_V_I(vec,1L));
inc(S_V_I(S_V_I(vec,1L),0L));
m_il_nv(in,S_V_I(vv,0L));
inc(S_V_I(S_V_I(vv,0L),0L)); /* neu */
m_il_v(in+1L,w); /* neu */
M_I_I(0L,S_V_I(w,0l)); /* neu */
M_I_I(0L,S_V_I(w,1l)); /* neu */
for (i=2L;i<S_V_LI(w);++i) M_I_I(S_V_II(w,i-1L)+i-1L,S_V_I(w,i)); /* neu */
nn=2L;
i=1L;
l=0L;
++l; /* neu */
/*while (next_vector_ogrenze(i,S_V_I(vv,0L),&l)==1L)*/
while (next_graph_ogrenze(i,S_V_I(vv,0L),&l,w,&vert)==1L)
{
/*printf("untersuchen ");println(S_V_I(vv,0L));*/
  i=untersuchen(0L,vv,ibasis,b,anz);
  if (i==1L)
  {
    /*printf("k=%d v=%d  ",l,vert);*/
    println(S_V_I(vv,0L));
    inc(vec);
    /*copy(S_V_I(vv,0L),S_V_I(vec,S_V_LI(vec)-1L));*/
    copy(S_V_I(vv,0L),S_V_I(vec,nn));
    ++nn;
    if (nn % 50L==0L) printf("%d\n",nn);
  }
}
freeall(b); freeall(anz); freeall(vv);
/* printf("%d\n",nn); */
}

INT list_all_reps(g,og,vec)
OP g,og,vec;
{
OP b=callocobject();
OP anz=callocobject();
OP vv=callocobject();
INT in,ibasis,i,l,nn;
o_grenze=S_I_I(og);
in=S_P_LI(S_V_I(g,0L));
sims_kette(g,in,b,anz,&ibasis);
m_il_v(ibasis+2L,vv);
for (i=0L;i<S_V_LI(vv);++i) m_il_v(in,S_V_I(vv,i));
m_il_v(1L,vec);
m_il_nv(in,S_V_I(vv,0L));
nn=1L;
i=1L;
l=0L;
while (next_vector_ogrenze(i,S_V_I(vv,0L),&l)==1L)
{
/*printf("untersuchen ");println(S_V_I(vv,0L));*/
  i=untersuchen(0L,vv,ibasis,b,anz);
  if (i==1L)
  {
    /*printf("k=%d v=%d  ",l,vert);*/
    println(S_V_I(vv,0L));
    /*inc(vec);
    copy(S_V_I(vv,0L),S_V_I(vec,S_V_LI(vec)-1L));*/
    ++nn;
    if (nn % 50L==0L) printf("%d\n",nn);
  }
}
freeall(b); freeall(anz); freeall(vv);
printf("%d\n",nn);
}

static INT next_kandidat(m,v)
OP v;
INT m;
/***********************************************************************
listet alle Vektoren der Laenge S_V_LI(v) auf mit Eintraegen in 
{0,...,m} auf.
Input m ist groesster Wert, der im VECTOR v auftreten darf
      v ein VECTOR mit INTEGER Eintraegen
Output v wird veraendert zu dem lexikografisch naechst groesseren
       VECTOR
Return 2L alle Vektoren aufgelistet
       1L noch nicht alle Vektoren aufgelistet
***********************************************************************/
{
int i,j,fertig;
i=S_V_LI(v)-1L;
fertig=0L;
while (fertig==0L && i >= 0L)
{ 
  if (S_V_II(v,i)<m)
  {
    inc(S_V_I(v,i));
    fertig=1L;
  }
  else
  {
    M_I_I(0L,S_V_I(v,i));
    --i;
  }
}
if (i<0L) return(2L); 
else return(1L);
}

static INT untersuchen(i,vv,basis,b,c)
int i,basis;
OP vv,b,c;
{
int j,k,l;
if (monoton_fallendp(S_V_I(vv,i),i)==1L) return(1L);
for (j=i; j<i+S_V_II(c,i); ++j)
{
  if (j==i)
  {
    for (k=i;k<S_V_LI(S_V_I(vv,i));++k) M_I_I(S_V_II(S_V_I(vv,i),k),S_V_I(S_V_I(vv,i+1L),k));
    k=vergleich(S_V_I(vv,0L),S_V_I(vv,i+1L),i);
    if (k==4L) return(0L);  
    else if ((k!=1L) && (i<basis))
    {
      k=untersuchen(i+1L,vv,basis,b,c);
      if (k==0L) return(0L);
    }
  }
  else
  {
    my_operate_perm_vector(i,S_M_IJ(b,i,j),S_V_I(vv,i),S_V_I(vv,i+1L));
    k=vergleich(S_V_I(vv,0L),S_V_I(vv,i+1L),i);
    if (k==4L) return(0L);  
    /* es gibt ein lexikogr. kleineres Element im Orbit  */
    else
    if ((k==2L) && (i<basis))
    {
      l=untersuchen(i+1L,vv,basis,b,c);
      if (l==0L) return(0L);
      /* es gibt ein lexikogr. kleineres Element im Orbit  */
    }
  }
}
return(1L);  
/* noch kein lexikogr. kleineres Element gefunden   */
}

static INT monoton_fallendp(a,j)
OP a;
INT j;
{
INT i;
i=j+1L;
if (S_V_II(a,j)==1L)
{
while ((i<S_V_LI(a)) && (S_V_II(a,i)==1L)) ++i;
}
for(;i<S_V_LI(a);++i) if (S_V_II(a,i)==1L) return(0L);
return(1L);
}

INT m_perm_2setsperm(a,b) OP a,b;
/* 140488 */
/* diese routine berechnet die induzierte permutation in n ueber 2
oder anders gesprochen:
berechnet die operation von pi aus S_n auf der identitaet 
in S_(n ueber 2) */
/* AK 170789 V1.0 */ /* AK 181289 V1.1 */ /* AK 150891 V1.3 */
{
OP c;
INT i,j,ni,nj,e=1L,z=1L;
INT erg = OK;
c = callocobject();
if (a == b) 
{ /* AK 291092 */
  *c = *b;
  C_O_K(b,EMPTY);
  erg += m_perm_2setsperm(c,b);
  erg += freeall(c);
  return erg;
}
erg += binom(S_P_L(a),cons_zwei,c);
/* c ist jetzt die laenge der ergebnis permutation */
/* c = n ueber 2 */
if (not EMPTYP(b)) /* AK 220891 */
  erg += freeself(b);
erg += b_ks_p(VECTOR,callocobject(),b);
erg += b_l_v(c,S_P_S(b));
/* die permutation ist nun initialisiert */
z=0L;
for(i=0L;i<S_P_LI(a);i++)
  for(j=0L;j<i;j++)
  {
    ni = S_P_II(a,i); nj = S_P_II(a,j);
    if (ni>nj) { e=ni; ni=nj; nj=e; };
    /* ni < nj ist ergebnis der permutation */
    /* nun nur noch die nummer zu bestimmen */
    /* sie ist e */
    e = ni+(nj-1L)*(nj-2L)/2L;
    /* e ist die nummer des neuen paars  */
    M_I_I(e,S_P_I(b,z)); 
    z++; 
  };
  /*println(b);*/
return erg;
}

static INT next_graph_ogrenze(i,v,n,w,vertices)
OP v,w;
INT i,*n,*vertices;
/***********************************************************************
Input i = INT, gibt an ob der zuletzt berechnete charakteristische 
      Vektor zu einem kanonischen Rep. gefuehrt hat (i=1L) oder 
      nicht (i=0L)
Input und
Output v = VECTOR, mit 0,1-Eintragungen, es wird lexikografisch der 
       naechstmoegliche Vektor (unter Verwendung der "Lerneffekte" 
       bestimmt, der eine maximal vorgegebene Anzahl von 1-en nicht 
       ueberschreitet)
       n = *INT (k+Anzahl der Einsen in w), wobei k die Dimension der 
       linearen Codes darstellt
Return 2L falls alle Kandidaten aufgelistet wurden
       1L sonst
***********************************************************************/
{
INT j,k=-1L,l,max=S_V_LI(v)-1L,m;
/*printf("i=%d n=%d vert=%d\n",i,*n,*vertices);
println(v);println(w);*/
for (j=max;j>=0L;--j)
  if (S_V_II(v,j)==1L) 
  {
    k=j; j=0L;
  }
if (k==max)
{
  lab030697:
  l=k;
  for (j=k-1L;j>=0L;--j)
    if (S_V_II(v,j)==1L)
    {
      l=j; j=0L;
    }
  if (l==k) return(2L);
  /*******************************************
  alle Moeglichkeiten aufgelistet
  ********************************************/
  for (j=*vertices;j>=1L;--j)
    if (l>=S_V_II(w,j-1L))
    {
      *vertices=j;
      if (l==S_V_II(w,*vertices-1L)) goto weiter;
      if (l+1L>=S_V_II(w,j))
      {
        for (m=S_V_II(w,*vertices-1L);m<l;++m)
	  if (S_V_II(v,m)==1L)  goto weiter;
      }
      else goto weiter;
      --(*vertices);
      k=S_V_II(w,*vertices);
      goto lab030697;
    }
  if (*vertices<=2L) return(2L);
  weiter:
  M_I_I(0L,S_V_I(v,l)); --(*n);
  ++l;
  if (S_V_II(v,l)==0L) 
  {
    M_I_I(1L,S_V_I(v,l));
    ++(*n);
  }
  if (l>=S_V_II(w,*vertices)) ++(*vertices);
  for (j=++l;j<S_V_LI(v);++j)
    if (S_V_II(v,j)==1L) 
    {
      M_I_I(0L,S_V_I(v,j));
      --(*n);
    }
  return(1L);
}
++k;
if ((i==1L) && (*n<o_grenze))
{
  if (k>=S_V_II(w,*vertices)) ++(*vertices);
  M_I_I(1L,S_V_I(v,k));
  ++(*n);
  return(1L);
}
else goto lab030697;
}

INT number_of_simple_graphs(a,d)
OP a,d;
/***********************************************************************
Input a = INTEGER
Output d = Anzahl der einfachen Graphen auf S_I_I(a) Knoten
Return OK oder ERROR;
***********************************************************************/
{
INT erg=OK;
OP b=callocobject();
OP c=callocobject();
if (S_O_K(a)!=INTEGER) WTO("number_of_simple_graphs",a);
if (!emptyp(d)) erg+=freeself(d);
erg+=zykelind_Sn(a,b);
erg+=zykelind_on_2sets(b,c);
erg+=polya_const_sub(c,cons_zwei,d);
erg+=freeall(b); erg+=freeall(c);
ENDR("number_of_simple_graphs");
}

INT numbers_of_simple_graphs(a,d)
OP a,d;
/***********************************************************************
Input a = INTEGER
Output d = VECTOR, S_V_I(d,i) ist die Anzahl der einfachen Graphen 
       auf i+1 Knoten
Return OK oder ERROR;
***********************************************************************/
{
OP hilf=callocobject();
INT i,erg=OK;
if (S_O_K(a)!=INTEGER) WTO("numbers_of_simple_graphs",a);
if (!emptyp(d)) erg+=freeself(d);
erg+=m_il_v(S_I_I(a),d);
for (i=0L,M_I_I(1L,hilf);i<S_I_I(a);++i,inc(hilf))
  erg+=number_of_simple_graphs(hilf,S_V_I(d,i));
erg+=freeall(hilf);
ENDR("numbers_of_simple_graphs");
}

INT numbers_of_connected_simple_graphs(a,v)
OP a,v;
/***********************************************************************
Input a = INTEGER
Output v = VECTOR, S_V_I(v,i) ist die Anzahl der einfachen 
       zusammenhaengenden Graphen auf i+1 Knoten
Return OK oder ERROR;
***********************************************************************/
{
OP help=callocobject();
OP part=callocobject();
OP sg=callocobject();
OP hilf=callocobject();
OP hilf1=callocobject();
OP hilf2=callocobject();
OP hilf3=callocobject();
INT i,j,erg=OK;
if (S_O_K(a)!=INTEGER) WTO("numbers_of_connected_simple_graphs",a);
if (!emptyp(v)) erg+=freeself(v);
erg+=numbers_of_simple_graphs(a,sg);
erg+=m_il_v(S_I_I(a),v);
erg+=M_I_I(1L,S_V_I(v,0L));
if (S_I_I(a)>1L) erg+=M_I_I(1L,S_V_I(v,1L));
i=2L;
for (M_I_I(3L,help);S_I_I(help)<=S_I_I(a);inc(help))
{
  erg+=m_i_i(0L,hilf3);
  first_part_EXPONENT(help,part);
  next(part,part);
  do
  {
  /*println(part);*/
  m_i_i(1L,hilf2);
  for (j=1L;j<S_PA_LI(part)-1L;++j)
    if (S_PA_II(part,j)>0L)
    {
      erg+=zykelind_Sn(S_PA_I(part,j),hilf);
      erg+=polya_const_sub(hilf,S_V_I(v,j),hilf1);
      erg+=mult_apply(hilf1,hilf2);
    }
  /*println(hilf2);*/
  erg+=add_apply(hilf2,hilf3);
  } while(next(part,part));
  /*println(hilf3);*/
  erg+=sub(S_V_I(sg,i),hilf3,S_V_I(v,i));
  ++i;
  /*println(v);*/
}
erg+=freeall(help); erg+=freeall(part); erg+=freeall(hilf); 
erg+=freeall(sg); erg+=freeall(hilf1); erg+=freeall(hilf2); 
erg+=freeall(hilf3);
ENDR("numbers_of_connected_simple_graphs");
}