### This program calculates the classical invariants of a binary quartic.
### More precisely, suppose we are given a homogeneous quartic polynomial
### in two variables f = f(x,y) and we ask for polynomials in the
### coefficients of f(x,y) that are invariant under SL(2,C) acting as
### a change of variables in the usual way.  These are the "classical
### Invariants" of f(x,y) and they form a ring under addition and
### multiplication in the obvious way.  It is a classical theorem that
### this ring is freely generated by two generators J and K, homogeneous
### of degrees 2 and 3 respectively.  It is J and K that are calculated by
### this program.  It then goes on to calculate j = 4 J^3/(4 J^3 - 24 K^2).
### The reason for this is that j is then easily seen to be invariant under
### GL(2,C) rather than just SL(2,C).  Equivalently, it is invariant under
### scaling of f(x,y) by any non-zero constant.  It turns out to be the
### classical j-invariant of the elliptic curve in the complex projective
### plane defined by the equation y^2 z^2 = f(x,y) (assuming, of course, that
### the curve is non-singular (this being equivalent to the denominator
### of j being non-zero)).  To see this it suffices to run the program,
### using for f(x,y) a non-singular quartic in one of the standard normal
### forms (e.g. Weierstrass) and check that the answer is the classical one
### for the chosen normal form.  [WARNING: There are two common normalisations
### for the j-invariant.  The one here follows K. Saito "Einfach-elliptische
### Singularitaeten" Invent. Math. 23 (1974) 289--325.  The other one is
### 1728 times this one.]  For convenience, several possible normal
### forms are included in the program.  Finally, the program calculates
### k = J^3/(6 K^2).  This is also an absolute invariant under
### the action of GL(3,C) and the one needed in the article "Continuous
### Invariants of Isolated Hypersurface Singularities".

### Here are several possible quartic polynomials, including some standard
### normal forms.  Uncomment one of these or type in your own.
#f:=x*(x-y)*(x-lambda*y)*y;
#f:=(x^3-27*c4*x*y^2-54*c6*y^3)*y;
#f:=(x^3+p*x*y^2+q*y^3)*y;
#f:=(4*x^3-g2*x*y^2-g3*y^3)*y;
#f:=x^4+t*x^2*y^2+y^4;
#f:=y^4+(e*y^2+x^2)^2;
#f:=A*x^4+4*B*x^3*y+6*C*x^2*y^2+4*DD*x*y^3+EE*y^4;
f:=-t/2*x^4+6*x^2*y^2-t/2*y^4;

### Here is an optional general change of coordinates.
### You can use this to check invariance.
#f:=subs({x=a*x+b*y,y=c*x+d*y},f);

### The program itself starts here.
with(linalg):
a1111:=diff(f,x$4)/24:
a1112:=diff(diff(f,x$3),y)/24:
a1121:=diff(diff(f,x$3),y)/24:
a1211:=diff(diff(f,x$3),y)/24:
a2111:=diff(diff(f,x$3),y)/24:
a1122:=diff(diff(f,x$2),y$2)/24:
a1212:=diff(diff(f,x$2),y$2)/24:
a1221:=diff(diff(f,x$2),y$2)/24:
a2112:=diff(diff(f,x$2),y$2)/24:
a2121:=diff(diff(f,x$2),y$2)/24:
a2211:=diff(diff(f,x$2),y$2)/24:
a1222:=diff(diff(f,x),y$3)/24:
a2122:=diff(diff(f,x),y$3)/24:
a2212:=diff(diff(f,x),y$3)/24:
a2221:=diff(diff(f,x),y$3)/24:
a2222:=diff(f,y$4)/24:
epsilon11:=0:
epsilon22:=0:
epsilon12:=1:
epsilon21:=-1:

for i to 2 do:
 for j to 2 do:
  for k to 2 do:
   for l to 2 do:
    b.i.j.k.l:=simplify(sum('sum(
               'a.i.j.p.q*epsilon.p.k*epsilon.q.l',
               'p'=1..2)',
               'q'=1..2)):
   od:
  od:
 od:
od:

BB:=array([[b1111,b1112,b1121,b1122],
           [b1211,b1212,b1221,b1222],
           [b2111,b2112,b2121,b2122],
           [b2211,b2212,b2221,b2222]]):
J:=factor(trace(BB^2));
K:=factor(trace(BB^3));
j:=factor(4*J^3/(4*J^3-24*K^2));
k:=factor(J^3/(6*K^2));