////////////////////////////////////////////////////////////////////////
// Verbose mode                                                       //
////////////////////////////////////////////////////////////////////////

declare verbose RelationsJacG2, 2;

////////////////////////////////////////////////////////////////////////
//             Generic code for evaluation of polynomials             //
////////////////////////////////////////////////////////////////////////

intrinsic '@'(S::[RngElt],f::RngMPolElt) -> RngElt
    {}
    require #S eq Rank(Parent(f)) :
	"Argument 1 must have length equal to " *
	"the rank of the parent of argument 2.";
    return Evaluate(f,S);
end intrinsic;

intrinsic '@'(S::[RngElt],f::FldFunRatElt) -> RngElt
    {}
    require #S eq Rank(Parent(f)) :
	"Argument 1 must have length equal to " *
	"the rank of the parent of argument 2.";
    return Evaluate(f,S);
end intrinsic;

// Hack to allow evaluation of rational functions.
intrinsic Evaluate(f::FldFunRatMElt, S::[RngElt]) -> RngElt
    {}
    return Evaluate(Numerator(f),S) / Evaluate(Denominator(f),S);;
end intrinsic;

////////////////////////////////////////////////////////////////////////
//  Code specific to solving for relations in the genus 2 Jacobian    //
////////////////////////////////////////////////////////////////////////

intrinsic JacobianWeights() -> SeqEnum
    {}
    wts1 := [ 4, 3, 3, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0];
    wts2 := [-4,-2,-3, 0,-1,-2, 2, 1, 0,-1, 4, 3, 2, 2, 1, 0];
    M := RSpace(Integers(),2);
    return [ M![wts1[i],wts2[i]] : i in [1..16] ],
	[ M![2,-i] : i in [0..6] ] cat [ M!0 ];
end intrinsic;

intrinsic MonomialsOfWeight(W::[RngIntElt] : Prune := true) -> SeqEnum
    {}
    require #W eq 2 : "Argument must be a sequence of two integers.";
    wts, cfs := JacobianWeights();
    w := Universe(wts)!W;
    wghts := [ [i,j,m,n] : i, j in [1..8], m, n in [1..16]
	| (cfs[i]+cfs[j]+wts[m]+wts[n]) eq w
	and i le j and m le n ];
    a, r := JacobianRelations(Rationals());
    F := Universe(a);
    K<f0,f1,f2,f3,f4,f5,f6> := BaseRing(F);
    ff := [f0,f1,f2,f3,f4,f5,f6,1];
    P := Universe(r);
    if not Prune then
	return
	    [ ff[w[1]]*ff[w[2]]*a[w[3]]*a[w[4]] : w in wghts ],
	    [ ff[w[1]]*ff[w[2]]*P.w[3]*P.w[4] : w in wghts ];
    end if;
    ainvs := [ ];
    terms := [ ];
    mons := {@ @};
    for w in wghts do
	x, y, i, j := Explode(w);
	m := P.i*P.j;
	if m notin mons then
	    Append(~ainvs,ff[x]*ff[y]*a[i]*a[j]);
	    Append(~terms,ff[x]*ff[y]*m);
	    Include(~mons,m);
	end if;
    end for;
    return ainvs, terms;
end intrinsic;

intrinsic MonomialWeight(m::RngMPolElt) -> SeqEnum
    {}
    I := [ ];
    E := Exponents(m);
    for i in [1..16] do
	if E[i] ge 1 then Append(~I,i); end if;
	if E[i] ge 2 then Append(~I,i); end if;
    end for;
    i, j := Explode(I);
    wts, cfs := JacobianWeights();
    wi1 := wts[i][1]; wj1 := wts[j][1]; 
    wi2 := wts[i][2]; wj2 := wts[j][2];
    c := Coefficients(m); assert #c eq 1;
    E := Exponents(Monomials(Numerator(c[1]))[1]);
    F := Exponents(Monomials(Denominator(c[1]))[1]);
    c1 := 2*&+E - 2*&+F;
    c2 := &+[ -(i-1)*(E[i]-F[i]) : i in [1..7] ];
    vprintf RelationsJacG2, 2 :
	"W = [%2o,%2o] + [%2o,%2o] + [%2o,%2o]\n",
	wts[i][1], wts[i][2], wts[j][1], wts[j][2], c1, c2;
    W := [ wi1+wj1+c1, wi2+wj2+c2 ];
    return W;
end intrinsic;

intrinsic MonomialRelationsOfWeight(W::[RngIntElt] :
    Relations := [], StrictWeight := true) -> SetIndx
    {}
    require #W eq 2 : "Argument must have length 2.";
    ainvs, mons := MonomialsOfWeight(W);
    V := LinearRelations(ainvs : Relations := Relations);
    WghtRels := {@ @};
    wts, cfs := JacobianWeights();
    for v in Basis(V) do
	f := &+[ v[i]*mons[i] : i in [1..#mons] ];
	k := Random([ i : i in [1..Degree(V)] | v[i] ne 0 ]);
	wghtf := MonomialWeight(v[k]*mons[k]);
	if not StrictWeight or W eq wghtf then
	    Include(~WghtRels,f);
	else
	    vprintf RelationsJacG2 : 
		"Found relation of weight [%2o,%2o] (not [%2o,%2o]):\n",
		wghtf[1], wghtf[2], W[1], W[2];
	    vprint RelationsJacG2 : f;
	end if;
    end for;
    return WghtRels;
end intrinsic;

/*
F<f0,f1,f2,f3,f4,f5,f6> := FunctionField(Rationals(),7 : Global);
P<X00,X01,X02,X03,X04,X05,X06,X07,
  X08,X09,X10,X11,X12,X13,X14,X15> := PolynomialRing(F,16 : Global);
a, r := JacobianRelations([f0,f1,f2,f3,f4,f5,f6]);
R<x1,x2,y1,y2> := PolynomialRing(F,4 : Global);
F1 := y1^2-(f6*x1^6 + f5*x1^5 + f4*x1^4 + f3*x1^3 + f2*x1^2 + f1*x1 + f0);
F2 := y2^2-(f6*x2^6 + f5*x2^5 + f4*x2^4 + f3*x2^3 + f2*x2^2 + f1*x2 + f0);
DomRels := [F1,F2];
I := ideal< R | DomRels >;
//time MonomialRelationsOfWeight([4,0] : Relations := DomRels);
//
TotRels := {@ @};
for i in [6..0 by -1] do
    if i eq 6 then jInt := [-4..-6 by -1];
    elif i eq 5 then jInt := [-2..-5 by -1];
    elif i eq 4 then jInt := [ 0..-4 by -1];
    elif i eq 3 then jInt := [ 2..-3 by -1];
    elif i eq 2 then jInt := [ 4..-2 by -1];
    elif i eq 1 then jInt := [ 3.. 1 by -1];
    elif i eq 0 then jInt := [ 6.. 2 by -1];
    end if;
    for j in jInt do
        printf "Weight: [%2o,%2o]\n", i, j;
        NewRels := MonomialRelationsOfWeight([i,j] : Relations := DomRels);
        print NewRels;
        TotRels join:= {@ Factorization(f)[1][1] : f in NewRels @};
        print "Total relations:", #TotRels;
    end for;
end for;
*/

////////////////////////////////////////////////////////////////////////
// Completely generic...
////////////////////////////////////////////////////////////////////////

intrinsic AlgebraicRelations(S::[FldFunRatElt] :
    DomainRelations := [],
    CodomainRelations := []) -> RngMPol
    {}
    m := #S;
    R := Universe(S);
    n := Rank(R);
    P := PolynomialRing(R,m+n);
    AssignNames(~P,
	[ "r" cat IntegerToString(i): i in [1..m]] cat 
	[ "x" cat IntegerToString(i): i in [1..n]]);
    vars := [ P.j : j in [m+1..m+n] ];
    rels := [
	Denominator(S[i])(vars)*P.i -
	Numerator(S[i])(vars) : i in [1..m] ];
    if #DomainRelations ne 0 then
	Q := Universe(DomainRelations);
	require Type(Q) eq RngMPol : 
	    "Parameter DomainRelations must consist of polynomials.";
	require Rank(Q) eq n : 
	    "Rank of universe of DomainRelations " *
	    "must equal universe of argument 1.";
	rels cat:= [ f([ P.j : j in [m+1..m+n] ]) : f in Relations ];
    end if;
    return EliminationIdeal(ideal< P | rels >,{P.i : i in [1..m]});
end intrinsic;

intrinsic LinearRelations(S::[FldFunRatElt] : Relations := [])
    -> ModTupFld
    {}
    for i in [1..#S] do
	den := Denominator(S[i]); 
	if den ne 1 then
	    S := [ den*f : f in S ];
	end if;
    end for;
    S := [ Numerator(f) : f in S ]; 
    if #Relations ne 0 then
	vprint RelationsJacG2 : "Putting functions in normal form.";
	P := Universe(Relations);
	require P eq Universe(S) : "Invalid Relations parameter.";
	I := ideal< Universe(S) | Relations >;
	S := [ NormalForm(f,I) : f in S ];
    end if;
    vprint RelationsJacG2 : "Extracting coefficients for matrix.";
    K := BaseRing(Universe(S));
    mons := &join[ SequenceToSet(Monomials(f)) : f in S ];
    vprint RelationsJacG2 : "Number of monomials:", #mons;
    vprintf RelationsJacG2 : "Matrix size %o x %o\n", #S, #mons;
    M := Matrix([ [ MonomialCoefficient(f,m) : m in mons ] : f in S ]);
    if GetVerbose("RelationsJacG2") ge 2 then
	for j in [1..#mons] do
	    "j =", j;
	    for i in [1..#S] do
		if M[i,j] ne 0 then
		    printf "M[%3o,%2o] = %o\n", i, j, M[i,j];
		end if;
	    end for;
	end for;
    end if;
    return Kernel(M);
end intrinsic;

