# File: http://www.math.fsu.edu/~hoeij/files/abel_H72
#
# The following code requires Maple V release 5. It calculates the polynomial Q(Z^3)
# in section 2 in my paper "The Minimum Polynomial of an Algebraic Solution of
# Abel's problem".
#
# Although the implementation is not yet complete, it could be used to do
# other examples as long as you make sure that the variable d has the
# correct value, the variable acc is large enough, and subs(x=0,P) is
# irreducible in Q[Z].
#
# Mark van Hoeij, January 2000.
#
# Feb 29: bug fix, the code only worked when P was monic.


# F must be monic in Z.
# p=point, specified by either infinity or an irreducible polynomial.
bound_point:=proc(F,Z,x,p)
	local n,i,G;
	options remember;

	n:=degree(F,Z);
	if n<1 then
		ERROR()
	elif p=infinity then
		G:=collect(subs(x=1/x,Z=Z*x,F*x^n),Z,Normalizer);
		G:=traperror(subs(x=0,G));
		if G=lasterror then
			if G=`division by zero` then
				infinity
			else
				ERROR(G)
			fi
		else
			G:=select(type,map(i -> -i[1],roots(G,Z)),rational);
			min(infinity,op(G))
		fi
	elif lcoeff(F,Z)<>1 then
		procname(collect(F/lcoeff(F,Z),Z,Normalizer),Z,x,p)
	else
		G:=collect(subs(Z=Z*diff(p,x)/p,F*p^n),Z,Normalizer);
		G:=evala(Expand(rem(numer(Normalizer(G)),p,x)));
		if G=0 then
			infinity
		else
			G:=select(type,map(i -> i[1],roots(G,Z)),rational);
			min(infinity,op(G))
		fi
	fi
end:

bound_i:=proc(F,Z,x,i)
	local b,D,j,P,N;
	options remember;
	if lcoeff(F,Z)<>1 then
		procname(collect(F/lcoeff(F,Z),Z,Normalizer),Z,x,i)
	else
		N:=bound_point(F,Z,x,infinity);
		if N=infinity then RETURN(N) fi;
		N:=-ceil(i*N);
		D:=map(j -> j[1],factors(denom(Normalizer(F)),indets([args],RootOf))[2]);
		P:=1;
		for j in D do
			b:=bound_point(F,Z,x,j);
			if b=infinity then RETURN(b) fi;
			b:=ceil(i*b);
			P:=P*j^b;
			N:=N-degree(j,x)*b
		od;
		[P,N]
	fi
end:


P:=
Z^9+30*x/(3*x^2+1)*Z^8+1/4*(1599*x^2+1)/(3*x^2+1)^2*Z^7+1/72*(223587*x^3+15*x^
2+441*x+5)/(3*x^2+1)^3*Z^6+1/1152*(17855019*x^4+4800*x^3+73962*x^2+1600*x-21)/
(3*x^2+1)^4*Z^5+1/3456*(178134957*x^5+119907*x^4+1290114*x^3+40062*x^2-1071*x+
31)/(3*x^2+1)^5*Z^4+1/248832*(28420622685*x^6+38312352*x^5+323581833*x^4+
12864960*x^3-524073*x^2+31392*x+67)/(3*x^2+1)^6*Z^3+1/497664*(80931106503*x^7+
191139507*x^6+1350862263*x^5+64704555*x^4-3559875*x^3+330129*x^2+1349*x-111)/(
3*x^2+1)^7*Z^2+1/47775744*(6449857313913*x^8+24395769216*x^7+150200784540*x^6+
8354150784*x^5-580059306*x^4+73858176*x^3+437148*x^2-72576*x+505)/(3*x^2+1)^8*
Z+1/644972544*(32112055562886*x^9+182329225533*x^8+1005364927404*x^7+
63398867400*x^6-5314028328*x^5+869144634*x^4+6682932*x^3-1669824*x^2+23490*x-\
127)/(3*x^2+1)^9;

# The following numbers depend on P.
d:=3;
acc:=8;

n:=degree(P,Z);
p:=algcurves[puiseux](P,x=0,Z,acc-1):
if nops(p)<>1 then
  ERROR("this preliminary implementation supposes that all Puiseux expansions at x=0 are conjugated")
fi;
r:=op(indets(p,RootOf));
p:=convert(series(p[1],x=0,acc-1),polynom):
p:=series(exp(int(d*p,x)),x=0,acc):
p:=evala(Expand( convert(p,polynom) )):
pps:=[]:
for i from 1 to n do
	B:=bound_i(P,Z,x,i*d);
	pp[i]:=evala(Expand(convert(series(`if`(i=1,p,p*pp[i-1]),x=0,acc),polynom)));
	pi:=convert(series(pp[i]/B[1],x=0,acc),polynom);
	C[i]:=add(c[i,j]*r^j,j=0..n-1);
	if nops(indets(C[1]))=1 then
		C[i]:=evala(Expand(C[1]^i))
	fi;
	pi:=evala(Trace(C[i]*pi));
	pps:=[op(pps),pi*B[1]];
	if nops(indets(C[1]))=1 then
		next
	fi;
	S:={seq(coeff(pi,x,j),j=B[2]+1..acc-1)};
	S:=solve(S,{seq(c[i,j],j=0..n-1)}); # solve linear equations.
	pps:=subs(S,pps);
	C[i]:=subs(S,C[i]);
od:

F:=Z^n+add(E.i*Z^(n-i),i=1..n);
R:=RootOf(F,Z):
F:=sort(subs(solve({seq(P.i - evala(Trace(R^i)),i=1..n)},
  {seq(E.i,i=1..n)}),F),Z);
Q:=collect(subs(Z=Z^d,seq(P.i=pps[i],i=1..n),F),Z,factor);

# Now substitute a suitable value for the remaining constant c[..,..]
# to obtain:
QQ:=
Z^27+168/(3*x^2+1)^5*Z^18+405*(x+1)/(3*x^2+1)^7*Z^15-636/(3*x^2+1)^10*Z^9+324*
(x+1)/(3*x^2+1)^12*Z^6-243/4*(x+1)^2/(3*x^2+1)^14*Z^3+8/(3*x^2+1)^15;

# The Galois group of QQ over C(x) has 3*72=216 elements, which can
# be verified in Maple V Release 6 (which will soon be available)
# with the following commands:
#
# G:=algcurves[monodromy](numer(normal(QQ)),x,Z,group);
# group[grouporder](G);