# This function takes a polynomial f(x,y) quartic in y and with constant
# leading coefficient and returns:
# - the reduced form fred := y^4 + p(x)y^2 + q(x)y + r(x);
# - the cubic resolvent res := y^3 - 2py^2 + (p^2-4r)y + q^2
# and its coefficient b1 := p^2-4r;
# - the coefficients g2, g3 in the reduced resolvent red := y^3 + g2y + g3.
#resolvent := proc(f)
# global fred, res, red, b1, g2, g3, d, p, q, r;
## d := factor(discrim(f, y)); print(%);
# fred := coeff(f, y^4);
# if degree(fred, x) <> 0 then error "erDegree" " Non-constant leading term"; fi;
# fred := f/fred;
# fred := simplify(subs(y=y-coeff(fred, y^3)/4, fred));
# p := factor(coeff(fred, y^2));
# q := factor(coeff(fred, y));
# r := factor(coeff(fred, y, 0));
# b1 := factor(p^2 - 4*r);
# res := y^3 - 2*p*y^2 + b1*y + q^2;
# red := simplify(subs(y=y+2*p/3, res));
# g2 := factor(coeff(red, y));
# g3 := factor(coeff(red, y, 0));
# return ;
#end:

reduced := proc(f)
 local res;
 res := coeff(f, y^4);
 if degree(res, x) <> 0 then error "erDegree" " Non-constant leading term"; fi;
 res := f/res;
 return sort(sort(collect(simplify(subs(y=y-coeff(res, y^3)/4, res)), [y, x]), x), y);
end:

resolvent := proc(f)
 local p, q, r;
 if (simplify(coeff(f, y^4)) <> 1) or (simplify(coeff(f, y^3)) <> 0) then
  error "erReduced" " Not a reduced polynomial";
 fi;
 p := simplify(coeff(f, y^2));
 q := simplify(coeff(f, y));
 r := simplify(coeff(f, y, 0));
 return sort(sort(y^3 - factor(2*p)*y^2 + factor(p^2 - 4*r)*y + factor(q)^2, x), y);
end:

## The $j$-invariant of a trigonal curve f(x, y)
j_invariant := proc(f)
 local a0, a1, a2;
 a0 := coeff(f, y, 0);
 a1 := coeff(f, y, 1);
 a2 := coeff(f, y, 2);
 return factor(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 the cubic resolvent of a tetragonal curve f(x, y)
J4 := f -> j_invariant(resolvent(reduced(f))):

## Analyze the $j$-invariant of a maximizing sextic
J := proc(f)
 local ff, m, j, j0, j1, jx, d, i, res, red, is_in;
 
 red := proc(f)
  local res, i;
  res := op(2, factors(evala(f)));
  return [seq([primpart(op(1, op(i, res))), op(2, op(i, res))], i=1..nops(res))];
 end;
 
 is_in := proc(f, m, list)
  local i;
  for i from 1 to nops(list) do
   if (f = op(1, op(i, list))) and (m = op(2, op(i, list))) then return true; fi;
  od;
  return false;
 end;
 
 j := evala(J4(f));
 j0 := red(j);
 j1 := red(j-1);
 d := red(1/discrim(f, y));
 jx := red(diff(j, x));
 res := 1;
 for i from 1 to nops(jx) do
  ff := op(1, op(i, jx));
  m := op(2, op(i, jx));
  if is_in(ff, m+1, d) then next; fi;
  if is_in(ff, m+1, j0) then next; fi;
  if is_in(ff, m+1, j1) then next; fi;
  res := res*ff^m;
 od;
 if degree(res, x) = 1 then
  j0 := [-coeff(res, x, 0)/coeff(res, x, 1)];
 else
  j0 := [solve(res, x)];
 fi;
 j0 := [seq(evala(subs(x=op(i, j0), j)), i=1..nops(j0))];
 return res, j0;
end:

## Automated processing of the j-invariant
## Extra parameters: epsilon, deg to be expressed in terms of j
JJ := proc(f)
 local res, eps, deg, i, t, eqn, sol;
 res := J(f)[2];
 print(res);
 print('points'=evalp(evalf(res)));
 if nops(res) <> 1 then
  lp(">> Wrong:", 'degree'=nops(res));
  return false;
 fi;
 if nargs = 1 then return; fi;
 eps := args[2];  deg := args[3];
 res := op(1, res);
 res := simplify(evala(add(t[i]*res^i, i=0..deg-1) - eps));
 eqn := [seq(coeff(res, eps, i), i=0..deg-1)];
 sol := solve(eqn, [seq(t[i], i=0..deg-1)]);
 if nops(sol) = 0 then lp(">> Wrong:", "smaller field"); return; fi;
 print(eps=subs(op(1, sol), add(t[i]*J^i, i=0..deg-1)));
 return;
end:

#########################################################################
# These are auxiliary functions for 'fibers'
#########################################################################

lp := proc() printf("%A %q\n", args); end:

## Sort a list of tables according to Re(l['val'])
sort_values := proc(list)
 local comp;
 comp := proc(l1, l2) return evalb(Re(l1['val']) < Re(l2['val'])); end;
 return sort(list, comp);
end:

## A more accurate evaluation
evaluate := e -> evalf(evala(e)):

## All roots of an irreducible polynomial f(t)
all_roots := proc(f, t)
 return seq(RootOf(subs(t=_Z, f), index=i), i=1..degree(f, t));
end:

## Evaluated sorted roots (of a reducible polynomial), separately real and complex
all_root_values := proc(f, t)
 local ff, roots, mult, rec, i, j, res_R, res_C;
 ff := op(2, factors(f));
 res_R := [];
 res_C := [];
 for i from 1 to nops(ff) do
  mult := op(2, op(i, ff));
  roots := [all_roots(op(1, op(i, ff)), t)];
  for j from 1 to nops(roots) do
   rec := table(['val'=evaluate(op(j, roots)), 'm'=mult, 'rt'=op(j, roots)]);
   if Im(rec['val']) <> 0 then
    res_C := [op(res_C), copy(rec)];
   else
    res_R := [op(res_R), copy(rec)];
   fi;
  od;
 od;
 return sort_values(res_R), sort_values(res_C);
end:

## The fiber of f(x, y) at x=x0
# If a third argument is given (a table), values are set
fiber := proc(f, x0)
 local res, dest, res_R, res_C, sz, i;
 res := all_root_values(subs(x=x0, f), y);
 res_R := [seq(op(i, res[1])['val'], i=1..nops(res[1]))];
 res_C := [seq(op(i, res[2])['val'], i=1..nops(res[2]))];
 res := [seq(op(i, res[1])['m'], i=1..nops(res[1]))];
 sz := sum(op(i, res), i=1..nops(res));
 if res = [2] then
 # check the position of the imaginary branches
  i := op(1, res_R) - Re(op(1, res_C));
  if i > 0 then res := [`-`, 2]; fi;
  if i < 0 then res := [2, `-`]; fi;
 fi;
 if nargs > 2 then
  dest := args[3];
  dest['sign'] := res;
  dest['R'] := res_R;
  dest['C'] := res_C;
  dest['size'] := sz;
 fi;
 return 'sign'=res, 'R'=res_R, 'C'=res_C, 'size'=sz;
end:

## Check that a sorted list of complex numbers splits into conjugate pairs
err_str := table([true=OK, false=Error]):
check_complex := proc(list)
 local i;
 if nops(list) mod 2 <> 0 then return false; fi;
 for i from 1 to nops(list)/2 do
  if Im(op(2*i-1, list) + op(2*i, list)) <> 0 then return false; fi;
 od;
 return true;
end:
 
printing := table(['fibers'=true,   # print fibers
                   'complex'=false,  # details on complex fibers
                   'values'=false,   # print critical values
                   'digits'=3,       # round up the values
                   'GLE'=false,     # print GLE input
                   'GAP'=true       # print GAP code
                   ]):
evalp := proc(val)
 local res;
 res := val;
 if (printing['digits'] > 0) and (res <> infinity) then
  res := evalf[printing['digits']](res);
 fi;
 return res;
end:

print_fibers := proc(rec)
 local ff, res_R, res_C, res_m, i, print_m, print_R, print_C, vals, value;

 vals := proc(rec)
  if printing['values'] = true then
   return err_str[check_complex(rec['C'])], rec['R'], rec['C'];
  fi;
  return err_str[check_complex(rec['C'])];
 end;
 
 value := rec -> evalp(rec['val']);
 
 print_m := proc(rec)
  lp("========", value(rec), rec['sign'], vals(rec));
  if rec['sign'] = [] then lp("%% info:", 'twist'=rec['twist']); fi;
 end;
 
 print_R := proc(rec)
  lp(">>>>", rec['m']=value(rec), rec['sign'], vals(rec));
  if rec['val'] = 0 then
   lp("%% info:", 'real'=rec['real']);
  fi;
 end;
 
 print_C := proc(rec) lp(">>>>", rec['m']=value(rec)); end;

 res_R := rec['R'];
 res_C := rec['C'];
 res_m := rec['mid'];
 lp("## Real fibers:", 'fibers'=nops(res_R));
 for i from 1 to nops(res_R) do
  print_m(op(i, res_m));
  print_R(op(i, res_R));
 od;
 if nops(res_m) > nops(res_R) then print_m(op(nops(res_m), res_m)); fi;
 lp("## Imaginary fibers:", 'fibers'=nops(res_C),
        err_str[check_complex([seq(op(i, res_C)['val'], i=1..nops(res_C))])]);
 if printing['complex'] = true then
  for i from 1 to nops(res_C) do print_C(op(i, res_C)); od;
 fi;
end:

print_GLE := proc(rec)
 global format;
 lp("#### GLE input:");
 format := format_GLE;
 printf("%A", cat(
        format['null'],
        seq(GLE_point(rec, i), i=1..nops(rec['R'])),
        format['null'],
        seq(sprintf(format['beta_line'], beta_braid(rec, i)), i=1..nops(rec['R'])),
        format['null'],
        seq(sprintf(format['gamma_line'], gamma_braid(rec, i)), i=1..nops(rec['mid']))
        ));
end:

print_GAP := proc(rec)
 global format;
 lp("#### GAP code:");
 format := format_GAP;
 printf("%A", cat(
        sprintf(format['beta_line'], seq(beta_braid(rec, i), i=1..nops(rec['R']))),
        sprintf(format['gamma_line'], seq(gamma_braid(rec, i), i=1..nops(rec['mid']))),
        sprintf(format['mu_line'], seq(mu_braid(rec, i), i=1..nops(rec['R']))),
        kappa_braid(rec),
        format['init'],
        format['process'][evalb(rec['torus']=true)],
        ""));
end:

invert := e -> simplify(subs(x=1/x, y=y/x^2, e)*x^8):

## Neighboring fibers
next_m := proc(rec, idx, shift)
 return op(((idx+shift-1) mod nops(rec['mid'])) + 1, rec['mid']);
end:

next_R := proc(rec, idx, shift)
 return op(((idx+shift-1) mod nops(rec['R'])) + 1, rec['R']);
end:

## The position of the `essential' singular point
position_sing := proc(fiber)
 local sg;
 sg := fiber['sign'];
 if (sg = [2, 1, 1]) or (sg = [2, `-`]) or (sg = [2, `*`]) then return 3; fi;
 if sg = [1, 2, 1] then return 2; fi;
 if (sg = [1, 1, 2]) or (sg = [`-`, 2]) or (sg = [`*`, 2]) then return 1; fi;
 return 0;
end:

## The position of the imaginary branch (assuming there is one)
position_img := proc(fiber)
 local sg;
 sg := fiber['sign'];
 # if there is an imaginary branch, return it
 if sg = [2, `-`] then return 1; fi;
 if sg = [`-`, 2] then return 3; fi;
 # two singular points: depends on which one has real branches
 if sg = [2, `*`] then
  if fiber['real'] then return 3; else return 1; fi;
 fi;
 if sg = [`*`, 2] then
  if fiber['real'] then return 1; else return 3; fi;
 fi;
 # otherwise, the imaginary branches are at the singular point
 return position_sing(fiber);
end:

# formats
format_GLE := table([
 'times' = "",
 'id' = "\\id",
 'beta' = "\\Gs_%d^{-%d}",
 'gamma' = "\\Gs_%d\\1\\Gs_%d",
 'tau' = "\\tau^{%d}",
 'beta_line' = " cbeta _next() \"%s\"\n",
 'gamma_line' = " cgamma _next()-1 \"%s\"\n",
 'point' = " cpt 0 %d \"%d\" 0\n",
 'null' = " _pos=0\n"
]):
format_GAP := table([
 'times' = "*",
 'id' = "\\1",
 'beta' = "s%d^-%d",
 'gamma' = "s%d^-1*s%d",
 'tau' = "tau^%d",
 'mu' = "s%d^%d",
 'kappa' = "K := %d;;\nkappa := a%d*a%d;;\n",
 'beta_line' = "beta := [%Q];;\n",
 'gamma_line' = "gamma := [%Q];;\n",
 'mu_line' = "braid := [%Q];;\n",
 'init' = "init_group(0: Simplified);\n",
 'process' = table([true = "qsize(g.1/g.2);\nsimplify(g.1/g.2);\n", # torus type
                    false = "Size(g);\n"])                          # not torus type
]):
format := format_GLE:

smart_cat := proc()
 local res, ins, i;
 res := args[1];
 for i from 2 to nargs do
  ins := "";
  if (res <> "") and (args[i] <> "") then ins := format['times']; fi;
  res := cat(res, ins, args[i]);
 od;
 return res;
end:

## Computing the braid at a singular fiber
beta_braid := proc(rec, idx)
 local f, f0, f1, s, s0, s1, mult, pos, str0, str;
 
 str := proc(s)
  if s > 0 then return smart_cat(sprintf(format['beta'], pos, s), str0); fi;
  if str0 <> "" then return str0; fi;
  return format['id'];
 end;
 
 f := op(idx, rec['R']); s := f['size']; mult := f['m'];
 f0 := next_m(rec, idx, 0); s0 := f0['size'];
 f1 := next_m(rec, idx, 1); s1 := f1['size'];
 pos := position_sing(f);
 if pos = 0 then return err_str[false]; fi;
 str0 := "";
 if f['rt'] = 0 then
  if pos = 2 then return err_str[false]; fi;
  mult := mult - 2;
  str0 := sprintf(format['beta'], 4-pos, 1);
 fi;
 if mult mod 2 = 0 then
  if s0 <> s1 then return err_str[false]; fi;
  return str(mult/2);
 fi;
 if s1 - s0 = 2 then s := 0
 elif s0 - s1 = 2 then s := 1
 else return err_str[false]; fi;
 s := s + (mult-1)/2;
 return str(s);
end:

## Computing the braid between singular fibers
gamma_braid := proc(rec, idx)
 local f, s, p0, p1, i;
 f := op(idx, rec['mid']); s := f['size'];
 if s = 4 then return format['id']; fi;
 if s = 0 then
  s := f['twist'];
  if s = 0 then return format['id']; fi;
  return sprintf(format['tau'], s);
 fi;
 p0 := position_img(next_R(rec, idx, -1));
 p1 := position_img(next_R(rec, idx, 0));
 if (p0 = 0) or (p1 = 0) then return err_str[false]; fi;
 if p0 = p1 then return format['id']; fi;
 if p0 < p1 then
  return smart_cat(seq(sprintf(format['gamma'], i+1, i), i=p0..p1-1));
 fi;
# return smart_cat(seq(sprintf(format['gamma'], i-1, i), i=p1+1..p0));
 return smart_cat(seq(sprintf(format['gamma'], -i-1, -i), i=-p0..-p1-1));
end:

## Computing the 'mu' braid at a singular fiber
mu_braid := proc(rec, idx)
 local f, mult, pos, str0;
 f := op(idx, rec['R']); mult := f['m'];
 pos := position_sing(f);
 if pos = 0 then return err_str[false]; fi;
 str0 := "";
 if f['rt'] = 0 then
  if pos = 2 then return err_str[false]; fi;
  mult := mult - 2;
  str0 := sprintf(format['mu'], 4-pos, 2);
 fi;
 return smart_cat(sprintf(format['mu'], pos, mult), str0);
end:

## Details on the distinguished fiber for GAP
kappa_braid := proc(rec)
 local i, f, pos;
 for i from 1 to nops(rec['R']) do
  f := op(i, rec['R']);
  if f['rt'] <> 0 then next; fi;
  pos := position_sing(f);
  if (pos <> 1) and (pos <> 3) then break; fi;
  return sprintf(format['kappa'], i, pos, pos+1);
 od;
 return err_str[false];
end:

## Singular point for GLE input
GLE_point := proc(rec, idx)
 local f, w, mu;
 f := op(idx, rec['R']); mu := f['m'] - 1;
 w := -1;
 if f['rt'] = 0 then w := 0; mu := mu - 2; fi;
 return sprintf(format_GLE['point'], w, mu);
end:

# we want very high precision!
Digits := 60:
distance := 4:
fibers := proc(expr)
 local new_expr, d, q, b1, ff, res_R, res_C, res_m, i, j, q0, q1;

# invert 'x'
 new_expr := invert(expr);
# collect the roots of 'q'
 d := reduced(new_expr);
 q := coeff(d, y, 1); q0 := q; q1 := q;
 b1 := coeff(resolvent(d), y, 1);
 # remove roots common with 'b1'
 while degree(gcd(q0, b1, 'd'), x) > 0 do q0 := d; od;
 ff := all_root_values(q0, x)[1];
 q := [];
 for i from 1 to nops(ff) do
 # select only those of odd multiplicity and with 'b1' > 0
  if (op(i, ff)['m'] mod 2 <> 0) and (evaluate(subs(x=op(i, ff)['rt'], b1)) > 0) then
   d := simplify(subs(x=u+op(i, ff)['rt'], q1));
   d := evaluate(coeff(d, u, op(i, ff)['m']));
   q := [op(q), table(['val'=op(i, ff)['val'], 'sign'=sign(d)])];
  fi;
 od;
# find the critical values
 d := simplify(discrim(new_expr, y));
 ff := all_root_values(d, x);
 res_R := ff[1];
 res_C := ff[2];
# compute the fibers at the real critical values
 for i from 1 to nops(res_R) do
  ff := op(i, res_R);
  fiber(new_expr, ff['rt'], ff);
  if ff['rt'] <> 0 then next; fi;
 # check if the branches at the node (0, 0) are real
  b1 := coeff(coeff(new_expr, y, 2), x, 0)*t^2 + coeff(coeff(new_expr, y, 1), x, 1)*t
        + coeff(coeff(new_expr, y, 0), x, 2);
  ff['real'] := evalb(evaluate(discrim(b1, t)) > 0);
 # position of the node (0, 0)
  for j from 1 to nops(ff['R']) do
   if op(j, ff['R']) = 0 then
    ff['sign'] := subsop(j=`*`, ff['sign']);
    break;
   fi;
  od;
 od;
# select intermediate points and compute their fibers
 ff := [op(res_R), op(res_C)];
 ff := sort([seq(Re(op(i, ff)['val']), i=1..nops(ff))]);
 res_m := [op(1, ff) - distance,
           seq((op(i, res_R)['val'] + op(i+1, res_R)['val'])/2, i=1..nops(res_R)-1),
           op(nops(ff), ff) + distance];
 res_m := [seq(table(['val'=op(i, res_m),
                      fiber(new_expr, convert(op(i, res_m), fraction))]), i=1..nops(res_m))];
# compute the twist parameters
 for i from 1 to nops(res_m) do
  ff := op(i, res_m);
  if ff['sign'] <> [] then next; fi;
  if i = 1 then q0 := -infinity; else q0 := op(i-1, res_R)['val']; fi;
  if i = nops(res_m) then q1 := infinity; else q1 := op(i, res_R)['val']; fi;
  b1 := 0;
  for j from 1 to nops(q) do
   if (op(j, q)['val'] > q0) and (op(j, q)['val'] < q1) then
    b1 := b1 + op(j, q)['sign'];
   fi;
  od;
  ff['twist'] := b1;
 od;
# the fiber at infinity
 if degree(d, x) < 24 then
  res_R := [op(res_R), table(['val'=infinity, 'm'=24-degree(d, x),
                             fiber(expr, 0), 'rt'=infinity])];
 fi;
 ff := table(['R'=res_R, 'C'=res_C, 'mid'=res_m, args[2..nargs]]);
# print the results
 if printing['fibers'] = true then print_fibers(ff); fi;
 if printing['GLE'] = true then print_GLE(ff); fi;
 if printing['GAP'] = true then print_GAP(ff); fi;
 return ff;
end:
