// This magma script implements the computation of the discriminant as described in Section 2.1 of the paper
// "A database of nonhyperelliptic genus 3 curves over Q", by Andrew V. Sutherland

// We order monomials in the variables x,y,z lexicogrpahically, which means exponent vectors are in reverse lex order
// x*x*x*x is the lexicographically least monomial, so <4,0,0> is the first exponent vector

E1:=[<0,0,1>,<0,1,0>,<1,0,0>];
E2:=[<0,0,2>,<0,1,1>,<0,2,0>,<1,0,1>,<1,1,0>,<2,0,0>];
E4:=[<4,0,0>,<3,1,0>,<3,0,1>,<2,2,0>,<2,1,1>,<2,0,2>,<1,3,0>,<1,2,1>,<1,1,2>,<1,0,3>,<0,4,0>,<0,3,1>,<0,2,2>,<0,1,3>,<0,0,4>];
E4tab:=AssociativeArray(E4);
for i:=1 to 15 do E4tab[E4[i]]:=i; end for;

// 15 coefficients of generic ternary quartic form are m1,...,m15 (indexed in lex monomial order)
coeff_ring<m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15>:=PolynomialRing(Integers(),15);
tqf_ring<x,y,z>:=PolynomialRing(coeff_ring,3);

// functions to convert exponent vectors to monomials and back (from Ed to Bd and back)
monomial:=function(e) return x^e[1]*y^e[2]*z^e[3]; end function;
exponent:=function(m) return <Degree(m,i):i in [1..3]>; end function;

// Returns the 15 coefficients of given ternary quartic form in lex monomial order
function tqf_coeffs(f)
    c:=[BaseRing(Parent(f))!0:i in [1..15]];
    for t in Terms(f) do c[E4tab[exponent(t)]]:=Coefficients(t)[1]; end for;
    return c;
end function;
    
// Given a vector of 15 coefficients in lex monomial order and a polynomial ring R in (at least) 3 variables,
// constructs the corresponding ternary quartic form in the variables R.1, R.2, R.3
function tqf_from_coeffs(c,R)
    return &+[c[i]*Evaluate(monomial(E4[i]),[R.1,R.2,R.3]):i in [1..15]];
end function;

// the three derivatives of a generic ternary quartic form
tqf_generic:=tqf_from_coeffs([coeff_ring.i:i in [1..15]],tqf_ring);
tqf_derivatives:=[Derivative(tqf_generic,tqf_ring.i):i in [1..3]];

// Evaluate the linear operator D:C[x]_2* --> C[x]_4 defined in Section 2.1 as the determinant of a 3x3 matrix (d=3)
// on a basis element identified by an exponent vector v in E2
function tqf_D(v)
    function Fi (i,v)
        r:=[tqf_ring!0,0,0];
        for t in Terms(tqf_derivatives[i]) do
            e:=exponent(t);
            if e[1] gt v[1] then j:=1; else j:= e[2] gt v[2] select 2 else 3; end if;
            r[j] +:= Coefficients(t)[1] * monomial(<e[k] - (k eq j select v[j]+1 else 0):k in [1..3]>);
        end for;
        return r;
    end function;
    return Determinant(Matrix([Fi(i,v):i in [1..3]]));
end function;

// Evaluates the linear operator T:C[x]_1^3 --> C[x]_4 defined in Section 2.1 as the product
// of a derivative identified by i with a monomial identified by and exponent vector w
function tqf_T(i,w)
    return tqf_coeffs(tqf_derivatives[i]*monomial(w));
end function;

// Build 15x15 matrix of the linear operator Phi as described in Section 2.1 of the paper
// The first 6 rows are values of tqf_D(v) with v identifies a monomial of degree 2
// The last 9 rows are values of tqf_T(i,j) where i identifies a derivative and j identifies a monomial of degree 1
Phi_matrix:=[tqf_coeffs(tqf_D(v)):v in E2];
for w in E1 do for i in [1..3] do Append(~Phi_matrix,tqf_T(i,w)); end for; end for;
Phi_matrix:=Matrix(Phi_matrix);
 
// Computes the discriminant of a ternary quartic form f by computing the determinant of the 15 x 15 matrix
// obtained by instantiating the discriminant matrix with the coefficients of f
function tqf_disc(f)
    // compute the discriminant Delta_4(f) via Propositions 2.5 and 2.7 as -4^7 * det Phi_matrix
    // The determinant is always divisible by 2^14, we remove this factor when computing the discriminant,
    // This ensures the content of the disriminant is 1, as required by Definition 2.1.
    return -ExactQuotient(Determinant(Evaluate(Phi_matrix,tqf_coeffs(f))),2^14);
end function;
