##############################################################################
# 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:
