/* dw.c  */

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

INT fixed_point_id(length,weight,fp)
	OP length,weight,fp; 
	/* fp ist ein random VECTOR der Laenge lenght vom Gewicht weight 
	 * also eine 0,1 VECTOR mit weight Eintragungen "1", die zufaellig
	 * verteilt sind */
{
	INT i,l=S_I_I(length);
	OP hilf=callocobject();
	if (S_O_K(fp)!=VECTOR) m_il_nv(l,fp);
	else if (S_V_LI(fp)!=l) m_il_nv(l,fp);
	else for (i=0L;i<l;++i) 
		if (!nullp(S_V_I(fp,i))) m_i_i(0L,S_V_I(fp,i));
	i=0L;
	while (i<S_I_I(weight)) {
		random_integer(hilf,cons_null,length);
		if (S_V_II(fp,S_I_I(hilf))==0L) {
			M_I_I(1L,S_V_I(fp,S_I_I(hilf)));
			++i;
		}
	}
	freeall(hilf);
	return(OK);
}

INT fixed_point_nonid(perm,nfixp,weight,fp)
	OP perm,nfixp,weight,fp; 
	/* perm opreriert uf dem Definitionsbereich von {0,1}-Funktionen, also 
	 * von charakteristischen Funktionen.
	 * Berechnet einen random Fixpunt der Permutation perm vom Gewicht
	 * weight, der als VECTOR fp weitergegeben wird.
	 * Dabei ist nfixp die Anzahl der Fixpunkte vom Gewicht weight
	 * von perm */
{
	OP zyktyp=callocobject();
	OP zz=callocobject();
	OP hilf=callocobject();
	OP part=callocobject();
	OP part_stored=callocobject();
	OP zsum=callocobject();
	OP pz=callocobject();
	if (S_O_K(nfixp)==INTEGER) {
		random_integer(zz,cons_null,nfixp);
		inc(zz);
	}
	else random_longint(zz,nfixp);
	m_i_i(0L,zsum);
	zykeltyp(perm,zyktyp);
	t_VECTOR_EXPONENT(zyktyp,zyktyp);
	bestimme_pz(zyktyp,perm,pz);
	first_part_EXPONENT(weight,part);
	do {
		if (comp_le(part,zyktyp)==0L) {
			prod_binom(zyktyp,part,hilf);
			add_apply(hilf,zsum);
			copy(part,part_stored);
		}
	} while ((gt(zz,zsum)) && (next(part,part)));
	bestimme_fixpunkt(perm,part_stored,zyktyp,pz,fp);
	freeall(zyktyp);
	freeall(zz);
	freeall(zsum); 
	freeall(hilf);
	freeall(part); 
	freeall(part_stored); 
	freeall(pz); 
} 

static INT bestimme_pz(a,g,pz)
	OP a,g,pz;
/* RH 031092 */ /* Entnommen und leicht abgeaendert aus di.c HF */
/*  pz[0][i]=j  i+1 liegt in einem Zykel der Laenge j
    pz[1][i]=k  i+1 liegt im k-ten Zykel (der Laenge j) */
{

	INT weiter,anz,length,beginn,i;
	OP anz_length = callocobject();
	OP besucht = callocobject();
	m_il_v(2L,pz);
	for (i=0L;i<=1L;++i) m_il_nv(S_P_LI(g),S_V_I(pz,i));
	anz = length = 0L; 
	beginn = 1L;
	m_il_nv(S_PA_LI(a),besucht);
	copy(besucht,anz_length);
	while (anz < S_PA_LI(a)) {
		length = 0L;
		weiter = beginn;
		/* Laenge des Zykels von beginn-1 bestimmen */
		do {
			length++;
			weiter = S_P_II(g,weiter-1L);
			M_I_I(1L,S_V_I(besucht,weiter-1L));
		} while (weiter != beginn);
		/* Information fuer beginn-1 schreiben */
		M_I_I(1L+S_V_II(anz_length,length),S_V_I(anz_length,length));
		M_I_I(length,S_V_I(S_V_I(pz,0L),beginn-1L));
		M_I_I(S_V_II(anz_length,length),S_V_I(S_V_I(pz,1L),beginn-1L));
		/* Information fuer die anderen Elemente des Zykels von beginn-1 schreiben */
		weiter = beginn;
		do {
			weiter = S_P_II(g,weiter-1L);
			M_I_I(length,S_V_I(S_V_I(pz,0L),weiter-1L)); 
			M_I_I(S_V_II(anz_length,length),S_V_I(S_V_I(pz,1L),weiter-1L));
		} while(weiter != beginn);
		/* Ein Element eines noch nicht besuchten Zykels suchen */
		for(i=0L;i<S_PA_LI(a);++i) {
			if(S_V_II(besucht,i) == 0L) {
				beginn = i+1L;
				break; 
			}
		}
		anz+= length;
	}
	freeall(anz_length);
	freeall(besucht);
	return OK;
}

static INT bestimme_fixpunkt(perm,part,zyktyp,pz,fp,weight)
	OP perm,part,zyktyp,pz,fp,weight;
{
	INT i,j,k,l=S_P_LI(perm),jj; 
	OP hilf=callocobject(); 
	if (S_O_K(fp)!=VECTOR) m_il_nv(l,fp);
	else if (S_V_LI(fp)!=l) m_il_nv(l,fp);
	else for (i=0L;i<l;++i) 
		if (!nullp(S_V_I(fp,i))) m_i_i(0L,S_V_I(fp,i));
	for (i=0L;i<S_PA_LI(part);++i) {
		if (S_PA_II(part,i)!=0L) {
			if (S_PA_II(part,i)==S_PA_II(zyktyp,i)) {
				for (j=0L;j<S_P_LI(perm);++j) { 
					if (S_V_II(S_V_I(pz,0L),j)==i+1L) 
						M_I_I(1L,S_V_I(fp,S_P_II(perm,j)-1L)); 
				}
			}
			else { 
				k=0L;
				l=0L;
				while (k<S_PA_II(part,i)) {
					random_integer(hilf,cons_null,S_PA_I(zyktyp,i));
					inc(hilf);
					for (j=0L;j<S_P_LI(perm);++j) {
						if ((S_V_II(S_V_I(pz,0L),j)==i+1L) && (S_V_II(S_V_I(pz,1L),j)==S_I_I(hilf))) {
							jj=S_P_II(perm,j)-1L; 
							if (S_V_II(fp,jj)==0L) {
								M_I_I(1L,S_V_I(fp,jj));
								l=1L;
							}
						}
					}
					k=k+l; 
					l=0L;
				}
			}
		}
	}
	freeall(hilf); 
}

INT berechne_fixpunkt_anzahl(perm,weight,fanz)
	OP perm,weight,fanz; 
{
	OP poly;
	OP cycle_a=callocobject();
	OP a=callocobject();
	OP b=callocobject();
	INT erg=OK;
	/*m_l_nv(S_P_L(perm),a);*/
	erg+=zykeltyp(perm,cycle_a);
	t_VECTOR_EXPONENT(cycle_a,cycle_a);
	m_skn_po(S_PA_S(cycle_a),cons_eins,NULL,a);
	erg+=co_polya_sub(a,S_PA_L(cycle_a),weight,b);
	poly=b;
	while (S_PO_N(poly) != NULL)
		poly=S_PO_N(poly);
	if (eq(S_PO_SI(poly,0L),weight))
		erg+=copy(S_PO_K(poly),fanz);
	else m_i_i(0L,fanz);
	erg+=freeall(cycle_a);
	erg+=freeall(a);
	erg+=freeall(b);
	ENDR("berechne_fixpunkt_anzahl");
}