##############################################################################
# Geiser involution by Alex Degtyarev
#
# This is a Maple implementation of E. Moody's formulas for the Geiser involution
# and my extension for the associated map P^2 --> P^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]
# Important: f1 must be singular at (0:0:1) !!!
##############################################################################

## 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
# w1 is singular at (0:0:1)
a1[1] := 0:
a1[2] := 0:

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

# 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]*(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):
 wy := subs(x=y, w1):

 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:
 gamma1 := -(a2[1]*y[1] + a2[2]*y[2]):
 C2 := - a2[2]*dv(B[1]-a2[1]*b1[1]*y[1]*y[3]^2, y[2])
  + a2[1]*dv(a2[2]*y[3]*(wy - b1[1]*y[1]^2*y[3]) -  B[3]*y[2], y[1]*y[2]):
 phi3 := -a2[1]*C[2] + y[3]*C2:
 psi3 := -a2[2]*C[1] + y[3]*C2:
# these are $\bar r'_1$, $\bar r'_3$
 rr[1] := B[1]*a2[2]^2 - B[2]*a2[1]*a2[2] + B[3]*a2[1]^2:
 rr[3] := -a2[2]*C[1] + a2[1]*C[2]:
# and these are $\bar r_1$, $\bar r_2$, $\bar r_3$
 Rr[1] := -a2[2]*rr[1]:
 Rr[2] := a2[1]*rr[1]:
 Rr[3] := a2[1]*a2[2]*rr[3]:
# the Bertini involution
 yb[1] := phi3*dv(a2[2]^2*phi3*wy + B[3]*rr[1], y[1]):
 yb[2] := psi3*dv(a2[1]^2*psi3*wy + B[1]*rr[1], y[2]):
 yb[3] := psi3*phi3*C2:
# the coefficients in psi6 = phi6 + s[2]*w1^2 + s[1]*w1*w2 + s[0]*w2^2:
 s[0] := 0:
 s[1] := a2[1]*c1[2]-a2[2]*c1[1]:
 s[2] := a2[2]*c2[1]-a2[1]*c2[2]:
# the fixed point set
 K := psi3*dv(-a2[1]*wy*y[3] + B[1]*y[1], y[2]) -  phi3*dv(-a2[2]*wy*y[3] + B[3]*y[2], y[1]):
# the coefficients in the equation for K^2:
 r[0] := 0:
 r[1] := a2[1]*b1[2]*c1[2]-a2[1]*b1[3]*c1[1]-a2[2]*b1[1]*c1[2]:
 r[2] := -a2[1]*b2[2]*c1[2]-a2[1]*b1[2]*c2[2]+a2[1]*b2[3]*c1[1]
         +a2[1]*b1[3]*c2[1]+a2[2]*b2[1]*c1[2]+a2[2]*b1[1]*c2[2]:
 r[3] := a2[1]*b2[2]*c2[2]-a2[1]*b2[3]*c2[1]-a2[2]*b2[1]*c2[2]:

 p[0] := b1[2]^2-4*b1[1]*b1[3]:
 p[1] := -8*a2[1]*c1[2]+4*a2[2]*c1[1]+4*b1[1]*b2[3]+4*b2[1]*b1[3]-2*b1[2]*b2[2]:
 p[2] := b2[2]^2-4*a2[2]*c2[1]-4*b2[1]*b2[3]+8*a2[1]*c2[2]:
 
 q[0] := 0:
 q[1] := -4*a2[1]*b1[1]*b1[3]*c1[2]+2*a2[1]*b1[2]^2*c1[2]-2*a2[2]*b1[1]*b1[2]*c1[2]
         -2*a2[1]*b1[2]*b1[3]*c1[1]+4*a2[2]*b1[1]*b1[3]*c1[1]:
 q[2] := -4*a2[1]^2*c1[2]^2+4*a2[1]*a2[2]*c1[1]*c1[2]+4*a2[1]*b1[1]*b1[3]*c2[2]
         +4*a2[1]*b1[1]*b2[3]*c1[2]+4*a2[1]*b2[1]*b1[3]*c1[2]-2*a2[1]*b1[2]^2*c2[2]
         -4*a2[1]*b1[2]*b2[2]*c1[2]+2*a2[1]*b1[2]*b1[3]*c2[1]+2*a2[1]*b1[2]*b2[3]*c1[1]
         +2*a2[1]*b2[2]*b1[3]*c1[1]+2*a2[2]*b1[1]*b1[2]*c2[2]+2*a2[2]*b1[1]*b2[2]*c1[2]
         +2*a2[2]*b2[1]*b1[2]*c1[2]-4*a2[2]*b1[1]*b1[3]*c2[1]-4*a2[2]*b1[1]*b2[3]*c1[1]
         -4*a2[2]*b2[1]*b1[3]*c1[1]:
 q[3] := -4*a2[1]*a2[2]*c1[1]*c2[2]+4*a2[1]*b2[2]*b1[2]*c2[2]-2*a2[2]*b2[1]*b1[2]*c2[2]
         -2*a2[1]*b1[2]*b2[3]*c2[1]+4*a2[2]*b2[1]*b1[3]*c2[1]+4*a2[2]*b1[1]*b2[3]*c2[1]
         -4*a2[1]*b2[1]*b1[3]*c2[2]-4*a2[1]*b1[1]*b2[3]*c2[2]-2*a2[2]*b1[1]*b2[2]*c2[2]
         -4*a2[1]*a2[2]*c2[1]*c1[2]-2*a2[1]*b2[2]*b1[3]*c2[1]-4*a2[1]*b2[1]*b2[3]*c1[2]
         -2*a2[2]*b2[1]*b2[2]*c1[2]-2*a2[1]*b2[2]*b2[3]*c1[1]+4*a2[2]*b2[1]*b2[3]*c1[1]
         +2*a2[1]*b2[2]^2*c1[2]+8*a2[1]^2*c2[2]*c1[2]:
 q[4] := 4*(a2[1]*c2[2] - b2[1]*b2[3])*s[2] - 2*b2[2]*r[3]:
 
## Just check the identities to make sure that everything has initialized properly
check_data := proc()
 print("The expansion of $\\psi_3$");
 print(simplify(psi3 - subs(x=y, phi3 + s[2]*w1 + s[1]*w2)));
 print("The expansion of $K^2$ (coefficientwise)");
 print(simplify(- K^2 - 4*phi3^3*wy + subs(x=y,
  + phi3^2*(p[2]*w1^2 + p[1]*w1*w2 + p[0]*w2^2)
  + phi3*(q[4]*w1^3 + q[3]*w1^2*w2 + q[2]*w1*w2^2 + q[1]*w2^3)
  + (r[3]*w1^2 + r[2]*w1*w2 + r[1]*w2^2)^2)));
 print("The identities for $r$, $\\phi$, and $\\psi$");
 print(simplify(y[1]*Rr[3] - y[3]*Rr[1] - a2[2]*gamma1*phi3));
 print(simplify(y[2]*Rr[3] - y[3]*Rr[2] + a2[1]*gamma1*psi3));
 print("The divisibility of $y_3 z_1 - y_1 z_3$ by $K$");
 print(simplify((y[3]*yb[1]-y[1]*yb[3])/phi3/K - a2[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])/psi3/K + a2[1]));
 print("All done");
end:
