/*   CF.c    */

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

INT inumber_of_fixed_points();
INT number_of_fixed_points();
INT inumber_of_cycles();
INT number_of_cycles();

INT CF();
INT CF_from_zykelind();
INT CF_right_sur();
INT CF_left_sur();
INT CF_left_right_sur();
INT CF_right_inj();
INT CF_left_inj();
INT CF_left_right_inj();

INT gen_arb();
INT gen_Sn();
INT gen_Cn();
INT gen_Dn();
INT gen_In();
INT gen_An();
INT gen_An_3();
INT gen_M24();
INT dir_sum_perm();
INT gen_dir_sum();
INT dir_prod_perm();
INT gen_dir_prod();
INT gen_on2sets();
INT generator_test();

INT t_kranz_plethysm();
OP s_kr_gi();
OP S_KR_GI();
INT s_kr_gii();
INT S_KR_GII();
OP s_kr_ij();
OP S_KR_IJ();
INT s_kr_iji();
INT S_KR_IJI();
INT gen_full_mon_group();
INT gen_plethysm();
INT gen_centralizer_permutation();
INT gen_stabilizer_partition();
INT gen_young_group();

INT zykelind_In();
INT degree_from_zykelind();
INT zykelind_M24();

INT Bell_numbers();
INT Bell_number();
INT Bell_number_recursive();
static INT bellrec();

INT polya_compl_sub();

static INT sum_ueber_alle_teiler();
INT comp_le();
INT prod_binom();

static INT first_unordered_part_into_atmost_k_parts();
static INT next_unordered_part_into_atmost_k_parts();

/* ****************************************************************

Some basic routines for enumeration under finite group actions.

******************************************************************* */

INT inumber_of_fixed_points(a)
OP a;
/* a eine PERMUTATION
   int Ergebnis ist Anzahl der Fixpunkte
*/
{
INT i,j;
j=0L;
for (i=0L;i<S_P_LI(a);++i)
{
  if (S_P_II(a,i)==i+1L) ++j;
}
return(j);
}

INT number_of_fixed_points(a,b)
OP a,b;
/* a eine PERMUTATION
   b = Anzahl der Fixpunkte
*/
{
INT i,j;
m_i_i(0L,b);
for (i=0L;i<S_P_LI(a);++i)
{
  if (S_P_II(a,i)==i+1L) inc(b);
}
}

INT number_of_cycles(a,b)
OP a,b;
/* a = PERMUTATION
   b = Anzahl der Zyklen von a 
*/
{
OP hilf=callocobject();
zykeltyp(a,hilf);
copy(s_pa_l(hilf),b);
freeall(hilf);
}

INT inumber_of_cycles(a)
OP a;
/* a = PERMUTATION
   gibt die Anzahl der Zyklen von a als integer Ergebnis weiter
*/
{
INT i;
OP hilf=callocobject();
zykeltyp(a,hilf);
i=s_pa_li(hilf);
freeall(hilf);
return(i);
}

/* *****************************************************************

Some methods for computing the Cauchy Frobenius Formula in different
situations. For instance from a set of generators of a permutation 
group, or from a cyle index. Furthermore there are routines for the
Cauchy Frobenius Formula for various group actions on the sets of 
functions f : X -> Y where the group actions on X an Y are given by
their cycle index.

******************************************************************** */

INT CF(a,b)
OP a,b;
/* a ist ein VECTOR, dessesn Komponenten PERMUTATIONen sind. Dies sind
   die Elemente der operierenden Gruppe.
   b ist Anzahl der G-Orbiten.
*/
{
INT i,j;
OP hilf=callocobject();
copy(cons_null,b);
for (i=0L;i<S_V_LI(a);++i)
{
  j=inumber_of_fixed_points(S_V_I(a,i));
  m_i_i(j,hilf);
  add_apply(hilf,b);
}
m_i_i(s_v_li(a),hilf);
ganzdiv(b,hilf,b);
freeall(hilf);
}

INT CF_from_zykelind(a,b)
OP a,b;
/* war auch Cauchy_Frobenius_fast  */
/* a ist ein Zyklenzeiger der Aktion von G auf X
   b ist Anzahl der G-orbiten auf X
*/
{
OP monom;
OP hilf=callocobject();
copy(cons_null,b);
monom=a;
while (monom!=NULL)
{
  mult(S_PO_K(monom),S_PO_SI(monom,0L),hilf);
  add_apply(hilf,b);
  monom=S_PO_N(monom);
}
freeall(hilf);
}

/*CF_polya_sit_zyk_sur(a,b,c)*/
INT CF_right_sur(a,b,c)
OP a,b,c;
{
INT i,j;
OP hilf=callocobject();
OP hilf1=callocobject();
OP hilf2=callocobject();
OP hilf3=callocobject();
OP monom;
monom=a;
copy(cons_null,hilf3);
for (i=0L;i<S_V_LI(S_PO_S(monom));++i) 
{
  m_i_i(i+1L,hilf); 
  mult_apply(S_V_I(S_PO_S(monom),i),hilf);
  add_apply(hilf,hilf3);
}
if (lt(hilf3,b)) return m_i_i(0L,c);
copy(cons_null,c);
while (monom!=NULL)
{
  for (j=1L;j<=S_I_I(hilf3);++j)
  {
    m_i_i(j,hilf);
    binom(b,hilf,hilf1);
    if (!emptyp(hilf2)) freeself(hilf2);
    sum_vector(S_PO_S(monom),hilf2);
    hoch(hilf,hilf2,hilf2);
    mult_apply(hilf1,hilf2);
    mult_apply(S_PO_K(monom),hilf2);
    sub(b,hilf,hilf1);
    if (evenp(hilf1)) add_apply(hilf2,c);
    else sub(c,hilf2,c);
  }
monom=S_PO_N(monom);
}
freeall(hilf); freeall(hilf1); freeall(hilf2); freeall(hilf3);
}

/*CF_action_from_left_zyk_sur(a,b,c)*/
INT CF_left_sur(a,b,c)
OP a,b,c;
{
INT i,j,k;
OP hilf=callocobject(); OP hilf1=callocobject(); OP hilf2=callocobject();
OP hilf3=callocobject(); OP hilf4=callocobject(); OP hilf5=callocobject();
OP hilf6=callocobject();
OP monom;
monom=a;
copy(cons_null,hilf6);
for (i=0L;i<S_V_LI(S_PO_S(monom));++i) 
{
  m_i_i(i+1L,hilf); 
  mult_apply(S_V_I(S_PO_S(monom),i),hilf);
  add_apply(hilf,hilf6);
}
if (GT(hilf6,b)) return m_i_i(0L,c);
copy(cons_null,c);
while (monom!=NULL)
{
  while (nullp(S_PO_SI(monom,S_V_LI(S_PO_S(monom))-1L))) dec(S_PO_S(monom));
  if (!emptyp(hilf)) freeself(hilf);
  sum_vector(S_PO_S(monom),hilf);
  for (j=1L;j<=S_I_I(hilf);++j)
  {
    M_I_I(j,hilf1);
    copy(cons_null,hilf5);
    first_unordered_part_into_atmost_k_parts(j,S_V_LI(S_PO_S(monom)),hilf2);
    do
    {
      copy(cons_eins,hilf3);
      for (k=0L;k<S_V_LI(hilf2);++k)
      {
        binom(S_V_I(S_PO_S(monom),k),S_V_I(hilf2,k),hilf4);
        mult_apply(hilf4,hilf3);
      }
      hoch(S_V_I(hilf2,0L),b,hilf4);
      mult_apply(hilf4,hilf3);
      add_apply(hilf3,hilf5);
      k=next_unordered_part_into_atmost_k_parts(hilf2);
    } while (k==1L);
    mult_apply(S_PO_K(monom),hilf5);
    sub(hilf,hilf1,hilf1);
    if (evenp(hilf1)) add_apply(hilf5,c);
    else sub(c,hilf5,c);
  }
monom=S_PO_N(monom);
}
freeall(hilf); freeall(hilf1); freeall(hilf2); 
freeall(hilf3); freeall(hilf4); freeall(hilf5); freeall(hilf6); 
}

/*CF_action_from_left_right_zyk_sur(a,b,c)*/
INT CF_left_right_sur(a,b,c)
OP a,b,c;
/* a = Gruppe, die auf dem Definitionsbereich operiert
   b = Gruppe, die auf dem Bildbereich operiert
   c = Ergebnis
*/
{
INT i,j,k,l;
OP hilf=callocobject(); OP hilf1=callocobject(); OP hilf2=callocobject();
OP hilf3=callocobject(); OP hilf4=callocobject(); OP hilf5=callocobject();
OP hilf6=callocobject(); OP hilf7=callocobject();
OP monom,monom1;
monom=a;
copy(cons_null,hilf6);
for (i=0L;i<S_V_LI(S_PO_S(monom));++i) 
{
  m_i_i(i+1L,hilf); 
  mult_apply(S_V_I(S_PO_S(monom),i),hilf);
  add_apply(hilf,hilf6);
}
monom1=b;
copy(cons_null,hilf7);
for (i=0L;i<S_V_LI(S_PO_S(monom1));++i) 
{
  m_i_i(i+1L,hilf); 
  mult_apply(S_V_I(S_PO_S(monom1),i),hilf);
  add_apply(hilf,hilf7);
}
if (lt(hilf6,hilf7)) return m_i_i(0L,c);
copy(cons_null,c);
while (monom!=NULL)
{
  monom1=b;
  while (monom1!=NULL)
  {
    while (nullp(S_PO_SI(monom1,S_V_LI(S_PO_S(monom1))-1L))) dec(S_PO_S(monom1));
    if (!emptyp(hilf)) freeself(hilf);
    sum_vector(S_PO_S(monom1),hilf);
    for (k=1L;k<=S_I_I(hilf);++k)
    {
      M_I_I(k,hilf1);
      copy(cons_null,hilf5);
      first_unordered_part_into_atmost_k_parts(k,S_V_LI(S_PO_S(monom1)),hilf2);
      do
      {
        copy(cons_eins,hilf3);
        for (l=0L;l<S_V_LI(hilf2);++l)
        {
          binom(S_V_I(S_PO_S(monom1),l),S_V_I(hilf2,l),hilf4);
          if (nullp(hilf4)) goto ref281094;
          mult_apply(hilf4,hilf3);
        }
        for (l=0L;l<S_V_LI(S_PO_S(monom));++l)
        {
          if (!nullp(S_V_I(S_PO_S(monom),l)))
          {
            sum_ueber_alle_teiler(l+1L,hilf2,hilf4);
            if (nullp(hilf4)) goto ref281094;
            hoch(hilf4,S_V_I(S_PO_S(monom),l),hilf4);
            mult_apply(hilf4,hilf3);
          }
        }
        add_apply(hilf3,hilf5);
        ref281094: 
        l=next_unordered_part_into_atmost_k_parts(hilf2);
      } while (l==1L);
      mult_apply(S_PO_K(monom),hilf5);
      mult_apply(S_PO_K(monom1),hilf5);
      sub(hilf,hilf1,hilf1);
      if (evenp(hilf1)) add_apply(hilf5,c);
      else sub(c,hilf5,c);
    }
    monom1=S_PO_N(monom1);
  }
  monom=S_PO_N(monom);
}
freeall(hilf); freeall(hilf1); freeall(hilf2); freeall(hilf3); 
freeall(hilf4); freeall(hilf5); freeall(hilf6); freeall(hilf7);
}

INT CF_right_inj(a,b,c)
OP a,b,c;
/* a Zyklenzeiger der Aktion auf Definitionsbereich
b Anzahl der Elemente im Bildbereich
c Anzahl der Bahnen injektiver Funktionen  */
{
OP hilf=callocobject();
OP hilf1=callocobject();
degree_from_zykelind(a,hilf);
if (lt(b,hilf)) 
{
  m_i_i(0L,c); goto frip_ende;
}
coeff_of_in(hilf,a,hilf1);
fakul(hilf,c);
mult_apply(hilf1,c);
binom(b,hilf,hilf1);
mult_apply(hilf1,c);
frip_ende: freeall(hilf);freeall(hilf1);
}

INT CF_left_inj(a,b,c)
OP a,b,c;
/* a = Zyklenzeiger von Gruppenaktion auf Bildbereich
b = Maechtigkeit des Definitionsbereichs
c = Anzahl der Klassen von Funktionen  */
{
INT i,j;
OP monom;
OP hilf=callocobject();
OP hilf1=callocobject();
OP hilf2=callocobject();
m_i_i(0L,c);
degree_from_zykelind(a,hilf);
if (lt(hilf,b)) goto frip_ende;
j=S_I_I(hilf);
for (i=S_I_I(b),copy(b,hilf);i<=j;++i,inc(hilf))
{
  monom=a;
  m_i_i(0L,hilf1);
  while (monom!=NULL)
  {
    if (S_PO_SII(monom,0L)==i) add_apply(S_PO_K(monom),hilf1);
    monom=S_PO_N(monom);
  }
  binom(hilf,b,hilf2);
  mult_apply(hilf2,hilf1);
  add_apply(hilf1,c);
}
fakul(b,hilf1);
mult_apply(hilf1,c);
frip_ende: freeall(hilf); freeall(hilf1); freeall(hilf2);
}

INT CF_left_right_inj(a,b,c)
OP a,b,c;
/* a = Gruppe, die auf dem Definitionsbereich operiert
   b = Gruppe, die auf dem Bildbereich operiert
   c = Ergebnis
*/
{
INT i;
OP monom,monom1;
OP hilf=callocobject();
OP hilf1=callocobject();
OP hilf2=callocobject();
degree_from_zykelind(a,hilf);
degree_from_zykelind(b,hilf1);
m_i_i(0L,c);
if (GT(hilf,hilf1)) goto frip_ende;
monom=a;
while (monom!=NULL)
{
  m_i_i(1L,hilf1);
  for (i=0L,m_i_i(1L,hilf);i<S_PO_SLI(monom);++i,inc(hilf))
  {
    if (!nullp(S_PO_SI(monom,i)))
    {
      hoch(hilf,S_PO_SI(monom,i),hilf2);
      mult_apply(hilf2,hilf1);
      fakul(S_PO_SI(monom,i),hilf2);
      mult_apply(hilf2,hilf1);
    }
  }
  monom1=b;
  m_i_i(0L,hilf2);
  while (monom1!=NULL)
  {
  /*println(S_PO_S(monom));
  println(S_PO_S(monom1));
  printf("%d\n",comp_le(S_PO_S(monom),S_PO_S(monom1)));*/
    if (comp_le(S_PO_S(monom),S_PO_S(monom1))==0L)
    {
      prod_binom(S_PO_S(monom1),S_PO_S(monom),hilf);
      /*println(hilf);*/
      /*mult_apply(hilf1,hilf);*/
    }
    else m_i_i(0L,hilf);
    mult_apply(S_PO_K(monom1),hilf);
    add_apply(hilf,hilf2);
    monom1=S_PO_N(monom1);
  }
  mult_apply(hilf1,hilf2);
  mult_apply(S_PO_K(monom),hilf2);
  add_apply(hilf2,c);
  monom=S_PO_N(monom);
}
frip_ende: freeall(hilf); freeall(hilf1); freeall(hilf2);
}

/* ***************************************************************

Various routines for computing the generators of some permutation 
groups. Furthermore the generators for induced actions can be computed.
In addition to these routines there is make_perm_paareperm in perm.c

******************************************************************** */

INT gen_arb(a)
OP a;
{
	INT i,j;
	OP hilf=callocobject();
	OP hil=callocobject();
	OP hilf1;
	printeingabe("Number of generators = ?");
	scan(INTEGER,hilf);
	m_l_v(hilf,a);
	printeingabe("All generators must have the same length!");
	printeingabe("Please enter the length of a permutation");
	scan(INTEGER,hil);
	for (i=0L;i<S_I_I(hilf);++i) {
		m_il_p(S_I_I(hil),S_V_I(a,i));
		hilf1=S_V_I(a,i);
		printf("Please enter the %d. generator \n",i+1L);
		for (j=0L;j<S_I_I(hil);++j) scan(INTEGER,S_P_I(hilf1,j));
	}
	freeall(hilf);
	freeall(hil);
	return(OK);
}

INT gen_Sn(a,b)
OP a,b;
/* b ist ein Erzeugendensystem von S_a. a=INTEGER, b=VECTOR.  */
{
	INT i;
	OP hilf;
	if (S_I_I(a)>2L) { /* two generators */
		m_il_v(2L,b);
		for (i=0L;i<2L;++i)
			m_il_p(S_I_I(a),S_V_I(b,i));
		hilf=S_V_I(b,0L); /* transposition */
  		M_I_I(2L,S_P_I(hilf,0L));
  		M_I_I(1L,S_P_I(hilf,1L));
  		for (i=2L;i<S_P_LI(hilf);++i) 
			M_I_I(i+1L,S_P_I(hilf,i));
		hilf=S_V_I(b,1L);
	}
	else { /* one generator */
		m_il_v(1L,b);
		m_il_p(S_I_I(a),S_V_I(b,0L));
		hilf=S_V_I(b,0L);
	}
	for (i=0L;i<S_I_I(a)-1L;++i)
		M_I_I(i+2L,S_P_I(hilf,i));
	M_I_I(1L,S_P_I(hilf,S_I_I(a)-1L)); /* a long cycle */
}

INT gen_Cn(a,b)
OP a,b;
/* b ist Erzeugendensystem von C_a  */
{
INT i;
OP hilf=callocobject();
m_il_v(1L,b);
m_il_p(S_I_I(a),hilf);
for (i=0L;i<S_P_LI(hilf)-1L;++i) M_I_I(i+2L,S_P_I(hilf,i));
M_I_I(1L,S_P_I(hilf,S_P_LI(hilf)-1L));
copy(hilf,S_V_I(b,0L));
freeall(hilf);
}

INT gen_Dn(a,b)
OP a,b;
/* b ist Erzeugendensystem von D_a  */
{
INT i;
OP hilf=callocobject();
m_il_v(2L,b);
m_il_p(S_I_I(a),hilf);
for (i=0L;i<S_P_LI(hilf)-1L;++i) M_I_I(i+2L,S_P_I(hilf,i));
M_I_I(1L,S_P_I(hilf,S_P_LI(hilf)-1L));
copy(hilf,S_V_I(b,0L));
for (i=0L;i<S_P_LI(hilf)-1L;++i) M_I_I(S_P_LI(hilf)-i,S_P_I(hilf,i));
copy(hilf,S_V_I(b,1L));
freeall(hilf);
}

INT gen_In(a,b)
OP a,b;
/* b ist Erzeugendensystem von I_a Identitaet auf a-elementiger Menge */
{
INT i;
OP hilf=callocobject();
m_il_v(1L,b);
m_il_p(S_I_I(a),hilf);
for (i=0L;i<S_P_LI(hilf);++i) M_I_I(i+1L,S_P_I(hilf,i));
copy(hilf,S_V_I(b,0L));
freeall(hilf);
}

INT gen_An(a,b)
OP a,b;
/* b ist Erzeugendensystem von A_a  */
{
INT i;
m_il_v(2L,b);
for (i=0L;i<2L;++i) m_il_p(S_I_I(a),S_V_I(b,i));
if (evenp(a))
{
  M_I_I(1L,S_P_I(S_V_I(b,0L),0L));
  for (i=1L;i<S_I_I(a)-1L;++i) M_I_I(i+2L,S_P_I(S_V_I(b,0L),i));
  M_I_I(2L,S_P_I(S_V_I(b,0L),S_I_I(a)-1L));
  M_I_I(3L,S_P_I(S_V_I(b,1L),0L));
  M_I_I(2L,S_P_I(S_V_I(b,1L),1L));
  for (i=2L;i<S_I_I(a)-1L;++i) M_I_I(i+2L,S_P_I(S_V_I(b,1L),i));
  M_I_I(1L,S_P_I(S_V_I(b,1L),S_I_I(a)-1L));
}
else
{
  for (i=0L;i<S_I_I(a)-1L;++i) M_I_I(i+2L,S_P_I(S_V_I(b,0L),i));
  M_I_I(1L,S_P_I(S_V_I(b,0L),S_I_I(a)-1L));
  M_I_I(3L,S_P_I(S_V_I(b,1L),0L));
  M_I_I(1L,S_P_I(S_V_I(b,1L),1L));
  for (i=2L;i<S_I_I(a)-1L;++i) M_I_I(i+2L,S_P_I(S_V_I(b,1L),i));
  M_I_I(2L,S_P_I(S_V_I(b,1L),S_I_I(a)-1L));
}
}

INT gen_An_3(a,b)
OP a,b;
/* b ist Erzeugendensystem von A_a bestehend aus allen 3-Zyklen  */
{
INT i,j,k,l;
OP hilf=callocobject();
if (le(a,cons_zwei))
{
m_il_v(1L,b);
m_il_p(S_I_I(a),hilf);
for (i=0L;i<S_P_LI(hilf);++i) M_I_I(i+1L,S_P_I(hilf,i));
copy(hilf,S_V_I(b,0L));
}
else
{
m_il_v(0L,b);
m_il_p(S_I_I(a),hilf);
for (i=0L;i<S_I_I(a)-2L;++i)
{
  for (j=i+1L;j<S_I_I(a)-1L;++j)
  {
    for (k=j+1L;k<S_I_I(a);++k)
    {
      for (l=0L;l<S_P_LI(hilf);++l)
      {
        if (l==i) M_I_I(j+1L,S_P_I(hilf,l));
        else if (l==j) M_I_I(k+1L,S_P_I(hilf,l));
        else if (l==k) M_I_I(i+1L,S_P_I(hilf,l));
        else M_I_I(l+1L,S_P_I(hilf,l));
      }
      inc(b);
      copy(hilf,S_V_I(b,S_V_LI(b)-1L));
    }
  }
}
}
freeall(hilf);
}

INT gen_M24(a)
OP a;
{
FILE *fp;
fp=fopen("data/m24.gen","r");
objectread(fp,a);
fclose(fp);
}

INT dir_sum_perm(a,b,c,d)
OP a,b,c,d;
{
INT i;
INT la=S_P_LI(a);
INT lb=S_P_LI(b);
m_il_p(la+lb,c);
copy(c,d);
for (i=0L;i<la;++i) 
{
  copy(S_P_I(a,i),S_P_I(c,i));
  M_I_I(i+1L,S_P_I(d,i));
}
for (;i<S_P_LI(c);++i) 
{
  M_I_I(i+1L,S_P_I(c,i));
  M_I_I(S_P_II(b,i-la)+la,S_P_I(d,i));
}
}

INT gen_dir_sum(a,b,c)
OP a,b,c;
{
if (emptyp(a)) copy(b,c);
else if (emptyp(b)) copy(a,c);
else
{
INT i,j;
OP speicher=callocobject();
OP hilf=callocobject();
OP hilf1=callocobject();
OP ax_e=callocobject();
init(BINTREE,speicher);
for (i=0L;i<S_V_LI(a);++i)
{
  for (j=0L;j<S_V_LI(b);++j)
  {
    dir_sum_perm(S_V_I(a,i),S_V_I(b,j),hilf,hilf1);
    ax_e=callocobject();copy(hilf,ax_e);
    insert(ax_e,speicher,NULL,NULL);
    ax_e=callocobject();copy(hilf1,ax_e);
    insert(ax_e,speicher,NULL,NULL);
  }
}
t_BINTREE_LIST(speicher,speicher);
t_LIST_VECTOR(speicher,c);
freeall(speicher);freeall(hilf);freeall(hilf1);
}
}

INT dir_prod_perm(a,b,c,d)
OP a,b,c,d;
{
INT i,j;
INT k=-1L;
INT la=S_P_LI(a);
INT lb=S_P_LI(b);
m_il_p(la*lb,c);
copy(c,d);
for (i=0L;i<la;++i) 
{
  for (j=0L;j<lb;++j)
  {
    ++k;
    M_I_I((S_P_II(a,i)-1L)*lb+j+1L,S_P_I(c,k));
    M_I_I(i*lb+S_P_II(b,j),S_P_I(d,k));
  }
}
}

INT gen_dir_prod(a,b,c)
OP a,b,c;
{
INT i,j;
OP speicher=callocobject();
OP hilf=callocobject();
OP hilf1=callocobject();
OP ax_e=callocobject();
init(BINTREE,speicher);
for (i=0L;i<S_V_LI(a);++i)
{
  for (j=0L;j<S_V_LI(b);++j)
  {
    dir_prod_perm(S_V_I(a,i),S_V_I(b,j),hilf,hilf1);
    ax_e=callocobject();copy(hilf,ax_e);
    insert(ax_e,speicher,NULL,NULL);
    ax_e=callocobject();copy(hilf1,ax_e);
    insert(ax_e,speicher,NULL,NULL);
  }
}
t_BINTREE_LIST(speicher,speicher);
t_LIST_VECTOR(speicher,c);
freeall(speicher);freeall(hilf);freeall(hilf1);
}

INT gen_on2sets(a,b)
OP a,b;
/* a ein Erzeugendensystem einer Permutationsgruppe
   b ein Erzeugendensystem der von a induzierten Permutationsgruppe
     auf der Menge aller 2-elementigen Mengen. */
{
INT i;
INT erg=OK;
if (generator_test(a)!=OK) 
   return error("gen_on2sets(a,b) a not a system of generators");
if (!emptyp(b)) freeself(b);
m_l_v(S_V_L(a),b);
for (i=0L;i<S_V_LI(a);++i) 
  /*m_perm_paareperm(S_V_I(a,i),S_V_I(b,i));*/
  m_perm_2setsperm(S_V_I(a,i),S_V_I(b,i));
}

INT generator_test(a)
OP a;
/* Testet ob a ein VECTOR von PERMUTATIONEN gleichen Grades ist. */
{
INT i,l;
if (S_O_K(a)!=VECTOR) return error("generator_test(a) a not VECTOR");
for (i=0L;i<S_V_LI(a);++i)
{
  if (S_O_K(S_V_I(a,i))!=PERMUTATION) 
     return error("generator_test(a) entries of a not PERMUTATIONs");
}
if (S_V_LI(a)>=1L) l=S_P_LI(S_V_I(a,0L));
for (i=1L;i<S_V_LI(a);++i)
{
  if (S_P_LI(S_V_I(a,i))!=l) 
     return error("generator_test(a) entries of a are not of the same degree");
}
return OK;
}

/* **************************************************************

Some special routines for handling wreath products KRANZ
and for the permutation representation in form of the
plethysm and exponentiation.

***************************************************************** */

INT t_kranz_plethysm(a,b)
OP a,b;
{
INT i,j,k=-1L;
INT n=S_P_LI(s_kr_g(a));
INT m=S_P_LI(s_kr_i(a,0L));
m_il_p(n*m,b);
for (j=0L;j<n;++j)
{
  for (i=0L;i<m;++i)
  {
    ++k;
    /*M_I_I((S_P_II(s_kr_g(a),j)-1L)*m+S_P_II(s_kr_i(a,S_P_II(s_kr_g(a),j)-1L),i),S_P_I(b,k));*/
    M_I_I((s_kr_gii(a,j)-1L)*m+s_kr_iji(a,s_kr_gii(a,j)-1L,i),S_P_I(b,k));
  }
}
}

OP s_kr_gi(a,i)
OP a;
INT i;
{
return(s_p_i(s_kr_g(a),i));
}

OP S_KR_GI(a,i)
OP a;
INT i;
{
return(S_P_I(s_kr_g(a),i));
}

INT s_kr_gii(a,i)
OP a;
INT i;
{
return(s_p_ii(s_kr_g(a),i));
}

INT S_KR_GII(a,i)
OP a;
INT i;
{
return(S_P_II(s_kr_g(a),i));
}

OP s_kr_ij(a,i,j)
OP a;
INT i,j;
{
return(s_p_i(s_kr_i(a,i),j));
}

OP S_KR_IJ(a,i,j)
OP a;
INT i,j;
{
return(S_P_I(s_kr_i(a,i),j));
}

INT s_kr_iji(a,i,j)
OP a;
INT i,j;
{
return(s_p_ii(s_kr_i(a,i),j));
}

INT S_KR_IJI(a,i,j)
OP a;
INT i,j;
{
return(S_P_II(s_kr_i(a,i),j));
}

INT gen_full_mon_group(a,b,c)
/* c ist ein VECTOR Objekt,
   die Eintragungen von c sind KRANZ Objekte
   c stellt ein Erzeugendensystem von G_a wr G_b dar
   wobei G_a und G_b Gruppen sind, die von a bzw. b (VECTORen)
   erzeugt werden. */
OP a,b,c;
{
INT la=S_P_LI(S_V_I(a,0L));
INT lb=S_P_LI(S_V_I(b,0L));
INT i,j,k,l;
OP einsa=callocobject();
OP einsb=callocobject();
OP borbits=callocobject();
m_il_v(0L,c);
k=-1L;
first_permutation(S_P_L(S_V_I(a,0L)),einsa);
first_permutation(S_P_L(S_V_I(b,0L)),einsb);
for (i=0L;i<S_V_LI(b);++i)
{
  k++;
  inc(c);
  init_kranz(S_V_I(c,k));
  m_il_v(lb,s_kr_v(S_V_I(c,k)));
  copy(S_V_I(b,i),s_kr_g(S_V_I(c,k)));
  for (j=0L;j<lb;++j) copy(einsa,s_kr_i(S_V_I(c,k),j));
}
allorbits(b,borbits);
for (l=0L;l<S_V_LI(borbits);++l)
{
  for (i=0L;i<S_V_LI(a);++i)
  {
    ++k;
    inc(c);
    init_kranz(S_V_I(c,k));
    m_il_v(lb,s_kr_v(S_V_I(c,k)));
    copy(einsb,s_kr_g(S_V_I(c,k)));
    for (j=0L;j<S_V_II(borbits,l)-1L;++j) 
	copy(einsa,s_kr_i(S_V_I(c,k),j));
    copy(S_V_I(a,i),s_kr_i(S_V_I(c,k),j));
    for (++j;j<lb;++j) 
	copy(einsa,s_kr_i(S_V_I(c,k),j));
  }
}
freeall(einsa); freeall(einsb); freeall(borbits);
}

INT gen_plethysm(a,b,c)
OP a,b,c;
/* c ist ein Erzeugendensystem des Kranzproduktes Ga wr Gb (operierend
in Form des Plethysmus) wobei Ga und Gb von a bzw. b erzeugt werden. */
{
INT i;
OP d=callocobject();
gen_full_mon_group(a,b,d);
m_l_v(S_V_L(d),c);
for (i=0L;i<S_V_LI(d);++i)
    t_kranz_plethysm(S_V_I(d,i),S_V_I(c,i));
freeall(d);
}

INT gen_centralizer_permutation(a,b)
OP a,b;
/* b ist ein Erzeugendensystem einer zu dem Stabilisator einer 
   Permutation (gegeben in Form ihrer zykelpartition (PARTITION))
   isomorphen Gruppe.
*/
{
INT i;
INT anfang=0L;
OP ga=callocobject();
OP gb=callocobject();
OP gc=callocobject();
OP gd=callocobject();
OP hilf=callocobject();
for (i=0L,M_I_I(1L,hilf);i<S_PA_LI(a);++i,inc(hilf))
{
  if (!nullp(S_PA_I(a,i)))
  { 
    gen_Cn(hilf,ga);
    gen_Sn(S_PA_I(a,i),gb);
    if (anfang==0L)
    {
      gen_plethysm(ga,gb,b);
      anfang=1L;
    }
    else
    {
      copy(b,gc);
      gen_plethysm(ga,gb,gd);
      gen_dir_sum(gc,gd,b);
    }
  }
}
freeall(ga); freeall(gb); freeall(gc); freeall(gd); freeall(hilf);
}

INT gen_stabilizer_partition(a,b)
OP a,b;
/* b ist ein Erzeugendensystem einer zu dem Stabilisator einer 
   Partition (von Typ EXPONENT) isomorphen Gruppe.
*/
{
INT i;
INT anfang=0L;
OP ga=callocobject();
OP gb=callocobject();
OP gc=callocobject();
OP gd=callocobject();
OP hilf=callocobject();
OP aa=callocobject();
if (S_PA_K(a)==VECTOR) t_VECTOR_EXPONENT(a,aa);
else copy(a,aa);
for (i=0L,M_I_I(1L,hilf);i<S_PA_LI(aa);++i,inc(hilf))
{
  if (!nullp(S_PA_I(aa,i)))
  { 
    gen_Sn(hilf,ga);
    gen_Sn(S_PA_I(aa,i),gb);
    if (anfang==0L)
    {
      gen_plethysm(ga,gb,b);
      anfang=1L;
    }
    else
    {
      copy(b,gc);
      gen_plethysm(ga,gb,gd);
      gen_dir_sum(gc,gd,b);
    }
  }
}
freeall(ga); freeall(gb); freeall(gc); freeall(gd); freeall(hilf); freeall(aa);
}

INT gen_young_group(a,b)
OP a,b;
{
OP aa,hilf,hilf1,hilf2;
INT i,j;
if (S_O_K(a)!=PARTITION) error("gen_young_group(a,b) a not a PARTITION ");
if (!emptyp(b)) freeself(b);
aa=callocobject();
hilf=callocobject();
hilf1=callocobject();
hilf2=callocobject();
if (S_PA_K(a)!=EXPONENT) t_VECTOR_EXPONENT(a,aa); 
else copy(a,aa);
for (i=0L,M_I_I(1L,hilf);i<S_PA_LI(aa);++i,inc(hilf))
{
  if (!nullp(S_PA_I(aa,i)))
  {
    gen_Sn(hilf,hilf1);
    for (j=1L;j<=S_PA_II(aa,i);++j) 
    {
      gen_dir_sum(b,hilf1,hilf2);
      copy(hilf2,b);
    }
  }
}
freeall(aa); freeall(hilf); freeall(hilf1); freeall(hilf2);
}

/* ***************************************************************

Various cycle index formulae. The standard cycle index formulae
can be found in zyk.c

******************************************************************** */

INT zykelind_In(a,b)
OP a,b;
{
INT erg=OK;
if (S_O_K(a)!=INTEGER) return error("zykelind_In(a,b) a not INTEGER");
if (S_I_I(a)<1L) return error("zykelind_In(a,b) a<1");
if (!emptyp(b)) erg+=freeself(b);
erg+=m_iindex_iexponent_monom(0L,s_i_i(a),b);
return(erg);
} 

INT degree_from_zykelind(a,b)
OP a,b;
/* Bestimmt den Grad eine Gruppenaktion aus deren Zyklenzeiger */
{
INT i;
OP hilf=callocobject();
OP hilf1=callocobject();
m_i_i(0L,b);
for (i=0L,m_i_i(1L,hilf);i<S_V_LI(S_PO_S(a));++i,inc(hilf))
{
  if (!nullp(S_PO_SI(a,i)))
  {
    mult(S_PO_SI(a,i),hilf,hilf1);
    add_apply(hilf1,b);
  }
}
freeall(hilf); freeall(hilf1);
}

INT zykelind_M24(a)
OP a;
{
FILE *fp;
fp=fopen("data/m24.zyk","r");
objectread(fp,a);
fclose(fp);
}

/* *****************************************************************

Combinatorial numbers and functions in SYMMETRICA:
The Bell numbers: both recursive and non recursive methods.

******************************************************************* */

INT Bell_numbers(a,b)
OP a,b;
/* a ein INTEGER Objekt
   b ist ein VECTOR Objekt, S_V_I(b,i) ist B(i) die i-te Bell Zahl
*/
{
OP hilf=callocobject();
INT i,j;
if (S_I_I(a)<0L) return error(" Bell_numbers(a,b)  a<0 ");
m_il_v(S_I_I(a)+1L,b);
stirling_second_tafel(a,hilf);
for (i=0L;i<=S_I_I(a);++i) 
{
  copy(cons_null,S_V_I(b,i));
  for (j=0L;j<=i;++j) add_apply(S_M_IJ(hilf,i,j),S_V_I(b,i));
}
freeall(hilf);
}

INT Bell_number(a,b)
OP a,b;
/* a ein INTEGER Objekt
   b Ergebnis die a-te Bell Zahl =S(a,0)+S(a,1)+...+S(a,a), 
   wobei S(i,j) die zweiten Stirling Zahlen sind.
*/
{
OP hilf=callocobject();
INT i;
if (S_I_I(a)<0L) return error(" Bell_number(a,b)  a<0 ");
copy(cons_null,b);
stirling_second_tafel(a,hilf);
for (i=0L;i<=S_I_I(a);++i) add_apply(S_M_IJ(hilf,S_I_I(a),i),b);
freeall(hilf);
}

INT Bell_number_recursive(a,b)
OP a,b;
{
OP rt=callocobject();
if (S_I_I(a)<0L) return error(" Bell_number_recursive(a,b)  a<0 ");
m_il_v(0L,rt);
bellrec(a,b,rt);
/*println(rt);*/
freeall(rt);
}

static INT bellrec(a,b,rt)
OP a,b,rt;
{
INT i;
if (le(a,S_V_L(rt))) return copy(S_V_I(rt,S_I_I(a)-1L),b);
else
{
  OP aa=callocobject();
  OP hilf=callocobject();
  OP hilf1=callocobject();
  OP hilf2=callocobject();
  copy(a,aa);dec(aa);
  m_i_i(1L,b);
  for (i=0L,copy(aa,hilf);i<S_I_I(a)-1L;++i,dec(hilf))
  {
    binom(aa,hilf,hilf1);
    bellrec(hilf,hilf2,rt);
    mult_apply(hilf2,hilf1);
    add_apply(hilf1,b);
  }
  inc(rt);
  copy(b,S_V_I(rt,S_V_LI(rt)-1L));
  freeall(hilf); freeall(hilf1); freeall(hilf2);freeall(aa);
}
}

/* ***************************************************************

Various methods of substitutions into a cycle index.
Furthermore there is polya_sub and polya_n_sub.
Special routines are available for multi-dimensional cycle indices.
See in the section where the cycle indices of the platonic solids
and the fulleren are computed.

****************************************************************** */

INT polya_compl_sub(a,b)
OP a,b;
/* einsetzung */ /* a ist polynom */ 
/* b wird ergebnis x_i ----> 0 falls i=1,3,5,7,9... und 2 falls i=2,4,6,8,...
*/
{
OP c,d;
INT i;
INT erg=OK;
if (S_O_K(a)!=POLYNOM) return error("polya_compl_sub(a,b) a not POLYNOM");
if (not EMPTYP(b)) erg+=freeself(b);
c=callocobject();
d=callocobject();
numberofvariables(a,c);
m_l_v(c,d);
for (i=0L;i<S_I_I(c);++i)
{
  if ((i+1L)%2L==1L) M_I_I(0L,S_V_I(d,i));
  else M_I_I(2L,S_V_I(d,i));
}
erg += eval_polynom(a,d,b); 
erg += freeall(c);
erg += freeall(d);
if (erg != OK)
return error("polya_compl_sub: error during computation");
return erg;
}

/* ****************************************************************

Some useful routines for zykelind.c 

****************************************************************** */

static INT sum_ueber_alle_teiler(i,a,b)
INT i;
OP a,b;
/* i ist eine INT Zahl
   a ist ein VECTOR Objekt
   b ist die Summe ueber a_j wobei j alle Teiler von i durchlaeuft
*/
{
INT j;
OP ii=callocobject();
OP hilf=callocobject();
OP hilf1=callocobject();
copy(cons_null,b);
M_I_I(i,ii);
alle_teiler(ii,hilf);
for (j=0L;j<S_V_LI(hilf);++j)
{
  if (S_V_II(hilf,j)<=S_V_LI(a)) 
  {
    mult(S_V_I(hilf,j),S_V_I(a,S_V_II(hilf,j)-1L),hilf1);
    add_apply(hilf1,b);
  }
}
freeall(ii);freeall(hilf);freeall(hilf1);
}

INT comp_le(a,b)
OP a,b;
/* a,b sind PARTITIONen (oder VECTORen) vom Typ EXPONENT
   comp_le ==0 <=> a[i]<=b[i] fuer alle
   comp_le ==1 sonst
*/
{
INT i;
if (S_O_K(a)==PARTITION) 
{
  for (i=0L;(i<S_PA_LI(a)) && (i<S_PA_LI(b)); ++i)
  {
    if (S_PA_II(a,i) > S_PA_II(b,i)) return(1L);
  }
  for (;i<S_PA_LI(a);++i)
  {
    if (S_PA_II(a,i)  != 0 ) return(1L);
  }
  return(0L);
}
else if (S_O_K(a)==VECTOR)
{
  for (i=0L;(i<S_V_LI(a)) && (i<S_V_LI(b)); ++i)
  {
    if (S_V_II(a,i)> S_V_II(b,i)) return(1L);
  }
  if (S_V_LI(a)>S_V_LI(b))
  {
    for (;i<S_V_LI(a);++i)
    {
      if (S_V_II(a,i) != 0 ) return(1L);
    }
  }
  return(0L);
}
}

OP binomspeicher[50][50];
INT prod_binom(a,b,c)
OP a,b,c;
/* a,b sin PARTITIONen vom Typ EXPONENT
   c = prod_i {a[i] choose b[i]}
*/
{
INT i;
m_i_i(1L,c);
/*printobjectkind(a);*/
if (S_O_K(a)==PARTITION)
{
  for (i=0L; (i<S_PA_LI(a)) && (i<S_PA_LI(b)); ++i)
  {
    if (S_PA_II(b,i) != 0) 
    {     
        //if (S_PA_II(a,i) == 0) return m_i_i(0,c);
        ////else if (S_PA_II(a,i) < S_PA_II(b,i) ) return m_i_i(0,c);
      //else 
       if (S_PA_II(b,i) == 1)
	      mult_apply(S_PA_I(a,i),c);
	else {
		if ((S_PA_II(a,i)<50) && (S_PA_II(b,i)<50))
		{
		if (binomspeicher[S_PA_II(a,i)][S_PA_II(b,i)]==NULL)
			{
			binomspeicher[S_PA_II(a,i)][S_PA_II(b,i)]=callocobject();
		        binom(S_PA_I(a,i),S_PA_I(b,i), binomspeicher[S_PA_II(a,i)][S_PA_II(b,i)]);
			}
	        mult_apply( binomspeicher[S_PA_II(a,i)][S_PA_II(b,i)] ,c);
		}
		else {
			OP hilf=callocobject();
		      binom(S_PA_I(a,i),S_PA_I(b,i),hilf);
		      mult_apply(hilf,c);
			freeall(hilf);
			}
		}
    }
  }
}
else if ((S_O_K(a)==VECTOR)||(S_O_K(a)==INTEGERVECTOR))
{
			OP hilf=callocobject();
  for (i=0L; (i<S_V_LI(a)) && (i<S_V_LI(b)); ++i)
  {
    if (S_V_II(b,i) != 0) 
    {
      binom(S_V_I(a,i),S_V_I(b,i),hilf);
      mult_apply(hilf,c);
    }
  }
			freeall(hilf);
}

//	print(a);printf(" ");
//	print(b);printf(" ");println(c);
}

static INT first_unordered_part_into_atmost_k_parts(s,k,a)
/*first_pa_into_atmost_k_parts(s,k,a)*/
/* Bestimmt die erste Darstellung von s als Summe von k Summanden >=0.
Diese Darstellung wird als VECTOR a weitergegeben.  */
OP a;
INT k,s;
{
int i;
INT erg=OK;
if (!emptyp(a)) freeself(a);
m_il_v(k,a);
for (i = 0L;i < k-1L; M_I_I(0L,S_V_I(a,i)) , ++i)
    ;
erg+=M_I_I(s,S_V_I(a,k-1L));
if (erg!=OK) error(" in computation of first_unordered_part_into_atmost_k_parts(s,k,a)  ");
return(erg);
}

static INT next_unordered_part_into_atmost_k_parts(a)
/*next_pa_into_atmost_k_parts(a)*/
OP a;
/* Berechnet den Nachfolger der Darstellung einer natuerlichen Zahl
als Summe von hoechstens k Summanden >=0.  Die natuerliche Zahl ist
dabei die Summe ueber alle Elemente von a, k ist die Laenge von a. */
{
int i;
if (S_O_K(a)!=VECTOR) return error("next_unordered_part_into_atmost_k_parts(a)  a not VECTOR");
i=s_v_li(a)-1L;
while( nullp(S_V_I(a,i)) && (i>=0L)) --i;
if (i<=0L)  return(2L); /* alle aufglistet */
copy(S_V_I(a,i),S_V_I(a,S_V_LI(a)-1L));
dec(S_V_I(a,S_V_LI(a)-1L));
inc(S_V_I(a,i-1L));
if (i<S_V_LI(a)-1L) M_I_I(0L,S_V_I(a,i));
return(1L); /*  kein Fehler aufgetreten  */
}