Code

ToCanonical:=proc(f) global i;
local Qpp,Qpq,Qqq,Qppp,Qppq,Qpqq,Qqqq,Qpppp,Qpppq,Qppqq,Qpqqq,Qqqqq,d,F,Fp, Fq, Fpp,
Fpq, Fqq, Fppp,Fppq, Fpqq,Fqqq,Fpppp,Fpppq,Fppqq,Fpqqq,Fqqqq,n,N1,N2,N3,
A,inv,p,q,dA,I1,I2,I3,mu,result,i4,i5;
with(linalg):
F:=f(p,q);
Fp:=diff(f(p,q),p);
Fq:=diff(f(p,q),q);
Fpp:=diff(Fp,p);
Fpq:=diff(Fp,q);
Fqq:=diff(Fq,q);
Fppp:=diff(Fpp,p);
Fppq:=diff(Fpq,p);
Fqqq:=diff(Fqq,q);
Fpqq:=diff(Fpq,q);
Fpppp:=diff(Fppp,p);
Fpppq:=diff(Fppq,p);
Fppqq:=diff(Fpqq,p);
Fpqqq:=diff(Fpqq,q);
Fqqqq:=diff(Fqqq,q);
n := 3;
N1:=(n-1)/n:
N2:=(n-1)*(n-2)/n^2:
N3:=(n-2)/n:

Qpp:=simplify(Fpp*F-N1*Fp^2);
Qpq:=simplify(Fpq*F-N1*Fp*Fq);
Qqq:=simplify(Fqq*F-N1*Fq^2);
Qppp:=simplify(Fppp*F^2-3*N3*Fp*Fpp*F+2*N2*Fp^3);
Qppq:=simplify(Fppq*F^2-N3*(Fpp*Fq+2*Fpq*Fp)*F+2*N2*Fp^2*Fq);
Qpqq:=simplify(Fpqq*F^2-N3*(Fqq*Fp+2*Fpq*Fq)*F+2*N2*Fp*Fq^2);
Qqqq:=simplify(Fqqq*F^2-3*N3*Fq*Fqq*F+2*N2*Fq^3);
d:=factor(expand(Qpq^2-Qpp*Qqq,power));
if evalb(d=0) then RETURN(`equivalent either to a binary form (case (9) or (10) or to p^3 (case (11))`)
else
I1:=factor(expand(Qppp*Qpqq*Qqq-Qppq^2*Qqq-Qppp*Qqqq*Qpq+Qppq*Qpqq*Qpq+Qppq*Qqqq*Qpp
    -Qpqq^2*Qpp,power))/factor(expand(d^2,power));
I2:=factor(expand(5*Qqqq^2*Qpp^3+5*Qppp^2*Qqq^3+36*Qppq^2*Qpq^2*Qqq+9*Qppq^2*Qqq^2*Qpp
    -4*Qppp*Qpq^3*Qqqq-30*Qppp*Qqq^2*Qppq*Qpq+24*Qppp*Qpq^2*Qqq*Qpqq+6*Qppp*Qqq^2*Qpqq*Qpp
    -6*Qppp*Qpp*Qqq*Qqqq*Qpq-54*Qppq*Qpp*Qqq*Qpqq*Qpq-36*Qppq*Qpq^3*Qpqq+36*Qpqq^2*Qpp*Qpq^2
     +9*Qpqq^2*Qpp^2*Qqq+24*Qppq*Qpp*Qpq^2*Qqqq+6*Qppq*Qpp^2*Qqq*Qqqq-30*Qpqq*Qpp^2*Qqqq*Qpq,
     power))/factor(expand(d^3,power));
I3:=factor(expand(Qppp^2*Qqqq^2-3*Qppq^2*
     Qpqq^2-6*Qppp*Qppq*Qpqq*Qqqq+4*Qppp*Qpqq^3+4*Qqqq*Qppq^3,power))/factor(expand(d^3,power));
i[1]:=simplify(10*(6*I1+1));
i[2]:=simplify(6* I2+126* I1+45* I3-10);
i[3]:=10*(9* I3+2);
if evalb(i[2] <>0) then A:=subs({p=1,q=1},simplify(normal(-27/40*i[1]^3/i[2]^2)))
else
 A:=`not defined`
fi;
if evalb(i[1]=90) and evalb(i[2]=270) and evalb(i[3]=180)  then RETURN([ p*q, ` case (8) `])
else
 if evalb(i[1]=0) and evalb(i[2]=0) and evalb(i[3]=0) then RETURN([p^2+q, ` case (6) `] )
 else
  if evalb(i[1]=0) and evalb(i[2]=0) then RETURN([         p^3-q^2, ` case (4) `])
  else
   i4 := simplify(normal(10 * (i[2])^2 - (i[1])^3));
   i5 := simplify(normal(i[1] * i[3]  - i[2] * (i[1] -30)));
   if evalb(evalb(i4 = 0) and evalb(i5 = 0)) then RETURN([ p^2+q^2+1, ` case (7) `])
   else
    if evalb(i[1]=0) then RETURN([ p^3+1-q^2, ` case (3) `])
    else
     if evalb(i[2]=0) then RETURN([ p^3+p-q^2, ` case (2) `])
     else
      if evalb(i4 = 0) then RETURN([p^2*(p+1)-q^2, ` case (5)         `])
      else
      RETURN ([ p^3+mu*p+1-q^2, `case (1) `, mu=root(A,3)]);
      fi
     fi
    fi
   fi
  fi
 fi
fi
fi;
end: