##!!! Input "common.txt" first
#read "/texts/papers/sextics/equations/Maple/common.txt"; read "/texts/papers/sextics/equations/Maple/A3A1A1.txt";
lp("\
#######################################################################","\n\
# Ramification locus A_3 + 2A_1 \n\
#######################################################################");
u1 := ((beta-alpha)*x[1]+(beta*alpha-2*alpha+1)*x[2])*x[3]^2
     +((beta*alpha*rho-2*alpha*rho+1)*x[1]*x[2]+(alpha-alpha*rho)*x[2]^2)*x[3]
     +(alpha*rho-alpha*rho^2)*x[1]*x[2]^2;
u2 := ((beta*alpha-2*beta+1)*x[1]+(alpha-beta)*x[2])*x[3]^2
     +((beta*alpha*rho-2*beta*rho+1)*x[1]^2+(-2*beta*rho+alpha+alpha*rho)*x[1]*x[2])*x[3]
     +(alpha*rho-beta*rho^2)*x[1]^2*x[2];
init_sextics(u1, u2);

lp("\
## Alternative section");
pp2 := subs(x[1]=1, op(2, solve([u1, u2], [x[2], x[3]])));
ss2 := simplify(alt_section(op(2, op(1, pp2)), op(2, op(2, pp2)))):
S2 := proc(t1, t2) return SS(ss2, t1, t2); end:
lp("\
## The corresponding singular cubic is 'k2*u1 + u2'");
k2 := 'k2':
simplify(subs(pp2, x[1]=1, diff(k2*u1 + u2, x[1]))):
k2 := factor(solve(%, k2));

lp("\
## Prohibited values of the parameters");
lp("\
# One of the cubics splits");
alpha=0, 'u1'=factor(subs(alpha=0, u1));
alpha=1, 'u1'=factor(subs(alpha=1, u1));
beta=0, 'u2'=factor(subs(beta=0, u2));
alpha=beta, 'u1'=factor(subs(alpha=beta, u1));
rho=0, 'u1'=factor(subs(rho=0, u1));
rho=1, 'u1'=factor(subs(rho=1, u1));
lp("\
# The third basepoint degenerates to [0,0,1]");
beta=1, 'pp2'=[1,1,infinity];
lp("\
# The third basepoint degenerates to [1,0,0]");
beta*rho=1, 'pp2'=[1,0,0];

lp("\
#######################################################################","\n\
# A set of singularities adjacent to A_10 + A_4 + A_3 \n\
#######################################################################");
lp("\
## Express 'beta', 'rho' in terms of 'alpha', 'mu', 'nu'");
section0 := the_factor(KS(t, 0), 2, t);
bad_factors := {alpha}:
solutions([subs(t=mu, section0), subs(t=nu, section0)], [beta, rho]);
br := op(1, %); ab := br:
beta0 := factor(subs(br, beta));
rho0 := factor(subs(br, rho));
bad_factors := {}:

lp("\
## Look for a section in the form 'y=a*(t-mu)^2' inflection tangent to 'K'");
a := get_a():

lp("\
## Make this section 'y=a*(t-mu)^2' tangent to 'L_2'");
lp("\
## Equations: ", 'Eq_P2()');
bad_factors := { all_factors(a_1) }:
alpha0 := factor(solve(Eq_P2(), alpha));
beta0 := factor(subs(alpha=alpha0, beta0));
rho0 := factor(subs(alpha=alpha0, rho0));
ab0 := [alpha=alpha0, beta=beta0, rho=rho0];
ab := ab0:
a := factor(subs(ab, a_1)):
lp("\
## Simplify the coefficients (takes time)");
M_simplified();

lp("\
## Invalid values:");
bad_factors := { all_factors(alpha0), all_factors(alpha0 - 1),
  all_factors(beta0), all_factors(beta0 - 1),
  all_factors(rho0), all_factors(rho0 - 1),
  all_factors(a), all_factors(alpha0 - beta0),
  all_factors(beta0*rho0 - 1)
 };
print2();

## Case 1
set_set("\
## The set of singularities", 14=A_12 + A_4 + A_3);
lp("\
## Equations: ", 'Eq_Ak(2)');
select_solution(analyze([Eq_Ak(2)], [mu, nu])):
print_solution();

## Case 2
set_set("\
## The set of singularities", 16=A_11 + 2*A_4);
lp("\
## Equations: ", 'Eq_Ak(1)', 'Eq_P3()');
alias(epsilon2=epsilon2):
select_solution(analyze([Eq_Ak(1), Eq_P3()], [mu, nu])):
print_solution(epsilon2=RootOf(50*_Z^4+300*_Z^3+685*_Z^2+720*_Z+302));
lp("\
## Find a real representation of this curve");
lp("\
## Guess quadratic subextensions");
ifactor(discrim(poly[idx], _Z));
alias(epsilon=epsilon, e2=e2):
set_alias(e2=RootOf(_Z^2-2)):
set_alias(epsilon=RootOf(_Z^2+15)):
factor(poly[idx], {e2, epsilon});
the_factors(factors(subs(_Z=t, poly[idx]), {epsilon, e2}), 1, t):
eps := solve(op(1, %), t):
epsilon2=eps;
snew := evala(subs(epsilon2=eps, found_s)):
dnew := evala(subs(epsilon2=eps, found_d)):
lp("\
## Place the projections of the two A_4 points to 'eps' and '-eps'");
ff := the_factors(factors(dnew), 1, x);
x1 := solve(op(1, ff), x); x2 := solve(op(2, ff), x);
sx := (x1 + x2)/2; cx := evala((x1 - sx)/epsilon);
evala(subs(x=cx*x+sx, dnew)):
lp("\
## The discriminant is real!");
dnew := sort(factor(evala(simplify(%/coeff(%, x, degree(%, x))))), [x, e2]);
snew := evala(simplify(subs(x=cx*x+sx, snew))):
lp("\
## Now, make the 'y'-coordinates of A_4 equal to 0");
y1 := solve(the_factor(subs(x=epsilon, snew), 1, y), y);
y2 := solve(the_factor(subs(x=-epsilon, snew), 1, y), y);
solve([y1=sy+cxy*epsilon, y2=sy-cxy*epsilon], [sy, cxy]):
suby := subs(op(1, %), [y=y+cxy*x+sy]);
snew := simplify(subs(suby, snew)):
lp("\
## Finally, make the coefficients of 'y^4' and 'y^3' equal to 1");
cy := coeff(coeff(snew, y, 3), x, 0);
print_set();
sub2 := [x=cx*x+sx, y=cy*y+cxy*x+sy];
rpoly[idx] := _Z^2-2;
rcurve[idx] := collect(primpart(evala(simplify(subs(y=cy*y, snew/cy^4)))), [y, x]);
dnew := discr(%):
'discr'=%;
'J'=evala(convert(j_invariant(dnew), radical));
snew := 'snew':
dnew := 'dnew':
print2();

lp("\
## Inflection tangency with 'K'");
lp("\
## Equations: ", 'Eq_A2(1, nu, gcd)');
alias(epsilon3=epsilon3):
alias(epsilon21=epsilon21):
alias(epsilon22=epsilon22):
mn := analyze([Eq_A2(1, nu, gcd)], [mu, nu]):
print2();

## Case 3
set_set("\
## Multiple root in 'L_2'", 21=A_10 + A_6 + A_3);
select_solution(mn):
print_solution(epsilon3=RootOf(_Z^2+7), complex);

## Case 4
set_set("\
## Multiple root not in 'L_2'", 25=A_10 + A_4 + A_3 + A_2);
select_solution(mn):
print_solution();

## Case t1
set_set("\
## Triple point:", t1=E_6 + A_10 + A_3);
select_solution(mn):
print_solution(epsilon21=op(2, op(1, found_mn)), radical);
tcurve[idx] := trigonal();
print2();

## Case t2
set_set("\
## Triple point:", t2=D_5 + A_10 + A_4);
select_solution(mn):
print_solution(epsilon22=op(2, op(1, found_mn)));
j_check(found_d, epsilon22, 3):
tcurve[idx] := trigonal();
print2();

lp("\
## Tangency with 'K' and with 'L_3'");
lp("\
## Equations: ", 'Eq_A1(1)', 'Eq_P3()');
alias(epsilon4=epsilon4):
alias(epsilon5=epsilon5):
alias(epsilon23=epsilon23):
alias(epsilon4=epsilon4, epsilon5=epsilon5):
mn := analyze([Eq_A1(1), Eq_P3()], [mu, nu]):
print2();

## Case 5
set_set("\
## Multiple root in 'L_2'", 23=A_10 + A_5 + A_4);
select_solution(mn):
print_solution(epsilon4=RootOf(7*_Z^2+22*_Z-2), radical):
# This representation looks much worse
##alias(epsilon4=RootOf(_Z^2-15));
##mn4 := [mu=-11/7-3/7*epsilon4, nu=-12/7-3/35*epsilon4];

## Case 6
set_set("\
## Multiple root not in 'L_2'", 24=A_10 + 2*A_4 + A_1);
select_solution(mn):
print_solution(epsilon5=RootOf(875*_Z^6+5375*_Z^5+13375*_Z^4+18025*_Z^3+14770*_Z^2+7180*_Z+1592)):
lp("\
## Find a real representation of this curve");
lp("\
## Guess quadratic subextensions");
ifactor(discrim(poly[idx], _Z));
alias(omega=omega, e5=e5):
set_alias(omega=RootOf(_Z^2+55)):
factor(poly[idx], omega);
lp("\
## Find a cubic subextensions using the sum of two conjugate roots");
ff := the_factors(factors(subs(_Z=t, poly[idx]), epsilon5), 1, t);
rt := solve(op(1, ff), t) + solve(op(2, ff), t);
eqn := subs(epsilon5=v, evala(add(x[i]*rt^i, i=0..3))):
mn := solve([seq(coeff(eqn, v, i), i=0..5), x[0]-1], [x[0], x[1], x[2], x[3]]);
rpoly[idx] := sort(primpart(subs(op(1, mn), add(x[i]*_Z^i, i=0..3))));
set_alias(e5=RootOf(rpoly[idx])):
'poly'=factor(poly[idx], {e5, omega});
the_factors(factors(subs(_Z=t, poly[idx]), {e5, omega}), 1, t):
eps := solve(op(1, %), t):
epsilon5=eps;
snew := evala(subs(epsilon5=eps, found_s)):
dnew := evala(subs(epsilon5=eps, found_d)):
lp("\
## Place A_1 to the origin and the parycenter of the A_4's to 1");
ff := op(2, factors(dnew, {e5, omega})):
cx := 0:
for i from 1 to nops(ff) do
 if op(2, op(i, ff)) = 2 then sx := solve(op(1, op(i, ff)), x); fi;
 if op(2, op(i, ff)) = 5 then cx := cx + solve(op(1, op(i, ff)), x); fi:
od;
i := 'i':
'sx'=sx;
cx := (cx - 2*sx)/2;
evala(collect(subs(x=x*cx+sx, dnew), x)):
lp("\
## The discriminant is real!");
dnew := factor(evala(%/coeff(%, x, degree(%, x))));
lp("\
## Place the y-coordinate of A_1 to the origin");
snew := evala(collect(subs(x=x*cx+sx, snew), [y, x])):
ff := op(2, factors(subs(x=0, snew))):
for i from 1 to nops(ff) do
 if op(2, op(i, ff)) = 2 then sy := solve(op(1, op(i, ff)), y); fi;
od;
i := 'i':
'sy'=sy;
snew := evala(collect(subs(y=y+sy, snew), [y, x])):
lp("\
## Make the coefficients of 'y^4' and 'y^3' equal to 1");
cy := coeff(coeff(snew, y, 3), x, 0);
snew := evala(collect(subs(y=y*cy, snew/cy^4), [y, x])):
lp("\
# and the coefficient of 'y^3*x' equal to 0");
cxy := -coeff(coeff(snew, y, 3), x, 1)/4;
snew := evala(collect(subs(y=y+cxy*x, snew), [y, x])):
print_set();
sub5 := [x=cx*x+sx, y=cy*y+cxy*x+sy];
'rpoly'=rpoly[idx];
rcurve[idx] := sort(collect(evala(snew), [y, x]), y);
dnew := discr(%):
'discr'=%;
'J'=simplify(j_invariant(dnew));
print2();
j_check(dnew, e5, 3):
snew := 'snew':
dnew := 'dnew':
print2();

lp("\
#######################################################################","\n\
# A set of singularities adjacent to A_7 + A_6 + A_4 \n\
#######################################################################");
a := 'a': lambda := 'lambda': lambda0 := 'lambda0': bad_factors := {}:
lp("\
## 'K' intersects 'L_3' at 't=lambda'");
L_3 := the_factor(KS(t, -S2(t, 1)), 2, t);
rho0 := factor(solve(subs(t=lambda, L_3), rho));
section2 := subs(rho=rho0, -S2(u+lambda, 1)) + a*u^2:
intr := the_factor(subs(rho=rho0, KS(u+lambda, section2)), 4, u):
for i from 0 to 4 do m[i] := factor(coeff(intr, u, i)): od:
i := 'i':
lp("\
## Make P_3 adjacent to A_6 (triple intersection with 'K')");
a_3 := factor(solve(m[0], a));
a := a_3:
ab := [rho=rho0]:
lp("\
## Simplify the coefficients");
M_simplified():

lp("\
## Invalid values:");
bad_factors := { all_factors(alpha), all_factors(alpha-1),
                 all_factors(beta), all_factors(beta-1),
                 all_factors(rho0), all_factors(rho0-1),
                 all_factors(alpha-beta),
                 all_factors(a_3),
                 # $\\lambda$ cannot equal $k2$
                 all_factors(lambda-subs(rho=rho0, k2)),
                 all_factors(beta*rho0-1) };

lp("\
## New printing routines:");
get_ab := mn -> [alpha=subs(mn, alpha), beta=subs(mn, beta), rho=subs(mn, rho0),
                 lambda=subs(mn, lambda)];
get_section := mn -> subs(get_ab(mn), simplify(subs(u=t-lambda, section2)));

lp("\
## Make the section tangent to 'L_2' (P_2 adjacent to A_4)");
Eq1 := simplified(discrim(section2 + subs(rho=rho0, S(u+lambda, 1)), u));
print2();

## Case 8
set_set("\
## The set of singularities", 30=A_8 + A_7 + A_4);
lp("\
## Make the section tangent to 'L_1' (P_1 adjacent to A_8)");
Eq0 := simplified(discrim(section2, u));
lp("\
## Precompute the resultant of 'Eq1' and 'Eq0'");
res0 := simplified(resultant(Eq1, Eq0, alpha));
lp("\
## Make P_3 adjacent to A_7 and precompute the resultant");
eqn0 := simplified(m[1]):
res1 := simplified(resultant(eqn0, Eq1, alpha));
select_solution(analyze([res0, res1, Eq0, Eq1, eqn0], [beta, lambda, alpha])):
print_solution(complex);

lp("\
## !!! Cheating !!!");
lp("\
## For the rest, ignore univariate factors of degree greater than 4");
eqn_filter := e -> evalb((degree(e, alpha) > 0) or (degree(e, lambda) > 0) or
                         (degree(e, beta) <= 4));

lp("\
## Keep 'Eq0', 'Eq1', and 'res0'");
lp("\
## A 'template' for the discriminant:");
d := discrim(add(n[i]*u^(i-1), i=1..4), u);
d := subs(n=m, d):
lp("\
## Compute the discriminant (takes time)");
eqn0 := simplified(d):
lp("\
## Precompute the resultants");
res1 := simplified(resultant(Eq1, eqn0, alpha)):
g01 := gcd(res0, res1):
res01 := simplified(resultant(res0/g01, res1/g01, lambda));
alias(epsilon7=epsilon7):
alias(epsilon9=epsilon9):
mn := analyze([res01, res0, res1, eqn0, Eq0, Eq1], [beta, lambda, alpha]):

## Case 9
set_set("\
## Multiple root not in 'L_1'", 31=A_8 + A_6 + A_4 + A_1);
select_solution(mn):
print_solution(epsilon9=op(2, op(1, found_mn))*9):
j_check(found_d, epsilon9, 3):
print2():

## Case 7
set_set("\
## Multiple root in 'L_1'", 27=A_9 + A_6 + A_4);
select_solution(mn):
print_solution(epsilon7=op(2, op(1, found_mn))/4);
lp("\
## The line 'x=infinity' is not a component; hence, the curve is irreducible");
j_check(found_d, epsilon7, 3):
print2():

## Case t3
set_set("\
## Triple point:", t3=D_9+A_6+A_4);
select_solution(mn):
print_solution();
tcurve[idx] := found_s;

## Case t4
set_set("\
## Triple point:", t4=D_5+A_8+A_6);
select_solution(mn):
print_solution(complex);
tcurve[idx] := trigonal();

lp("\
#######################################################################","\n\
# Singularities with A_7 and inflection tangency \n\
#######################################################################");
lp("\
## Equations for inflection tangency (perfect cube):");
lp("\
# ", 'eqn1 = 3*m[2]*m[4] - m[3]^2');
eqn1 := simplified(3*m[2]*m[4] - m[3]^2):
lp("\
# ", 'eqn2 = 3*m[1]*m[3] - m[2]^2');
eqn2 := simplified(3*m[1]*m[3] - m[2]^2):
lp("\
# ", 'eqn3 = 9*m[1]*m[4] - m[2]*m[3]');
eqn3 := simplified(9*m[1]*m[4] - m[2]*m[3]):
lp("\
## Eliminate 'lambda'");
lp("\
# ", 'resultant(Eq1, eqn1, lambda)');
res1 := simplified(resultant(Eq1, eqn1, lambda)):
lp("\
# ", 'resultant(Eq1, eqn3, lambda)');
res3 := simplified(resultant(Eq1, eqn3, lambda)):
lp("\
# ", 'resultant(Eq1, eqn2, lambda)');
res2 := simplified(resultant(Eq1, eqn2, lambda)):
lp("\
## Resultants with respect to 'alpha' and their 'gcd' (takes time)");
g13 := gcd(res1, res3);
res13 := simplified(resultant(res1/g13, res3/g13, alpha));
g12 := gcd(res1, res2);
res12 := simplified(resultant(res1/g12, res2/g12, alpha));
g23 := gcd(res2, res3);
res23 := simplified(resultant(res2/g23, res3/g23, alpha));
res := simplified(gcd(gcd(res13, res12), res23));
lp("\
## Separate quadratic factors with real roots\n\
## Otherwise, solutions are lost!");
res := op(2, factors(res)):
res_2 := 1: res_4 := 1:
for i from 1 to nops(res) do
 temp := op(1, op(i, res)):
 if (degree(temp) = 2) and (discrim(temp, beta) > 0) then
  res_2 := res_2*primpart(temp, beta):
  next:
 fi:
 res_4 := res_4*primpart(temp, beta):
od:
i := 'i': temp := 'temp':

## Case 10
set_set("\
## The set of singularities", 35=A_7 + A_6 + A_4 + A_2);
alias(epsilon10=epsilon10):
'res_2'=res_2;
select_solution(analyze([res_2, res1, res2, res3, Eq1, eqn1, eqn2, eqn3], [beta, alpha, lambda])):
print_solution(epsilon10=op(2, op(1, found_mn)), radical);

alias(epsilon11=epsilon11):
'res_4'=res_4;
mn := analyze([res_4, res1, res2, res3, Eq1, eqn1, eqn2, eqn3], [beta, alpha, lambda]):

## Case 11
set_set("\
## The set of singularities", 34=A_7 + 2*A_6);
select_solution(mn):
print_solution(epsilon11=op(2, op(1, found_mn)));
'J'=simplify(evala(evalc(convert(j_invariant(found_d), radical))));
lp("\
## This curve is defined over", sqrt(-7));
lp("\
## However, the equations are too bad!!");
print2();

## Case t5
set_set("\
## Triple point:", t5=E_6+A_7+A_6);
select_solution(mn):
print_solution(complex);
tcurve[idx] := trigonal(complex);
