//Funktionsbezeichner (Hilfe)
proc fcthelp ()
{
print("tupel_summe (matrix t1, matrix t2)	// Gibt Vielfachensumme von 2 Tupeln zurueck");
print("tupel_sum (matrix t1, matrix t2)	// Gibt Summe von 2 Tupeln zurueck (als Tupel)");
print("tupel_diff (matrix t1, matrix t2)	// Gibt Differenz von 2 Tupeln zurueck (als Tupel)");
print("");
print("syzg (matrix Q, matrix G)	//Gibt zurueck, ob Q eine Syzygie von G ist");
print("s_poly (poly p, poly g)	//Gibt S-Polynom von p und g zurueck");
print("sij_poly (int i, int j, ideal F)	// Gibt S_ij zurueck (Tupel, bzw. Matrix)");
print("lt_tupel (ideal F)	// Gibt Tupel der Leitterme des -Ideals- F zurueck");
print("lt_tupel_2 (matrix F)	// Gibt Tupel der Leitterme des -Tupels- F zurueck");
print("ideal_tupel (ideal F)	// Gibt Tupel der Polynome des Ideals F zurueck");
print("");
print("max_mg (matrix Q, ideal F)	// Gibt den Multigrad von Q bez. F zurueck");
print("f_hom (matrix Q, ideal F)	// Prueft, ob Q F-Homogen ist und schreibt ggf. den Multigrad u aus");
print("");
print("grad_zu_monom (intvec v)	// Wandelt einen Multigrad in das zugehoerige Monom um");
print("ci_X (matrix S, int i, int j, ideal F)	// Gibt den Term zurück, mit dem S_ij multipliziert wird ( c_i*X^(u-v_ij) )");
print("zerlegung (matrix S, matrix S_ij, int i, int j, ideal F)	// Zerlegt S in S(i) = S(i-1) - u(i)");
print("s_zerlegung (matrix S, ideal F)	// Zerlegt S in S = u(1)+u(2)+...+u(k)");
print("kgv_tupel (matrix LTF)	// Bildet das kgV von allen Monomen im Tupel");
print("hom (ideal F)	// Gibt ein F-homogenes Tupel zurueck, welches Syzygie von LT(F) ist");
print("");
print("lambda (matrix Q, ideal F)	// Funktion lambda (Diagramm)");
print("glambda (matrix Q, ideal F)	// Funktion Glambda (Diagramm)");
print("M (matrix Q, ideal F)	// Funktion M (Diagramm)");
print("");
print("qij (int i, int j, ideal F)	// Bildet Tupel Q_ij, mit dem sich lambda(S_ij)=sum(q_i*f_i) speziell erzeugen laesst");
print("sij_urbild (int i, int j, ideal F)	// Bildung von Urbildern der S_ij bez. M");
print("poly_mal_tupel (poly p, matrix T)	// Multipliziert ein Polynom mit jedem Element eines Tupels");
print("s_urbild (matrix S, ideal F)	// Bildung von Urbildern der hom. Syz. von LT(F) bez. M");
print("s_test (int n, ideal F)	//Test mit n homogenen Syzygien von LT(F)");
};


// Gibt Vielfachensumme von 2 Tupeln zurueck
proc tupel_summe (matrix t1, matrix t2)
{
int m = size(t1);
poly s = 0;
for (int k=1; k<=m; k=k+1)
{
	s = s + t1[1,k]*t2[1,k];
}
return(s);
};


// Gibt Summe von 2 Tupeln zurueck (als Tupel)
proc tupel_sum (matrix t1, matrix t2)
{
int m = size(t1);
poly s = 0;
for (int k=1; k<=m; k=k+1)
{
	t1[1,k] = t1[1,k] + t2[1,k];
}
return(t1);
};


// Gibt Differenz von 2 Tupeln zurueck (als Tupel)
proc tupel_diff (matrix t1, matrix t2)
{
int m = size(t1);
poly s = 0;
for (int k=1; k<=m; k=k+1)
{
	t1[1,k] = t1[1,k] - t2[1,k];
}
return(t1);
};


//Gibt zurueck, ob Q eine Syzygie von G ist
proc syzg (matrix Q, matrix G)
{
poly sum = tupel_summe(Q,G);
if (sum == 0)
{
	printf("(%s) ist Syzygie von (%s)", Q, G);
} else
{
	
	printf("(%s) ist keine Syzygie von (%s)", Q, G);
	printf("Summe: %s", sum);
}
};


//Gibt S-Polynom von p und g zurueck
proc s_poly (poly p, poly g)
{
poly kgV = normalize(lcm(leadmonom(p),leadmonom(g)));
poly s = (kgV/leadmonom(p))*p - (leadcoef(p)/leadcoef(g))*(kgV/leadmonom(g))*g;
return(s);
};


// Gibt S_ij zurueck (Tupel, bzw. Matrix)
proc sij_poly (int i, int j, ideal F)
{
int m = size(F);
matrix S[1][m];
poly kgV = normalize(lcm(leadmonom(F[i]),leadmonom(F[j])));
for (int k=1; k<=m; k=k+1)
{
	S[1,k]=0;
	if (k == i) 
	{
		S[1,k]=(kgV/leadmonom(F[i]));
	}
	if (k == j)
	{
		S[1,k]= -(leadcoef(F[i])/leadcoef(F[j]))*(kgV/leadmonom(F[j]));
	}
}
return(S);
};


// Gibt Tupel der Leitterme des -Ideals- F zurueck
proc lt_tupel (ideal F)
{
int m = size(F);
matrix T[1][m];
for (int k=1; k<=m; k=k+1)
{
	T[1,k]=lead(F[k]);
}
return(T);
};


// Gibt Tupel der Leitterme des -Tupels- F zurueck
proc lt_tupel_2 (matrix F)
{
int m = size(F);
matrix T[1][m];
for (int k=1; k<=m; k=k+1)
{
	T[1,k]=lead(F[1,k]);
}
return(T);
};


// Gibt Tupel der Polynome des Ideals F zurueck
proc ideal_tupel (ideal F)
{
int m = size(F);
matrix T[1][m];
for (int k=1; k<=m; k=k+1)
{
	T[1,k]=F[k];
}
return(T);
};


// Gibt den Multigrad von Q bez. F zurueck
proc max_mg(matrix Q, ideal F)
{
matrix TF = ideal_tupel(F);
int m = size(Q);
poly s = 0;
intvec u = leadexp(0);
for (int k=1; k<=m; k=k+1)
{
	if (Q[1,k] != 0)
	{
		if (u < leadexp(Q[1,k]*TF[1,k])) 
		{
			u = leadexp(Q[1,k]*TF[1,k]);
		}
	}
}
return(u);
};


// Prueft, ob Q F-Homogen ist und schreibt ggf. den Multigrad u aus 
proc f_hom (matrix Q, ideal F)
{
matrix TF = ideal_tupel(F);
int m = size(Q);
poly s = 0;
intvec u = max_mg(Q,F);
for (int k=1; k<=m; k=k+1)
{
	if (Q[1,k] != 0)
	{
		if (Q[1,k] != lead(Q[1,k]))
		{
			printf("kein Term in Q an Stelle: %s", k);
			return(0);
		}
		if (u != leadexp(Q[1,k]*TF[1,k]))
		{
			printf("falscher Multigrad an Stelle: %s", k);
			return(0);
		}
	}
}

	printf("F-homogen vom Multigrad: %s", u);

return(1);
};


// Wandelt einen Multigrad in das zugehoerige Monom um
proc grad_zu_monom (intvec v)
{
int m = size(v);
poly q = 1;
for (int k=1; k<=m; k=k+1)
{
	q = q*var(k)^v[k];
}
return(q);
};


// Gibt den Term zurück, mit dem S_ij multipliziert wird ( c_i*X^(u-v_ij) )
proc ci_X (matrix S, int i, int j, ideal F)
{
poly kgV = normalize(lcm(leadmonom(F[i]),leadmonom(F[j])));
intvec u = leadexp(F[i]*S[1,i]);
number c_i = leadcoef(S[1,i]);
poly mon = grad_zu_monom(u)/kgV;
poly p = c_i*mon;
return(p);
}


// Zerlegt S in S(i) = S(i-1) - u(i)
proc zerlegung (matrix S, matrix S_ij, int i, int j, ideal F)
{
int m = size(S);
int l=1;
matrix tupel[1][m];
poly p = ci_X(S,i,j,F);

	printf("S(%s) = (%s)", i-1, S);
	printf("");

	for (l=1; l<=m; l=l+1)
	{
		tupel[1,l] = p*S_ij[1,l];
	}
	printf("     - (%s)   = u(%s) =  %s*S_%s%s = q_%s%s*S_%s%s", tupel, i, p, i, j, i, j, i, j);
	printf("");

S[1,i]= S[1,i] - p*S_ij[1,i];
S[1,j]= S[1,j] - p*S_ij[1,j];

	printf("     = (%s)   = S(%s)", S, i);
	printf("");
	printf("");

return(S);
};


// Zerlegt S in S = u(1)+u(2)+...+u(k)
proc s_zerlegung (matrix S, ideal F)
{
printf("-->S-Zerlegung------------------------------------start");
print("");
int m = size(S);
matrix t;
for (int k=1; k<=m-1; k=k+1)
{
	S = zerlegung(S, sij_poly(k,m,F), k, m, F);
}
for (k=1; k<=m-1; k=k+1)
{
	t = sij_poly(k,m,F);
	printf("S_%s%s = (%s)", k, m, t);
}
print("");
string str = "S = q_"+string(1)+string(m)+"*S_"+string(1)+string(m);
for (k=2; k<=m-1; k=k+1)
{
	str = str + " + q_"+string(k)+string(m)+"*S_"+string(k)+string(m);
} 
print(str);
print("");
printf("-->S-Zerlegung--------------------------------------end");
print("");
};


// Bildet das kgV von allen (Leit-)Monomen im Tupel
proc kgv_tupel (matrix LTF)
{
LTF = lt_tupel_2(LTF);
int s = size(LTF);
poly r = LTF[1,s];
for (int k=1; k<=s-1; k=k+1)
{
	r = normalize(lcm(r,LTF[1,k]));
}
return(r);
};


// Gibt ein F-homogenes Tupel zurück, welches Syzygie von LT(F) ist
proc hom (ideal F)
{
matrix LTF = lt_tupel(F);
intvec u = leadexp(kgv_tupel(LTF));
int s = size(u);
int t = size(LTF);
matrix S[1][t];
int rand = 0;
number a = 0;
for (int k=1; k<=s; k=k+1)
{
	u[k] = u[k]+random(0,5);
}
for (k=1; k<=t; k=k+1)
{
	S[1,k] = grad_zu_monom(u - leadexp(LTF[1,k]));
}
for (k=1; k<=t; k=k+1)
{
	if (k != t) 
	{
		rand = random(-10,10);
		a = leadcoef(F[k])*rand + a;
		S[1,k] = S[1,k]*rand;
	} else
	{
		S[1,k] = S[1,k]*(-a/leadcoef(F[k]));
	}
}
return(S);
};


// Funktion lambda (Diagramm)
proc lambda (matrix Q, ideal F)
{
poly s = tupel_summe(Q, ideal_tupel(F));
return(s);
};

// Funktion glambda (Diagramm)
proc glambda (matrix Q, ideal F)
{
poly s = tupel_summe(Q, lt_tupel(F));
return(s);
};


// Funktion M (Diagramm)
proc M (matrix Q, ideal F)
{
int m = size(Q);
matrix Qneu[1][m];
for (int k=1; k<=m; k=k+1)
{
	if (Q[1,k] != 0 and leadexp(Q[1,k]*F[k]) == max_mg(Q,F))
	{
		Qneu[1,k] = lead(Q[1,k]);
	} else
	{
		Qneu[1,k] = 0;
	}
}
return(Qneu);
};


// Bildet Tupel Q_ij, mit dem sich lambda(S_ij)=sum(q_i*f_i) speziell erzeugen laesst
proc qij (int i, int j, ideal F)
{
int s = size(F);
matrix q[1][s];
poly g = lambda(sij_poly(i,j,F),F);
list dv = division(g, F);
for (int k=1; k<=s; k=k+1)
{
	q[1,k] = dv[1][k,1];
}
return(q);
};


// Bildung von Urbildern der S_ij bez. M
proc sij_urbild (int i, int j, ideal F)
{
printf("-->S_%s%s-Urbild------------------------------------start", i, j);
print("");
matrix S = sij_poly(i,j,F);
matrix Q = qij(i,j,F);
matrix TF = ideal_tupel(F);
matrix t;
poly p = lambda(S,F);

	printf("S_%s%s  = (%s)", i, j, S);
	print("");
	printf(" *( F = (%s) )", TF);
	print("");
	printf("      = sum(s_i*f_i) = lambda(S_%s%s) = %s", i, j, p);
	print("");
	print("");
	printf("Q_%s%s  = (%s)", i, j, Q);
	printf("");
	printf(" *( F = (%s) )", TF);
	print("");
	printf("      = sum(q_i*f_i) = %s", p);
	printf("");
	print("");

matrix SQ = tupel_diff(S, Q);

	printf("S_%s%s - Q%s%s = (%s)", i, j, i, j, SQ);
	print("");

p = lambda(SQ,F);

	printf("lambda(S_%s%s-Q%s%s) = %s", i, j, i, j, p);
	print("");
	printf("S_%s%s = (%s)", i, j, S);
	print("");

t = M(SQ, F);

	printf("M(S_%s%s-Q%s%s) = (%s)", i, j, i, j, t);
	printf("");

if (p == 0 and t == S)
{
	print("");
	printf("=> lambda(S_%s%s-Q%s%s) liegt in syz(F) und S_%s%s-Q%s%s ist Urbild von S_%s%s bez. M",i,j,i,j,i,j,i,j,i,j);
	print("");
}
printf("-->S_%s%s-Urbild--------------------------------------end", i, j);
print("");
};


// Multipliziert ein Polynom mit jedem Element eines Tupels
proc poly_mal_tupel (poly p, matrix T)
{
int m = size(T);
for (int k=1; k<=m; k=k+1)
{
	T[1,k] = p*T[1,k];
}
return(T);
};


// Bildung von Urbildern der hom. Syz. von LT(F) bez. M
proc s_urbild (matrix S, ideal F)
{
printf("-->S-Urbild--------------------------------------------------start");
printf("-->von S = (%s)", S);
print("");
int m = size(S);
int l = 0;
matrix q;
matrix t;
list SQ;

l=f_hom(S,F);
print("");
syzg(S,lt_tupel(F));
print("");
s_zerlegung(S,F);
print("");
print("");

for (int k=1; k<=m-1; k=k+1)
{
	SQ[k] = tupel_diff(sij_poly(k,m,F),qij(k,m,F));
	sij_urbild(k,m,F);
	print("");
} 
print("");
for (k=1; k<=m-1; k=k+1)
{
	SQ[k] = poly_mal_tupel(ci_X(S,k,m,F),SQ[k]);
	printf("-> Multipliziere (S_%s%s - Q_%s%s) mit c_%s*X^(u-v_%s%s) = q_%s%s", k, m, k, m, k, k, m, k, m);
	print("");
}
for (k=2; k<=m-1; k=k+1)
{
	SQ[1] = tupel_sum(SQ[1],SQ[k]);
}
print("");
printf("sum(qij*(S_ij-Q_ij) = (%s)", SQ[1]);
print("");
poly p = lambda(SQ[1], F);
printf("-> in lambda: Summe aller qij*(S_ij-Q_ij)*fi = %s", p);
print("");
printf("sum(qij*S_ij) = S = (%s)", S);
print("");
q = M(SQ[1], F);
printf("M(sum(qij*(S_ij-Q_ij))) = (%s)", q);
print("");
if (p == 0 and q == S) {
	print("");
	printf("=> sum(qij*(S_ij-Q_ij)) liegt in syz(F) und ist Urbild von S bez. M");
	print("");
}
printf("-->S-Urbild----------------------------------------------------end");
print("");
};


//Test mit n zufälligen homogenen Syzygien von LT(F)
proc s_test (int n, ideal F)
{
printf("-->Test mit %s homogenen Syzygien---------------------------------------start", n);
print("");
matrix hom_s;
for (int k=1; k<=n; k=k+1)
{
	hom_s = hom(F);
	s_urbild(hom_s, F);
	printf("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
} 
printf("-->Test mit %s homogenen Syzygien-----------------------------------------end", n);
print("");
};

LIB "poly.lib";

ring r=0,(x,y,z),lp;
ideal F= x2y+y, x3z+2, y6z2 + z2;
F=groebner(F);
matrix S = hom(F);