# A. Galligo and M. van Hoeij
#
# Version 0.01,  Jan 2007.

# Weaknesses in current version:
#
#   Program is slow, and if epsilon, the size of the errors in the
#   coefficients, is not small then the result may depend on the
#   projection used (i.e. interchanging x, y may change the result).  A
#   more sophisticated look at the monodromy data would give better
#   quality results even for relatively large epsilon (in particular,
#   if we have 3 nearby branchpoints, the code should consider triples
#   of branchpoints, instead of only pairs which it does at the moment).
#
#   Reconstruction of factors from the monodromy information is not yet
#   implemented, so the only information obtained inside the program is
#   the number of factors and the degree and genus of each factor (and of
#   this information, the only thing returned in the output is whether or
#   not the monodromy was transitive).


# Some macro's needed for the programs below:
macro(
	invperm = group[invperm],
	orbit = group[orbit]
):

mp := proc(a,b)
	options remember;
	if nargs=2 then
		group[mulperms](b,a)
	else
		mp(a,mp(args[2..-1]))
	fi
end:

# Input: M, the output of algcurves[monodromy].
# Note:  M[1] = basepoint.
#        M[2] = fiber above the basepoint
#        M[3] = list of branchpoints with permutations.
#
# p,q: two integers 1 < p < q < nops(M[3]) representing the
# branchpoints M[3][p] and M[3][q].
#
# Return true if the shortest path around p,q has permutation = [].
#
MatchingPermutations := proc(M,p,q)
	local g,i;
	g := [];
	for i in PointsInTriangle(M,p,q) do
		g := mp(M[3,i,2], g)
	od;
	evalb( [] = mp(invperm(g), M[3,q,2], g, M[3,p,2]) )
end:

# Find the points inside the triangle basepoint-p-q
#
PointsInTriangle := proc(M,p,q)
	options remember;
	local L,t,i;
	# The line through p and q:
	L := t * M[3,p,1] + (1-t) * M[3,q,1];
	[seq(`if`(Im(solve(L=M[1],t))*Im(solve(L=M[3,i,1],t)) > 0,
	i, NULL), i=p+1..q-1)]
end:

# M = monodromy
# P = the subset of {1,2,..nops(M[3])} that we keep
# Check if the monodromy is transitive with those points P.
#
TransitiveGroup := proc(M, P)
	local G,i;
	G := permgroup(nops(M[2]), {seq(M[3,i,2], i=P)});
	evalb( nops(orbit(G, 1)) = nops(M[2]) )
end:

ConjugatePerm := proc(s, g)
	mp(invperm(g), s, g)
end:

TransitiveMonodromy := proc(F, x, y)
	local Fr, M, S, p, q, d, P, M_upd, i, TR, per, j, u;
	# Round coefficients to integers.
	Fr := primpart(collect(expand(10^5*F),{x,y},round));

	M := algcurves[monodromy](Fr,x,y);

	# For some reason this goes slow without the line "if d > ..."\
	# even though it should be a trivial computation:
	S := NULL;
	for p to nops(M[3])-1 do for q from p+1 to nops(M[3]) do
		d := abs(M[3,p,1]-M[3,q,1]);
		if d > 0.3 then next fi;
		if MatchingPermutations(M,p,q) then
			S := S, [p,q, d]
		fi
	od od:

	S := sort([S], (i,j) -> i[3]< j[3]);
	nops(S);

	# Now we can start throwing away grouped branchpoints, and we'll
	# stop as soon as the group becomes reducible.

	P := {seq(i,i=1..nops(M[3]))};
	M_upd := M: # M with updated permutations

	for i in S while TransitiveGroup(M_upd,P) do
		if not (member(i[1],P) and member(i[2],P)) then
			# already threw i[1] or i[2] away, can't do it again:
			next
		fi;
		P := P minus {i[1],i[2]};
		TR := PointsInTriangle(M,i[1],i[2]);
		per := M_upd[3,i[1],2];
		for j from i[1]+1 to i[2]-1 do if member(j,P) then
			if member(j, TR) then
				per := ConjugatePerm(per, invperm(M_upd[3,j,2]))
			else
				u := ConjugatePerm( M_upd[3,j,2], per);
				M_upd := subsop(3 = subsop(j =
				  [M_upd[3,j,1], u], M_upd[3]), M_upd)
			fi
		fi od;
		if per <> M_upd[3,i[2],2] then
			error nonmatch
		fi;
	od:

	TransitiveGroup(M_upd,P)
end: