##############################################################################
# Bertini involution by Alex Degtyarev
#
# This is a Maple implementation of E. Moody's formulas for the Bertini involution
# and my extension for the associated map P^2 --> Sigma_2 (mainly, the equation
# of the ramification locus and the relation between the two sections).
#
# - Ethel I. Moody, Notes on the Bertini involution, Bull. Amer. Math. Soc. 49
#   (1943), 433-436
# - Alex Degtyarev, The Bertini involution, unpublished
#   http://www.fen.bilkent.edu.tr/~degt/papers/papers.htm
##############################################################################

##############################################################################
# Call init_cubics(f1, f2); to initialize the coefficients with two particular
# cubic equations f1, f2 in x[1], x[2], x[3]
##############################################################################

## These are to be used as variables:
x := 'x':  y := 'y':
## And these are the coefficients
a1 := 'a1':  b1 := 'b1':  c1 := 'c1':
a2 := 'a2':  b2 := 'b2':  c2 := 'c2':   # the primed versions

## Service functions
# coefficient in f of x[1]^i*x[2]^j*x[3]^k
cf := proc(f, x, i, j ,k)
 return coeff(coeff(coeff(f, x[1], i), x[2], j), x[3], k);
end:
# same simplified
scf := proc(f, x, i, j, k) return simplify(cf(f, x, i, j, k)); end:
# same factored
fcf := proc(f, x, i, j, k) return factor(cf(f, x, i, j, k)); end:

# all coefficients in a homogeneous degree n polynomial f in x
cfs := proc(f, x, n)
 return seq(seq(cf(f, x, i, j, n-i-j), j=0..n-i), i=0..n);
end:

# swapping primed and not primed letters
swap := proc(expr)
 local tempa, tempb, tempc, temp;
 temp := subs(a1=tempa, b1=tempb, c1=tempc, expr);
 temp := subs(a2=a1, b2=b1, c2=c1, temp);
 return subs(tempa=a2, tempb=b2, tempc=c2, temp);
end:
# {e} = e + e'
sswap := e -> e + swap(e):

# the brackets [.]_? as in \cite{Moody}
dv := proc(f, d) return simplify(f/d); end:

# to be given: two cubics f1, f2 in x[1], x[2], x[3]
# initialize w1, w2 and coefficients a_i, b_i, c_i
init_cubics := proc(f1, f2)
 global a1, b1, c1, a2, b2, c2, w1, w2;
 w1 := f1; w2 := f2;
 a1[1] := cf(w1, x, 1, 0, 2);  a2[1] := cf(w2, x, 1, 0, 2);
 a1[2] := cf(w1, x, 0, 1, 2);  a2[2] := cf(w2, x, 0, 1, 2);
 b1[1] := cf(w1, x, 2, 0, 1);  b2[1] := cf(w2, x, 2, 0, 1);
 b1[2] := cf(w1, x, 1, 1, 1);  b2[2] := cf(w2, x, 1, 1, 1);
 b1[3] := cf(w1, x, 0, 2, 1);  b2[3] := cf(w2, x, 0, 2, 1);
 c1[1] := cf(w1, x, 2, 1, 0);  c2[1] := cf(w2, x, 2, 1, 0);
 c1[2] := cf(w1, x, 1, 2, 0);  c2[2] := cf(w2, x, 1, 2, 0);
 return ;
end:

# original expressions for the two cubics
 w1 := x[3]^2*(a1[1]*x[1] + a1[2]*x[2]) +
       x[3]*(b1[1]*x[1]^2 + b1[2]*x[1]*x[2] + b1[3]*x[2]^2) +
       (c1[1]*x[1]^2*x[2] + c1[2]*x[1]*x[2]^2):
 w2 := x[3]^2*(a2[1]*x[1] + a2[2]*x[2]) +
       x[3]*(b2[1]*x[1]^2 + b2[2]*x[1]*x[2] + b2[3]*x[2]^2) +
       (c2[1]*x[1]^2*x[2] + c2[2]*x[1]*x[2]^2):
       
# just in case:
p := 'p': q := 'q': r := 'r': s := s:
## These polinomials are defined as functions for easy substitution:
P := proc(t1, t2) return p[2]*t1^2 + p[1]*t1*t2 + p[0]*t2^2; end:
Q := proc(t1, t2)
 return q[4]*t1^4 + q[3]*t1^3*t2 + q[2]*t1^2*t2^2 + q[1]*t1*t2^3 + q[0]*t2^4;
end:
R := proc(t1, t2)
 return r[3]*t1^3 + r[2]*t1^2*t2 + r[1]*t1*t2^2 + r[0]*t2^3;
end:
S := proc(t1, t2) return s[2]*t1^2 + s[1]*t1*t2 + s[0]*t2^2; end:
## The homogeneous relation for K^2:
KK := proc(y, t1, t2)
 return - 4*y^3 + y^2*P(t1, t2) + y*Q(t1, t2) + R(t1, t2)^2;
end:
## The affine equation of the ramification locus
KS := proc(x, y) return KK(y, x, 1); end:

 for i from 1 to 2 do A[i] := simplify(subs(x=y, a1[i]*w2 - a2[i]*w1)); od:
 for i from 1 to 3 do B[i] := simplify(subs(x=y, b1[i]*w2 - b2[i]*w1)); od:
 for i from 1 to 2 do C[i] := simplify(subs(x=y, c1[i]*w2 - c2[i]*w1)); od:
 i := 'i':
 
 gamma4 := y[1]*A[1] + y[2]*A[2]:
 kappa := a1[1]*b2[1] - a2[1]*b1[1]:
 C5 := A[2]*dv(B[1] + kappa*y[1]*y[3]^2, y[2]) +
  dv(A[1] - kappa*y[1]^2*y[3], y[2])*dv(A[2]*y[3]+B[3]*y[2], y[1]) +
  kappa*B[3]*y[1]*y[3]:
 phi6 := simplify(A[1]*C[2] + y[3]*C5):
 psi6 := simplify(A[2]*C[1] + y[3]*C5):
# these are $r'_1$, $r'_3$ in \cite{Moody} (there is no r'_2)
 rr[1] := B[1]*A[2]^2 - B[2]*A[1]*A[2] + B[3]*A[1]^2:
 rr[3] := A[2]*C[1] - A[1]*C[2]:
# and these are $r_1$, $r_2$, $r_3$
 Rr[1] := A[2]*rr[1]:
 Rr[2] := -A[1]*rr[1]:
 Rr[3] := A[1]*A[2]*rr[3]:
# the Bertini involution
 yb[1] := phi6*dv(A[2]^2*phi6 + B[3]*rr[1], y[1]):
 yb[2] := psi6*dv(A[1]^2*psi6 + B[1]*rr[1], y[2]):
 yb[3] := psi6*phi6*C5:
# the coefficients in psi6 = phi6 + s[2]*w1^2 + s[1]*w1*w2 + s[0]*w2^2:
 s[0] := a1[2]*c1[1] - a1[1]*c1[2]:
 s[1] := (a1[1]*c2[2]+a2[1]*c1[2])-(a1[2]*c2[1]+a2[2]*c1[1]):
 s[2] := swap(s[0]):   # will not work if a1, ..., c2 are evaluated
# the fixed point set
 K := psi6*dv(A[1]*y[3] + B[1]*y[1], y[2]) -  phi6*dv(A[2]*y[3] + B[3]*y[2], y[1]):
# the coefficients in the equation for K^2:
 r[0] := -a1[1]*b1[2]*c1[2]+a1[1]*b1[3]*c1[1]+a1[2]*b1[1]*c1[2]:
 r[1] :=
   (a1[1]*b1[2]*c2[2]+a1[1]*b2[2]*c1[2]+a2[1]*b1[2]*c1[2])
  -(a1[1]*b1[3]*c2[1]+a1[1]*b2[3]*c1[1]+a2[1]*b1[3]*c1[1])
  -(a1[2]*b1[1]*c2[2]+a1[2]*b2[1]*c1[2]+a2[2]*b1[1]*c1[2]):
 r[2] := -swap(r[1]):   # will not work if a1, ..., c2 are evaluated
 r[3] := -swap(r[0]):   # will not work if a1, ..., c2 are evaluated

 p[0] := b1[2]^2-4*a1[2]*c1[1]-4*b1[1]*b1[3]+8*a1[1]*c1[2]:
 p[1] := -8*(a1[1]*c2[2]+a2[1]*c1[2])+4*(a1[2]*c2[1]+a2[2]*c1[1])
  +4*(b1[1]*b2[3]+b2[1]*b1[3])-(b1[2]*b2[2]+b2[2]*b1[2]):
 p[2] := swap(p[0]):   # will not work if a1, ..., c2 are evaluated
 
 q[0] := 4*(a1[1]*c1[2] - b1[1]*b1[3])*s[0] + 2*b1[2]*r[0]:
 q[1] :=
   4*(a1[1]*a1[1]*c1[2]*c2[2]+a1[1]*a1[1]*c2[2]*c1[2]
      +a1[1]*a2[1]*c1[2]*c1[2]+a2[1]*a1[1]*c1[2]*c1[2])
  -4*(a1[1]*a1[2]*c1[1]*c2[2]+a1[1]*a1[2]*c2[1]*c1[2]
      +a1[1]*a2[2]*c1[1]*c1[2]+a2[1]*a1[2]*c1[1]*c1[2])
  -4*(a1[1]*b1[1]*b1[3]*c2[2]+a1[1]*b1[1]*b2[3]*c1[2]
      +a1[1]*b2[1]*b1[3]*c1[2]+a2[1]*b1[1]*b1[3]*c1[2])
  +2*(a1[1]*b1[2]*b1[2]*c2[2]+a1[1]*b1[2]*b2[2]*c1[2]
      +a1[1]*b2[2]*b1[2]*c1[2]+a2[1]*b1[2]*b1[2]*c1[2])
  -2*(a1[2]*b1[1]*b1[2]*c2[2]+a1[2]*b1[1]*b2[2]*c1[2]
      +a1[2]*b2[1]*b1[2]*c1[2]+a2[2]*b1[1]*b1[2]*c1[2])
  -2*(a1[1]*b1[2]*b1[3]*c2[1]+a1[1]*b1[2]*b2[3]*c1[1]
      +a1[1]*b2[2]*b1[3]*c1[1]+a2[1]*b1[2]*b1[3]*c1[1])
  +4*(a1[2]*b1[1]*b1[3]*c2[1]+a1[2]*b1[1]*b2[3]*c1[1]
      +a1[2]*b2[1]*b1[3]*c1[1]+a2[2]*b1[1]*b1[3]*c1[1]):
 q[2] :=
  -4*(a1[1]*a1[1]*c2[2]*c2[2]+a1[1]*a2[1]*c1[2]*c2[2]+a2[1]*a1[1]*c1[2]*c2[2]
      +a1[1]*a2[1]*c2[2]*c1[2]+a2[1]*a1[1]*c2[2]*c1[2]+a2[1]*a2[1]*c1[2]*c1[2])
  +4*(a1[1]*a1[2]*c2[1]*c2[2]+a1[1]*a2[2]*c1[1]*c2[2]+a2[1]*a1[2]*c1[1]*c2[2]
      +a1[1]*a2[2]*c2[1]*c1[2]+a2[1]*a1[2]*c2[1]*c1[2]+a2[1]*a2[2]*c1[1]*c1[2])
  +4*(a1[1]*b1[1]*b2[3]*c2[2]+a1[1]*b2[1]*b1[3]*c2[2]+a2[1]*b1[1]*b1[3]*c2[2]
      +a1[1]*b2[1]*b2[3]*c1[2]+a2[1]*b1[1]*b2[3]*c1[2]+a2[1]*b2[1]*b1[3]*c1[2])
  -2*(a1[1]*b1[2]*b2[2]*c2[2]+a1[1]*b2[2]*b1[2]*c2[2]+a2[1]*b1[2]*b1[2]*c2[2]
      +a1[1]*b2[2]*b2[2]*c1[2]+a2[1]*b1[2]*b2[2]*c1[2]+a2[1]*b2[2]*b1[2]*c1[2])
  +2*(a1[1]*b1[2]*b2[3]*c2[1]+a1[1]*b2[2]*b1[3]*c2[1]+a2[1]*b1[2]*b1[3]*c2[1]
      +a1[1]*b2[2]*b2[3]*c1[1]+a2[1]*b1[2]*b2[3]*c1[1]+a2[1]*b2[2]*b1[3]*c1[1])
  +2*(a1[2]*b1[1]*b2[2]*c2[2]+a1[2]*b2[1]*b1[2]*c2[2]+a2[2]*b1[1]*b1[2]*c2[2]
      +a1[2]*b2[1]*b2[2]*c1[2]+a2[2]*b1[1]*b2[2]*c1[2]+a2[2]*b2[1]*b1[2]*c1[2])
  -4*(a1[2]*b1[1]*b2[3]*c2[1]+a1[2]*b2[1]*b1[3]*c2[1]+a2[2]*b1[1]*b1[3]*c2[1]
      +a1[2]*b2[1]*b2[3]*c1[1]+a2[2]*b1[1]*b2[3]*c1[1]+a2[2]*b2[1]*b1[3]*c1[1]):
 q[3] := swap(q[1]):   # will not work if a1, ..., c2 are evaluated
 q[4] := swap(q[0]):   # will not work if a1, ..., c2 are evaluated
 
## Extra coefficients s[i] corresponding to another BASEPOINT [1, z2, z3]
# Returns in the form table([0=..., 1=..., 2=...]);
#ss[0] := s[0]*z[1]^2+a1[2]*c1[2]*z[1]*z[2]+(a1[2]*b1[2]-a1[1]*b1[3])*z[1]*z[3]
#        +a1[2]*b1[3]*z[2]*z[3]+a1[2]^2*z[3]^2:
#ss[2] := swap(ss[0]):
#ss[1] := s[1]*z[1]^2+(a1[2]*c2[2]+a2[2]*c1[2])*z[1]*z[2]
#        +(-a1[1]*b2[3]+a1[2]*b2[2]-a2[1]*b1[3]+a2[2]*b1[2])*z[1]*z[3]
#        +(a1[2]*b2[3]+a2[2]*b1[3])*z[2]*z[3]+2*a1[2]*a2[2]*z[3]^2:
alt_section := proc(z2, z3)
 return table([
    0 = s[0]+(a1[2]*c1[2]*z2+(a1[2]*b1[2]-a1[1]*b1[3])*z3)
        +a1[2]*b1[3]*z2*z3+a1[2]^2*z3^2,
    2 = s[2]+(a2[2]*c2[2]*z2+(a2[2]*b2[2]-a2[1]*b2[3])*z3)
        +a2[2]*b2[3]*z2*z3+a2[2]^2*z3^2,
    1 = s[1]
        +((-a2[2]*c1[2]-a1[2]*c2[2])*z2
         +(-a2[2]*b1[2]+a2[1]*b1[3]-a1[2]*b2[2]+a1[1]*b2[3])*z3)
        +(-a1[2]*b2[3]-a2[2]*b1[3])*z2*z3-2*a1[2]*a2[2]*z3^2
 ]);
end:
SS := proc(ss, t1, t2) return ss[2]*t1^2 + ss[1]*t1*t2 + ss[0]*t2^2; end:
 
## Just check the identities to make sure that everything has initialized properly
check_data := proc()
 local K0, n, i;
 print("The expansion of $\\psi_6$");
 print(simplify(psi6 - subs(x=y, phi6 + s[2]*w1^2 + s[1]*w1*w2 + s[0]*w2^2)));
 print("The expansion of $K^2$ (coefficientwise)");
 K0 := -K^2 + subs(x=y, KK(phi6, w1, w2));
 for n from 0 to 18 do print('n'=n, seq(scf(K0, y, i, n-i, 18-n), i=0..n)); od;
 print("The identities for $r$, $\\phi$, and $\\psi$");
 print(simplify(y[1]*Rr[3] - y[3]*Rr[1] + A[2]*gamma4*phi6));
 print(simplify(y[2]*Rr[3] - y[3]*Rr[2] - A[1]*gamma4*psi6));
 print("The divisibility of $y_3 z_1 - y_1 z_3$ by $K$");
 print(simplify((y[3]*yb[1]-y[1]*yb[3])/phi6/K + A[2]));
 print("The divisibility of $y_3 z_2 - y_2 z_3$ by $K$");
 print(simplify((y[3]*yb[2]-y[2]*yb[3])/psi6/K - A[1]));
 print("All done");
end:

##############################################################################
# End of file Bertini.txt
##############################################################################

##############################################################################
# Common stuff for computing maximizing sextics
##############################################################################

## Printing routines
ssets := table():
lp := proc() printf("%A %Q\n", args); end:
lp3 := proc(s1, s2, s3)
 printf("%A\n%A\n%A\n", s1, s2, s3);
end:
print1 := proc()
 lp("## --------------------------------------------------------------------");
end:
print2 := proc()
 lp("## ====================================================================");
end:
set_set := proc(s, set)
 global sset, idx, ssets;
 sset := op(2, set);
 idx := op(1, set);
 ssets[idx] := sset;
 lp("#######################################################################");
 lp(s, idx=sset);
end:
print_set := proc()
 print1();
 lp("## The set of singularities:", idx=sset);
end:
set_alias := proc(value)
 print(value);
 alias(value):
end:

## These may change
get_section_plain := mn -> subs(mn, parabola):
get_section := get_section_plain:

get_ab_plain := mn -> subs(mn, ab):
get_ab := get_ab_plain:

## Printing the solution
print_solution := proc()
 global sub_mn, sub_ab, curve, section, poly, found_s, found_d, complexify;
 local res, rad, Ji;
 rad := false;
 print_set();
 sub_mn[idx] := found_mn;
 curve[idx] := found_s;
 if nargs > 0 then
  if args[nargs] = complex then
   complexify := true;
  fi;
  if args[nargs] = radical then
   rad := true;
  fi;
  if type(args[1], `=`) then
   set_alias(args[1]);
   poly[idx] := op(1, op(2, args[1]));
   print('poly'=poly[idx]);
   print('discr'=ifactor(discrim(poly[idx], _Z)), 'galois'=galois(poly[idx]));
  fi;
 fi;
 if complexify then
  sub_mn[idx] := convert(sub_mn[idx], radical);
  found_s := get_sextic(sub_mn[idx]);
  curve[idx] := found_s;
  found_d := discr(curve[idx]);
 end;
 res := sub_mn[idx];
 print([mu,nu]=sub_mn[idx]);
 print('sextic'=curve[idx]);
 print('discr'=found_d);
 print_set();
 if rad then res := convert(res, radical); fi;
 sub_ab[idx] := simplify(get_ab(res));
 Ji := simplify(evala(j_invariant(found_d)));
 if complexify then
  sub_ab[idx] := simplify(evalc(sub_ab[idx]));
  Ji := simplify(evalc(Ji));
 fi;
 if rad then
  sub_ab[idx] := simplify(evala(sub_ab[idx]));
  Ji := simplify(evala(convert(Ji, radical)));
 fi;
 section[idx] := factor(simplify(get_section(res)));
 complexify := false;
 print([mu,nu]=res);
 print('sub_ab'=sub_ab[idx]);
 print('section'=section[idx]);
 print('J'=Ji);
 print2();
end:

## Initialization with two cubics "u1", "u2"
init_sextics := proc(u1, u2)
 local cx, cy;
 global bu1, bu2, bphi, sub;
 init_cubics(u1, u2);
# More convenient coordinates for the final results
 cx := simplify(cf(u1, x, 0, 2, 1)/cf(u1, x, 1, 2, 0));
 cy := simplify(cf(u2, x, 2, 0, 1)/cf(u2, x, 2, 1, 0));
 sub := [x[1]=1-cx*x, x[2]=y-cy*x, x[3]=x];
 print('sub'=sub);
 bu1 := collect(simplify(subs(sub, u1)), [y, x]);
 bu2 := collect(simplify(subs(sub, u2)), [y, x]);
 bphi := collect(simplify(subs(sub, subs(y=x, phi6))), [y, x]):
 print2();
## print('bu1'=bu1);
## print('bu2'=bu2);
end:

##############################################################################
# Printing the results
##############################################################################

## The pull-back of a section "section" (in the form a*t^2 + b*t + c)
#  "ab" is in the form [alpha=..., beta=..., ...]
Result := false:
normalize := true:
complexify := false:
sextic := proc(ab, section)
 local r, cc, su1, su2, ssub, i;
 global Result;
 Result := true;
 ssub := e -> simplify(subs(ab, e));
 su1 := ssub(bu1);
 su2 := ssub(bu2);
 r := simplify(-ssub(bphi) + ssub(coeff(section, t, 2))*su1^2 +
               ssub(coeff(section, t, 1))*su1*su2 + ssub(coeff(section, t, 0))*su2^2);
 if normalize then
  cc := coeff(r, y^4);
  if cc <> 0 then r := simplify(r/cc); fi;
 fi;
 if complexify then r := evalc(r); fi;
 r := sort(sort(collect(r, [y, x]), [x, alias()]), y);
 return r;
end:

## A section in the form a*(t-mu)^2
parabola := a*(t-mu)^2:

## Further substitutions in "sextic(...)" for deep singularities
#  "mn" is in the form [mu=..., nu=...]
#  It is assumed that global "ab" is stored
AS := proc(mn)
 local ab0, p0;
 ab0 := simplify(get_ab(mn));
 p0 := simplify(get_section(mn));
 return sextic(ab0, p0);
end:

## Sometimes this results in a better looking expression for the sextic
beautify := expr -> collect(op(1, op(1, op(2, factors(expr)))), [y, x]):

#Singularities := [A_1, A_2, A_3, A_4, A_5, A_6, A_7, A_8, A_9, A_10,
# A_11, A_12, A_13, A_14, A_15, A_16, A_17, A_18, A_19]:
str_wrong := ">> Wrong:":
str_OK := "OK:":
str_triple := ">> Triple point":

## Factorization of the discriminant
## Also returnt 'sset_out' and information 'last_t' about the last triple point
discr := proc(expr)
 local d, deg, f, i, n, F, res, validate;
 global Result, sset_out, last_t;
 
 validate := proc(p, m)
  local f, ff, i, x0, y0, cf;
  
  cf := proc(a, b) return simplify(coeff(coeff(f, x, a), y, b)); end;
  
  if m = 1 then return 0; fi;
  x0 := RootOf(subs(x=_Z, p));
  f := simplify(subs(x=x0, expr));
  if degree(f, y) < deg then
   Result := false;
   lp(str_wrong, "degenerate fiber");
   return `?`;
  fi;
  f := op(2, factors(f));
  ff := sort([seq(op(2, op(i, f)), i=1..nops(f))]);
  if 4 in ff then
   Result := false;
   lp(str_wrong, "4-fold tangent/triple point");
   return `?`;
  fi;
  if 3 in ff then
   for i from 1 to nops(f) do
    if op(2, op(i, f)) = 1 then next; fi;
    y0 := solve(op(1, op(i, f)), y);
    f := simplify(subs(x=x+x0, y=y+y0, expr));
    if (cf(0, 1) <> 0) or (cf(1, 0) <> 0) then
     Result := false;
     lp(str_wrong, "3-fold tangent");
     return 0;
    fi;
    if (cf(0, 2) <> 0) or (cf(1, 1) <> 0) or (cf(2, 0) <> 0) then
     Result := false;
     lp(str_wrong, "vertical tangent");
     return `?`;
    fi;
    lp(str_triple);
    f := cf(0, 3)*u^3 + cf(1, 2)*u^2 + cf(2, 1)*u + cf(3, 0);
    f := op(2, factors(f));
    last_t := [true, x0, y0, copy(f)];
    if (nops(f) = 0) or ((nops(f) = 1) and (op(2, op(1, f)) = 3)) then
     return degree(p, x)*E_||(m-2);
    fi;
    return degree(p, x)*D_||(m-2);
   od;
  fi;
  if ff = [2, 2] then
   Result := false;
   lp(str_wrong, "double tangent");
   return `?`;
  fi;
  return degree(p, x)*A_||(m-1);
 end;
 
 F := proc(a, b) return op(2, a) > op(2, b); end;
 
 Result := true;
 sset_out := 0;
 last_t := [false];
 deg := degree(expr, y);
 if degree(expr) < 6 then
  Result := false;
  lp(str_wrong, "reducible curve");
 fi;
 if nops(op(2, factors(expr))) > 1 then
  Result := false;
  lp(str_wrong, "reducible curve");
 fi;
 d := simplify(discrim(expr, y));
 if d = 0 then
  lp(str_wrong, "non-reduced curve");
  Result := false;
  return 0;
 fi;
 f := degree(d, x);
 if deg = 4 then
  sset_out := A_||(25 - f);
 else
  lp(str_triple);
  Result := false;
  n := degree(coeff(expr, y^3), x);
  if n >= 2 then sset_out := D_4;
  elif n = 1 then
   if f = 18 then sset_out := D_5; else sset_out := D_||(22 - f); fi;
  else
   if f = 18 then sset_out := E_6; else sset_out := E_||(23 - f); fi;
  fi;
 fi;
 f := sort(op(2, factors(d)), F);
 F := [seq([op(2, op(i, f)), degree(op(1, op(i, f)), x)], i=1..nops(f))];
 sset_out := sset_out + add(validate(op(op(i, f))), i=1..nops(f));
 res := str_wrong;
 if Result then
  res := str_OK;
#  if sset_out <> sset then Result := false; res := str_wrong; fi;
 fi;
 lp(res, 'sset'=sset_out, 'degree'=degree(d, x), F);
 n := nops(f);
 return mul(sort(collect(op(1, op(i, f)), x), [alias()])^op(2, op(i, f)), i=1..n);
end:

## Convert 'found_s' to a trigonal curve
trigonal := proc()
 local e, cx, i;
 e := op(4, found_t);
 cx := 0;
 for i from 1 to nops(e) do
  if op(2, op(i, e)) > 1 then
   cx := solve(op(1, op(i, e)), u);
   break;
  fi;
 od;
 e := simplify(subs(x=x+op(2, found_t), y=y+op(3, found_t), found_s));
 e := simplify(subs(y=y+cx*x, e));
 e := simplify(subs(y=1/y, x=x/y, e)*y^6);
 if degree(coeff(e, y^3), x) = 0 then e := simplify(e/coeff(e, y^3)); fi;
 if (nargs > 0) and (args[1] = complex) then e := evalc(e); fi;
 return sort(sort(collect(e, [y, x]), x), y);
end:

##############################################################################
# Simplifying and solving equations
##############################################################################

## The sequence of all factors of "expr" (to be placed to 'bad_factors')
all_factors := proc(expr)
 local ff, i;
 ff := op(2, factors(expr));
 return seq(primpart(op(1, op(i, ff))), i=1..nops(ff));
end:

# the list/set of bad factors (to be removed from equations)
bad_factors := []:
# a possible extra filter for 'simplifyed': remove 'expr' if 'false' is returned
return_true := expr -> true:
eqn_filter := return_true:

## Equation "eqn" with all bad factors removed
simplified := proc(eqn)
 local ff, i, p, good;
 good := proc(f)
  if op(2, f) <= 0 then return 1; fi;
  p := primpart(op(1, f));
  if (p in bad_factors) or (-p in bad_factors) then return 1; fi;
  if eqn_filter(p) then return p; fi;
  return 1;
 end:
 ff := simplify(eqn);
 if ff = 0 then return 0; fi;
 ff := op(2, factors(ff));
 return mul(good(op(i, ff)), i=1..nops(ff));
end:

## Takes the output of 'factors' and returns the factors of degree "d" in "t"
the_factors := proc(f, d, t)
 local ff, res, i;
 ff := op(2, f);
 res := [];
 for i from 1 to nops(ff) do
  if degree(op(1, op(i, ff)), t) = d then
   res := [op(res), collect(primpart(op(1, op(i, ff))), t)];
  fi;
 od;
 return res;
end:

## Extracting the (first) factor of a given degree "d" in "t"
the_factor := proc(expr, d, t)
 local res;
 res := the_factors(factors(expr), d, t);
 if nops(res) = 0 then return 1; fi;
 return op(1, res);
end:

## Validate a solution against 'bad_factors'
validate := proc(mn)
 local i;
 for i from 1 to nops(bad_factors) do
  if simplify(subs(mn, op(i, bad_factors))) = 0 then
#   print('solution'=mn, 'violated'=op(i, bad_factors));
   return i;
  fi;
 od;
 return 0;
end:

## Solve and return valid solutions only
solutions := proc(eqn, vars)
 local sol, res, i;
 sol := solve(eqn, vars);
 res := [];
 for i from 1 to nops(sol) do
  if validate(op(i, sol)) = 0 then res := [op(res), op(i, sol)]; fi;
 od;
 return res;
end:

## Construct a sextic from the data found
# this may change
get_sextic := AS:

## Analyze a list of solutions
analyze := proc()
 local sols, sol, s, d;
 
 sol := proc(item)
  print('solution'=item);
  s := get_sextic(item);
  d := discr(s);
  return table(['mn'=item, 'ss'=s, 'dd'=d, 'sset'=sset_out, 'tt'=last_t]);
 end;
 
 sols := args[1];
 if nargs > 1 then sols := solutions(args); fi;
 print2();
 lp("## Solutions found:", nops(sols));
 return [seq(sol(op(i, sols)), i=1..nops(sols))];
end:

## Selects a solution with a given set of singularities
select_solution := proc(list)
 local i;
 global found_mn, found_s, found_d, found_t;
 for i from 1 to nops(list) do
  if op(i, list)['sset'] <> sset then next; fi;
  found_mn := op(i, list)['mn'];
  found_s := op(i, list)['ss'];
  found_d := op(i, list)['dd'];
  found_t := op(i, list)['tt'];
  return true;
 od;
 return false;
end:

##############################################################################
# A few particular equations
##############################################################################

## Coefficients of M_4
M_section := proc()
 global m;
 local e, i;
 e := the_factor(subs(ab, KS(u+mu, a*u^2)), 4, u);
 for i from 0 to 4 do m[i] := factor(coeff(e, u, i)); od;
 return ;
end:

M_simplified := proc()
 global m, ab;
 local i;
 for i from 0 to 4 do m[i] := factor(subs(ab, m[i])): od:
 return ;
end:

## Compute 'a' from the equation 'm_0' = 0
get_a := proc()
 global m, a, a_1;
 a := 'a';
 M_section();
 a_1 := factor(solve(simplified(m[0]), a));
 print('a0'=a_1);
 return a_1;
end:

## M_{4-k}
M := k -> the_factor(add(m[i]*u^(i-k), i=k..4), 4-k, u):

## An extra A_1 point
Eq_A1 := k -> simplified(discrim(M(k), u)):

## An extra A_2 point
Eq_A2 := proc(k)
 local res, e, e1, e2, r1, r2, r3;
 e := M(k); e1 := diff(e, u); e2 := diff(e1, u);
 r1 := simplified(discrim(e, u));
 r2 := simplified(resultant(e, e2, u)/m[4]);
 r3 := simplified(resultant(e1, e2, u)/m[4]);
 res := [];
 # precompute resultants
 if nargs > 1 then
  res := [simplified(resultant(r1, r2, args[2])), simplified(resultant(r1, r3, args[2])),
          simplified(resultant(r2, r3, args[2]))];
  if (nargs > 2) and (args[3] = 'gcd') then
   res := [simplified(gcd(gcd(op(1, res), op(2, res)), op(3, res)))];
  fi;
 fi;
 return op(res), r1, r2, r3;
end:

## Extra degeneration of P_1
Eq_Ak := k -> seq(simplified(m[i]), i=0..k):

## Degeneration (by 1) of P_2 (tangency to y = -S(x))
Eq_P2 := proc()
 return simplified(discrim(the_factor(subs(ab, parabola + S(t, 1)), 2, t), t));
end:

## Degeneration (by 1) of P_3 (tangency to y = -S2(x))
## Assuming S2 defined!
Eq_P3 := proc()
 return simplified(discrim(the_factor(subs(ab, parabola + S2(t, 1)), 2, t), t));
end:

## The condition that a cubic polynomial in "t" is a perfect cube
perfect_cube := proc(expr)
 local c, i;
 c := seq(coeff(expr, t, i), i=0..3);
 return [simplified(3*c[2]*c[4] - c[3]^2), simplified(3*c[1]*c[3] - c[2]^2),
         simplified(9*c[1]*c[4] - c[2]*c[3])];
end:

##############################################################################
# The j-invariant
##############################################################################

## Reduced form of a degree 4 polynomial in 'x'
reduced := proc(f)
 local res;
 res := coeff(f, x^4);
 res := f/res;
 return sort(collect(simplify(subs(x=x-coeff(res, x^3)/4, res)), x), y);
end:

## RThe cubic resolvent of a reduced degree 4 polynomial in 'x'
resolvent := proc(f)
 local p, q, r;
 p := simplify(coeff(f, x^2));
 q := simplify(coeff(f, x));
 r := simplify(coeff(f, x, 0));
 return sort(x^3 - factor(2*p)*x^2 + factor(p^2 - 4*r)*x + factor(q)^2, x);
end:

## The $j$-invariant of a degree 3 polynomial in 'x'
J3 := proc(f)
 local a0, a1, a2, res;
 res := f/coeff(f, x^3);
 a0 := coeff(res, x, 0);
 a1 := coeff(res, x, 1);
 a2 := coeff(res, x, 2);
 return simplify(4/27*(-a2^2+3*a1)^3/(-a2^2*a1^2+4*a1^3+4*a2^3*a0+27*a0^2-18*a0*a1*a2));
end:

## The $j$-invariant of a degree 4 polynomial in 'x'
J4 := f -> J3(resolvent(reduced(f))):

## Just an invariant of a degree 5 polynomial
J5 := proc(f)
 local res;
 res := f/coeff(f, x, 5);
 res := subs(x=x-coeff(res, x, 4)/5, res);
 return simplify(coeff(res, x, 3)^10/discrim(res, x));
end:

## The $j$-invariants of the factors of 'd' of degree 3 and 4
j_invariant := proc(d)
 local ff, i, res;
 ff := op(2, factors(d));
 res := 1;
 for i from 1 to nops(ff) do
  if op(2, op(i, ff)) = 1 then res := res*op(1, op(i, ff)); fi;
 od;
 if degree(res, x) = 3 then return J3(res); fi;
 if degree(res, x) = 4 then return J4(res); fi;
 if degree(res, x) = 5 then return J5(res); fi;
 return 0;
end:

## Check that 'epsilon' (of degree 'deg') can be expressed in terms of
## the j-invariant of 'discr' (when not obvious)
j_check := proc(discr, epsilon, deg)
 local Ji, i, t, eqn, sol;
 Ji := simplify(j_invariant(discr));
 Ji := simplify(evala(add(t[i]*Ji^i, i=0..deg-1) - epsilon));
 eqn := [seq(coeff(Ji, epsilon, i), i=0..deg-1)];
 sol := solve(eqn, [seq(t[i], i=0..deg-1)]);
 if nops(sol) = 0 then return false; fi;
 print(epsilon=subs(op(1, sol), add(t[i]*J^i, i=0..deg-1)));
 return true;
end:
