/******************************************************************

 Latest version
 Last modified: 01/12/2017

 Works in k=FiniteField(p^d), as well as in characteristic 0.
 If k is finite and there are not enough elements inside k
 to compute the partition, this algorithm extends the base field k.

 When char(k)=0, this algorithm does not reduce the equations of the system
 modulo a prime number, only modulo the maximal ideal of some valuation ring.


 To use this code, define the rational function. For instance:

    > k:=GF(3);
    > F<t>:=FunctionField(k);
    > f:=(t^10 + 2*t^8 + 2*t^4 + t^2)/(t^12 + 2*t^6 + 1);

 Then load this file:

    > load Decompose;

 Finally, compute a (single) complete decomposition:

    > Decompose(f,k: all_decomps:=false);

 The output, in this case, should be:

	[*
  	     t^2,
             t/(t^2 + 1),
	     (2*t^3 + t^2 + 2*t)/(t^3 + 2)
	*]

 Which means that 

	f = (2*t^3+t^2+2*t)/(t^3+1) o t/(t^2+1) o t^2

 over k=GF(3) and each component is indecomposable. To compute all
 non-equivalent complete decompositions, simply type:

     > Decompose(f,k);

 The output, in this case, should be

	[* [*
	    t^2,
	    t/(t^2 + 1),
	    (2*t^3 + t^2 + 2*t)/(t^3 + 2)
	*], [*
	    t/(t^2 + 2),
	    t^2,
	    t/(t^3 + 1)
	*], [*
	    t/(t^2 + 1),
	    t^3 + 2*t,
	    t^2
	*], [*
	    t/(t^2 + 1),
	    t^2 + t,
	    t^3 + t^2
	*], [*
	    t/(t^2 + 1),
	    t^2 + 2*t,
	    t^3 + t^2
	*], [*
	    t/(t^2 + 1),
	    t^2,
	    t^3 + t^2 + t
	*] *] 

********************************************************************/


forward normalize, LurothGen, right, normalize_decomp, leftfactor, AllSubf, mods, deg;


Decompose:=function(T,k: all_decomps:=true)
	// Input  : Rational Function T in k(t).
	// Output : All non-equivalent complete decompositions of T.

	// Prints some basic info
	info:=false;
	info1:=false;

	// Timings present the timings for several parts of the algorithm.
	timings:=false;

	// Show presents more basic information.
	show:=false;

	char:=Characteristic(k);	

	// Fields with less than small_fin_field elements are considered
	// small, so we extend the base field to compute the partitions.
	// See Remark 30 on Function Decompositions with Principal Subfields.

	small_fin_field:=46327;
	if char gt 0 then
		d:=#k;
		if d lt small_fin_field then 
			new_k:=true;
			exp:=Floor(Log(d, small_fin_field));
			K, h:=ext<k|exp>;
			if info1 then "New base field     :", K; end if;
		else
			K:=k;
			new_k:=false;
			h:=1;
		end if;
	else 
		new_k:=false; 
		K:=k; 
		h:=1;
	end if;

	F<t>:=FunctionField(k);
	P<x>:=PolynomialRing(F);
	P2:=PolynomialRing(k);
	H2x:=hom<F->P2|P2.1>;
	PP:=PolynomialRing(Integers());
	Hx:=hom<F->P|x>;

	t_ini:=Cputime();	

	fn:=F!Eltseq(Numerator(T));
	fd:=F!Eltseq(Denominator(T));
	fnx:=Hx(fn);
	fdx:=Hx(fd);

	f:=(F!fn)/(F!fd);
	d:=deg(f);

	MonicMinPol:=function(Gn, Gd, B)
		//  Input: Rational Function g in k(t)
		// output: The minimal polynomial of the extension k(t)/k(g(t))

		gnx:=Hx(Gn);
		gdx:=Hx(Gd);

		if Degree(gnx) gt Degree(gdx) then
			return (Gn/Gd), gnx-(Gn/Gd)*gdx, IdentityMatrix(K,2)*B;
		elif Degree(gdx) gt Degree(gnx) then
			return (Gd/Gn), gdx-(Gd/Gn)*gnx, Matrix(K,[[0,1],[1,0]])*B;
		else
			gn, gd, A:=normalize(Gn/Gd, k);
			return $$(gn, gd, A*B);
		end if;
	end function;

	// If char gt 0 and Derivative(mp) eq 0 then f can be decompose as
	// f = \tilde(f) o x^(p^k), where p is the characteristic of the field
	// and k is some positive integer.

	if char gt 0 then
		count:=0;
		while Derivative(fnx) eq 0 and Derivative(fdx) eq 0 do
			count+:=1;
			cn:=Eltseq(fnx); fnx:=P![cn[char*(i)+1] : i in [0..Degree(fnx)/char]];
			cd:=Eltseq(fdx); fdx:=P![cd[char*(i)+1] : i in [0..Degree(fdx)/char]];
			//f:=Evaluate(fnx,t)/Evaluate(fdx,t);
		end while;	
	
		if count gt 0 then 
			"Warning: The extension defined by this", (fd eq 1) select "polynomial" else 
			"rational function", "is not separable.";
			"This means that the", (fd eq 1) select "polynomial" else "rational function", "f is of the form g o x^",char^count," for some", (fd eq 1) select "polynomial" else "rational function", "g.";
			"We shall only compute the decompositions of g.";
		end if;	
	end if;

	f, Phi_f, A:=MonicMinPol(Evaluate(fnx,t), Evaluate(fdx,t), IdentityMatrix(k,2));
	fd:=Denominator(f); 
	fn:=Numerator(f);

	if show then 
		"Minimal Polynomial:  ", Phi_f;	
	end if;

	t_fact:=Cputime();
	Fact:=Factorization(Phi_f);
	Fact_orig:=[i[1]: i in Fact];

	if #Fact eq 2 then return [*f*]; end if;

	for i:=2 to #Fact do
		if Fact_orig[i] eq x-t then
			Fact_orig[i]:=Fact_orig[1];
			Fact_orig[1]:=x-t;
		end if;
	end for;

	Fact_orig1:=Fact_orig;

	F<t>:=FunctionField(K);
	P<x>:=PolynomialRing(F);
	P2:=PolynomialRing(K);
	H2x:=hom<F->P2|P2.1>;
	Hx:=hom<F->P|x>;

	if show then "Factorization ", Fact_orig; end if;

	r:=#Fact;

	if info1 then "Number of Factors  :", r; end if;

	if timings then
		"Time Factorization :", Cputime(t_fact);
	end if;


	if char ne 0 then
		p:=char; q:=#k;
	else
		p:=NextPrime(23); q:=p;
	end if;


	// Choosing a good k(t)-ideal:

	t_ideal:=Cputime();

	Px:=PolynomialRing(K);	
	count:=0;
	flag:=false;
	Deg:=1;
	poly_tries:=[];
	repeat
		count+:=1;
		if char gt 0 and count gt 30*Deg+ 5 then
			// increases the degree of the poly defining the ideal
			// count:=1;
			 Deg+:=1;
		 end if;

		ideal_poly:=Px!PP!RandomIrreduciblePolynomial(GF(p),Deg);

		if not ideal_poly in poly_tries then
			Append(~poly_tries, ideal_poly);
		else
			continue;
		end if;

		if not IsIrreducible(ideal_poly) then continue; end if;
		fd_red:=H2x(fd) mod P2!ideal_poly;
		fn_red:=H2x(fn) mod P2!ideal_poly;

		if fd_red eq 0 then continue; end if;

		if char ne 0 and fn_red/fd_red in k then continue; end if;

		if Degree(ideal_poly) ne 1 then
			FF<alpha>:=ExtensionField<K,alpha|ideal_poly>;
			FFx<X>:=PolynomialRing(FF);
			halpha:=hom<F->FF|alpha>;
			halpha2:=hom<FF->Px|Px.1>;
		else
			FF:=K;
			FFx<X>:=PolynomialRing(FF);
			halpha:=hom<F->FF|FF!(PP.1-ideal_poly)>;
			halpha2:=hom<FFx->FFx| Px.1>;
		end if;
		
		Phi_f_mod_I:=FFx![halpha(F!cc): cc in Coefficients(Phi_f)];;

		g:=GCD(Phi_f_mod_I, Derivative(Phi_f_mod_I));

		// it seems faster to compute the gcd of phi and phi' to
		// determine if phi is separable then using the built-in
		// function IsSeparable(phi);

		if Degree(g) ne 0 then flag:=false; else flag:=true; end if;

	until flag;

	assert(flag eq true);

	if timings then "Time k(t)-G. ideal :", Cputime(t_ideal); end if;
	if info then    "Degree ideal poly  :", Degree(ideal_poly); end if;

	// If the base field k is too small to compute the partitions, convert
	// the factorization of Phi_f over k(t) to K(t), where K is an extension of k.

	ChangeUniverse(~Fact_orig,P);


Poly:=function(Fact, computed_c)
	// input:  The factorization of Phi_f.
	//	   computed_c: set of elements c in K already used.
	// output: The polynomials \tilde{Q}_jc=P_jc(x)-h_jc(t)*L_c(x), j=1..r, used to compute
	//         the systems \tilde{S}_i, i=1..r;

	r:=#Fact;
	flag:=true;
	Cs:=computed_c;
	repeat
		flag:=false;
		try
			// Choose element c inside K:
			not_already_tried:=false;
			
			repeat
				if char eq 0 then
					c:=#Cs;
					if not c in Cs then not_already_tried:=true; Append(~Cs, c); end if;
				else
					c:=Random(K);
					if not c in Cs then not_already_tried:=true; Append(~Cs, c); end if;
				end if;
			until not_already_tried;

			if info then "c", c, count+1; end if;

			h_jc:=func<|[Evaluate(Derivative(Fact_orig[i]),c)/Evaluate(Fact_orig[i], c): i in [1..r]]>();
	
			// compute modulo good ideal
			h_jc_hat:=[halpha(F!h_jc[i]): i in [1..r]];
	
			// If the algorithm reached this point it means that the elements in h_jc_hat are in O_p
			flag:=true;
		catch e
			flag:=false;
			if info then "Bad random point chosen: some f_i(c) is not invertible mod I!"; end if;
		end try;
	until flag;

	assert(flag);

	HH:=hom<P2->FFx|FFx.1>;

	Nj:=[P2!Eltseq(Numerator(h_jc[i])): i in [1..r]];
	Dj:=[P2!Eltseq(Denominator(h_jc[i])): i in [1..r]];

	L_c:=LCM(Dj);
	P_jc:=[P2!(L_c/Dj[i])*Nj[i]: i in [1..r]];	

	H1:=[(FFx![FF!c: c in Coefficients(P_jc[i])]): i in [1..r]];
	H2:=FFx!(L_c);

	Q_jc:=[(H1[i]-H2*h_jc_hat[i]): i in [1..r]];

return [*Q_jc, Cs*];
end function;



System:=function(Q_jc, ind, F_ind_hat)
	// Input: index ind.
	// 	  Polynomial Q_jc used to construct the system S_ind.
	//	  The image F_ind_hat of F_ind modulo a good ideal P.
	// Output:A system S_ind which (probably) gives us the partition P_ind.

	R_jc:=[FFx!(Q_jc[j] mod F_ind_hat): j in [1..r]];
	
	// Take coefficients in x
	di:=Degree(F_ind_hat);
	if not di eq 1 then 
		Rest:=[[Coefficient(R_jc[j],i): j in [1..r]]: i in [0..di-1]];
	else
		Rest:=[R_jc];
	end if;

	Si:={ };

	for rest in Rest do
		Q_hat:=[rest[l]: l in [1..r]];

		if Degree(ideal_poly) ne 1 then
			coefsx:=[halpha2(l): l in Q_hat];
			max_deg,_:=Max([Degree(coefsx[l]): l in [1..#Q_hat]]);
			M:={[Coefficient(coefsx[l],s): l in [1..#Q_hat] ] : s in [0..max_deg]};
		else
			M:={[Q_hat[l]: l in [1..#Q_hat] ]};
		end if;
				
		if M eq {} then continue; end if;
		Si:=Si join M;
	end for;

return Si;
end function;



Check_sol:=function(Sol, ind)
	// Input: Basis of solutions of a system of squations.
	// output: The partitions P_ind of the principal subfields L_ind
	// 	   or false.

	F<t>:=FunctionField(k);
	P<x>:=PolynomialRing(F);
	P2:=PolynomialRing(k);
	H2x:=hom<F->P2|P2.1>;
	Hx:=hom<F->P|x>;

	if Degree(ideal_poly) ne 1 then
		FF<alpha>:=ExtensionField<k,alpha|ideal_poly>;
		FFx<X>:=PolynomialRing(FF);
		halpha:=hom<F->FF|alpha>;
		halpha2:=hom<FF->Px|Px.1>;
	else
		FF:=k;
		FFx<X>:=PolynomialRing(FF);
		halpha:=hom<F->FF|FF!(PP.1-ideal_poly)>;
		halpha2:=hom<FFx->FFx| Px.1>;
	end if;

	// This is a simple check to verify that the sum of the solutions
	// is the vector (1,1,...,1);

	for i:=1 to Nrows(Sol) do
		if #({Sol[i][j]: j in [1..Ncols(Sol)]} diff {0}) ne 1 then return false; end if;
	end for;

	for i:=1 to Ncols(Sol) do
		c:=0;
		for j:=1 to Nrows(Sol) do 
			a:=Sol[j,i];
			if char eq 0 then
				if not a in {0,1,-1} then return false; end if;
			else
				if not mods(a,p,q) in {0,1,-1} then return false; end if;
			end if;
			c+:=Abs(Integers()!mods(a,p,q));
		end for;
		if c ne 1 then
			return false; 
		end if;
	end for;

	Fi_hat:=FFx![halpha(ci): ci in Coefficients(Fact_orig1[ind])];
		
	// This check uses the factorization reduced mod p
	t_ver_modp:=Cputime();
	nr:=Nrows(Sol);
	for i:=1 to nr do
		S:={s : s in [1..Ncols(Sol)] | Abs(Integers()!mods(Sol[i,s],p,q)) eq 1};
		g:=P!&*[Fact_orig1[j]: j in S];
		g:= g/LeadingCoefficient(g);
		C:=Coefficients(g);
		for c in C do

			if deg(c) eq 0 then continue; end if;
			if 1 in S and not assigned(gener_n) then 
				// We get a Luroth generator of L_ind for free!
				if #S eq 1 then 
					gener_n:=t; gener_d:=1; 
				else
					gener_n:=Numerator(c); gener_d:= Denominator(c); 
				end if;
			end if;
			cn:=Numerator(c);
			cd:=Denominator(c);
			anx:=FFx!Eltseq(cn);
			ca:=halpha(c);
			adx:=FFx!Eltseq(cd);
			a:=(FFx!(anx-ca*adx))mod Fi_hat;
			if not a eq 0 then 
				if info then "Not correct partition - too refined(mod p check)"; end if;
				return false;
			end if;
		end for;
	end for;
	t_ver_modp:=Cputime(t_ver_modp);

	// Finally, compute the partition-vector.
	vi:=[0: i in [1..r]];
       	for k:=1 to Nrows(Sol) do
       	   	for l:=1 to r do
               		if Abs(Integers()!mods(Sol[k,l],p,q)) eq 1 then 
				ii:=l; break l; 	
			end if;
            	end for;
		for l:=1 to r do
                       	if Abs(Integers()!mods(Sol[k,l],p,q)) eq 1 then 
				vi[l]:=ii; 
			end if;
               	end for;
        end for;
A:=[*vi, gener_n/gener_d*];
return A;
end function;


MyPartitions:=function(Fact)
	// Input:  The factorization Fact of the minimal polynomial Phi_f
	//         of the extension k(t)/k(f(t)).
	// Output: The partitions P1,...,Pr of the principal subfields of this extension.

	t_total:=Cputime();

	Fi_hat:=[FFx![halpha(ci): ci in Coefficients(Fact[index])]: index in [1..r]];

	all_part:=false;
	partitions:=[*[*[0],0*]: i in [1..r]*];
	partitions[1]:=[*[i: i in [1..r]],t*];
	count:=0;
	part:={i: i in [2..r]};
	Sys:=[*{}: i in [1..r]*];

	t_poly_t:=0;
	t_sys:=0;
	t_sol:=0;
	t_check:=0;
	computed_c:=[];
	repeat
		if count gt d^10 then return []; end if;

		tt:=Cputime();
		A1:=Poly(Fact_orig, computed_c);

		poly:=A1[1]; computed_c:=A1[2];
		t_poly_t+:=Cputime(tt);

		count+:=1;

		for j in part do
			ind:=j;
			if info then "Computing Partition ", ind; end if;
			tt:=Cputime();
			Sys2:=System(poly, ind, Fi_hat[ind]);
			t_sys+:=Cputime(tt);
			

			tt:=Cputime();
			Sys[ind]:=Sys[ind] join Sys2;
			M:=Transpose(Matrix(K,SetToSequence(Sys[ind])));
			Sol:=NullspaceMatrix(M);
			t_sol+:=Cputime(tt);

			tt:=Cputime();
			s:=Check_sol(Sol, ind);

			t_check+:=Cputime(tt);

			if not s cmpeq false then
				partitions[ind][2]:=s[2];
				partitions[ind][1]:=s[1];
				Exclude(~part,j);
			end if;
		end for;
	until part eq {};

	if info1 then "dp, number of c's  :", Degree(ideal_poly), #computed_c; end if;
	if timings then	"Time Partitions    :", Cputime(t_total); end if;
	if timings then "Time Polynomials   :   ",  t_poly_t; end if;
	if timings then "Time System        :   ",  t_sys; end if;
	if timings then "Time to Solve Syst.:   ",  t_sol; end if;
	if timings then "Time to Check Sol. :   ",  t_check; end if;

return partitions;
end function;




	F<t>:=FunctionField(k);
	P<x>:=PolynomialRing(F);
	P2:=PolynomialRing(k);
	H2x:=hom<F->P2|P2.1>;
	PP:=PolynomialRing(Integers());
	Hx:=hom<F->P|x>;




////////////////////////////////// Computing Partitions ///////////////////////////////////////



	t_part:=Cputime();	
	Parts:=MyPartitions(Fact_orig1);

	if info then
		"Number of Partitions: ", #Parts;
	end if;

	if show then
		"Partitions: ", Parts;
	end if;




///////////////////////////////// Intersecting Subfields ////////////////////////////////////////


	t_inter:=Cputime();
	Subf:=AllSubf(Parts);

	if timings then
		"Time Intersections :", Cputime(t_inter);
	end if;

	if info1 then 
		"Number of Subfields:",#Subf;
	end if;

	if show then
		"Subfields: ", Subf;
	end if;


////////////////// For each subfield, find a Luroth generator /////////////////////////////////


	t_luroth:=Cputime();

	Subf_list:=[{i:i in [1..r]| s[1][i] eq 1} : s in Subf];

	PSubf_list:=[{{i:i in [1..r]} | i : i in [1..r] | P[1][i] eq 1} : P in Parts];

	for i:=1 to #Subf_list do
	
		if #Subf_list[i] eq 1 then Subf[i][2]:=t; continue i; end if;
		if #Subf_list[i] eq r then Subf[i][2]:=f; continue i; end if;

		s:=Position(PSubf_list, Subf_list[i]);
		if s ne 0 then
			Subf[i][2]:=Parts[s][2];
			 continue i; 
		end if;
		g:=&*[Fact_orig1[k]: k in Subf_list[i]];
		G:=LurothGen(g,k);
		if G ne [] then
			gn:=F!G[1]; gd:=F!G[2];
		else
			Error("wrong reconstruction");
		end if;
		Subf[i][2]:=gn/gd;
	end for;

	if timings then "Time Find L. Gen.  :", Cputime(t_luroth); end if;



////////////////////////////////// Computing Maximal Chains ///////////////////////////////////////	


	t_maxc:=Cputime();

	S:=[{{i:i in [1..r]} | i : i in [1..r] | P[1][i] eq 1} : P in Subf];
	S:={s : s in S};

	max_c:=function(v, S)
		// Input:   Set of sets S and v, a subset of S.
		// Output:  The maximal chains of sets in S containing v.

		T2:={};
		I:=#v;

		while I le r  do
			I+:=1;
			T2:=T2 join {s : s in S | #s eq I and v subset s and {t : t in T2 | t subset s} eq {}};
		end while;

		if T2 eq {} then return [[v]]; end if;

		chains:=[];
		for s in T2 do
			if v subset s and v ne s then	
				S := S diff {s};
				a:=$$(s, S);
				for c in a do Append(~chains, [v] cat c); end for;	
			end if;
		end for;
		return chains;
	end function;


	Chains:=[];
	T:={};

	for i:=2 to r-1 do
		T:=T join {s : s in S | #s eq i and {t : t in T | t subset s} eq {}};
	end for;
	de_norm:=function(g,A)
		return (F!A[1][1]*g+F!A[1][2])/(F!A[2][1]*g+F!A[2][2]);
	end function;



	// This piece of code returns only one decomposition of f. This is helpful
	// if one wants only a single decomposition but f has several decompositions.
	
	if not all_decomps then
		t_time:=Cputime();
		T:=[t : t in T];
		a:=max_c(T[1],S)[1];
		l:=#a;
		gen:=[];
		for j:=1 to l do
			s:=Position(Subf_list, a[j]);
			G:=Subf[s][2];
			Append(~gen, [F!Numerator(G), F!Denominator(G)]);
		end for;
		a:=gen;
		c1:=Evaluate(a[1][1],t)/Evaluate(a[1][2],t);
		dec:=[* c1 *];
		for j:=2 to #a do
			c1:=leftfactor(a[j], a[j-1],k);
			Append(~dec, c1);
		end for;

		dec:=normalize_decomp(dec,k);
		dec[#dec]:=de_norm(dec[#dec],A^-1);
		if timings then "Time left comp.    :",  Cputime(t_time); end if;
	return dec;
	end if;



	for v in T do
		a:=max_c(v, S);
		for s in a do Append(~Chains,s); end for;
	end for;

	if timings then
		"Time Maximal Chains:", Cputime(t_maxc);
	end if;

	if info then
		"Number of maximal chains: ", #Chains;
	end if;


///////////////////////// Writing Maximal chains with luroth generators  /////////////////////////////

	chains:=[];
	for c in Chains do
		l:=#c;
		gen:=[];
		for j:=1 to l do
			s:=Position(Subf_list, c[j]);
			G:=Subf[s][2];
			Append(~gen, [F!Numerator(G), F!Denominator(G)]);
		end for;
		Append(~chains, gen);
	end for;	


//////////////////////// Computing Left Components  ///////////////////////////////////////


	dec_chains:=[];
	t_dec_c:=Cputime();
	for c in chains do
		c1:=Evaluate(c[1][1],t)/Evaluate(c[1][2],t);
		dec:=[* c1 *];
		for j:=2 to #c do
			c1:=leftfactor(c[j], c[j-1],k);
			Append(~dec, c1);
		end for;
		Append(~dec_chains, dec);
	end for;

	if timings then
		"Time Left Comp.    :", Cputime(t_dec_c);
	end if;


//////////////////////// Normalizing the decomposition chains ///////////////////////////////////////

	t_norm:=Cputime();
	Norm_Decompositions:=[* *];

	for d in dec_chains do
		norm_dec_chain:=normalize_decomp(d,k);
		norm_dec_chain[#norm_dec_chain]:=de_norm(norm_dec_chain[#norm_dec_chain],A^-1);
		// The line above, if uncommented, "corrects" the output, so that it gives
		// the decomposition of the input rational function f and not normalize(f).
		Append(~Norm_Decompositions, norm_dec_chain);
	end for;

	if timings then
		"Time to Normalize  :", Cputime(t_norm);
	end if;

return Norm_Decompositions;
end function;






/////////////////////////////////////// Additional functions /////////////////////////////////////////////


mods:=function(a,p,q)
	if p eq 0 then return a; end if;
        if p eq q then 
		a:=Integers()!a;
	        b:=Integers()!(a mod q);
		if b le Floor(q/2) then
			return b;
	        else 
			return (b-q);
		end if;
	else
		if a in GF(p) then return $$(a,p,p); end if;
		return a;
	end if;	
end function;


deg:=function (g)
	dn:=Degree(Numerator(g));
	dd:=Degree(Denominator(g));
	if dn ge dd then 
		return dn; 
	else 
		return dd; 
	end if;
end function;



normalize :=function(f,k)

	// This function was implemented by [Mohamed Ayad, Peter Fleischmann]

	// Input : function f in K(x)\ K
	// Output: N:= the normal form of f in left GL_2 - orbit
	//       : M:= 2 x 2 matrix with N = M o f.

	local f_d,f_n,p_0,q_0,N,p_m,q_n,m,n,M,U;
	F<t>:=FieldOfFractions(PolynomialRing(k));
	P<x>:=PolynomialRing(F);

	if f eq 0 then return 0; end if;
	if Degree(f) eq 0 then return 1,1,1; end if;
	M := Matrix(k, 2, 2, [1,0,0, 1]);
	N:=f;
	p:=Numerator(N); q:=Denominator(N); m:=Degree(p); n:=Degree(q);
	p_0:=Coefficient(p,0); q_0:=Coefficient(q,0);

	// first we want p_0=0;
        if not p_0 eq 0 then // if p_0=0 skip this, so here p_0 <> 0

                if (q_0 eq 0) then
                        N:=Evaluate(1/t,N); M:=Matrix(k, 2, 2, [0,1,1,0]);
                        p:=Numerator(N); q:=Denominator(N);
                        m:=Degree(p); n:=Degree(q);
                        p_0:=Coefficient(p,0); q_0:=Coefficient(q,0);

                else // now p_0*q_0 is nonzero
                        u:=(-q_0/p_0)*t+1;
                        U := Matrix(k, 2, 2, [-q_0/p_0,1,0,1]);
                        N:=Evaluate(u,N); M:=U*M;
                        p:=Numerator(N); q:=Denominator(N);
                        m:=Degree(p); n:=Degree(q);
                end if;
        end if; // now p_0 = 0 and q_0 <>0

        if m eq n then
                p_m:= Coefficient(p,m); q_n:= Coefficient(q,m);
                u:= t/(-q_n/p_m*t + 1); U:=Matrix(k, 2, 2,[1,0,-q_n/p_m,1]);
                N:=Evaluate(u,N); M:=U*M; p:=Numerator(N); q:=Denominator(N);
                m:=Degree(p); n:=Degree(q);
        end if; // now p_0 = 0 and q_0 <>0;
                // (p,q)=1, (which is invariant under left GL_2 -action)
                // p and q have different degrees.

        p_m:= Coefficient(p,m); q_n:= Coefficient(q,n);
        u:= (q_n/p_m)*t; U:=Matrix(k, 2, 2,[q_n/p_m,0,0,1]); // make p,q monic
        N:=Evaluate(u,N); M:=U*M; p:=Numerator(N); q:=Denominator(N);
        m:=Degree(p); n:=Degree(q);

        if m lt n then // remove q_m
        q_n:=Coefficient(q,m); u:= t/(-q_n*t + 1);
        U:=Matrix(k, 2, 2,[1,0,-q_n,1]);
        N:=Evaluate(u,N); M:=U*M;
        end if;

return Numerator(N), Denominator(N), M;
end function;


LurothGen:=function(g,k)
	// Given a polynomial g \in K(t)[x], returns any coefficient not in K.

	char:=Characteristic(k);
	F<t>:=FieldOfFractions(PolynomialRing(k));
	C:=Coefficients(g);
	for c in C do
		if not c in k then
			//gn, gd, _:=normalize(F!c,k);;
			gn:=Numerator(c); gd:=Denominator(c);
			return [gn, gd];
		end if;
	end for;
end function;



right := function(f,M,k)

	// This function was implemented by [Mohamed Ayad, Peter Fleischmann]

	// right action of GL_2 on K(x)
	// Input:  f in FF, M matrix in C^2;
	// Output: f^M;
	F<t>:=FieldOfFractions(PolynomialRing(k));
	//f:=F!T;
	P<x>:=PolynomialRing(F);
	if Determinant(M) eq 0 then 
		printf "Matrix M is not invertible\n";
		return(0);
	end if;
	u:=(M[1,1]*t+M[1,2])/(M[2,1]*t+M[2,2]);
return(Evaluate(f,u));
end function;


normalize_decomp := function(h_list,k)

	// This function was implemented by [Mohamed Ayad, Peter Fleischmann]

	//normalizes a decomposition list
	//ie. all factors except possibly the "leftmost" (i.e. last in the list) 
	//are normalized
	//Input: h_list
	F<t>:=FieldOfFractions(PolynomialRing(k));
	//f:=F!T;
	P<x>:=PolynomialRing(F);

	M := Matrix(k, 2, 2, [1,0,0, 1]); n_list:=[**];
	length:=#h_list; i:=1;
	while i lt length do
	    zn, zd,M:=normalize(right(h_list[i],M^-1,k),k);
	    Append(~n_list,Evaluate(zn,t)/Evaluate(zd,t));
	    i:=i+1;
	end while;
	Append(~n_list,right(h_list[length],M^-1,k));
return n_list;
end function;


leftfactor := function(f,h,k)

	// This function was implemented by [Mohamed Ayad, Peter Fleischmann]

	// Input: f=p/q and h=r/s in FF, where h is right-factor of f.
	// Output: g in FF with f = g o h,
	local g;
	F<t>:=FieldOfFractions(PolynomialRing(k));
	//f:=F!T;
	P<x>:=PolynomialRing(F);
	K<w>:=PolynomialRing(k);
	g:=0;
	h_n:= h[1]; h_d:=h[2];
	f_n:=f[1]; f_d:=f[2];
	ff_n:=hom< K -> P | P.1 >(f_n); ff_d:=hom< K -> P | P.1 >(f_d);
	//"ff_m", ff_n;
	hh:= hom< K -> P | P.1 >(h_n) -t*hom< K -> P | P.1 >(h_d);
	//"hh", hh;
	//hh:= K2!h_n - K2.1*K2!h_d;
	B:=ff_n mod hh; D:=ff_d mod hh;
	//"B", B; "D", D;
	g:=(LeadingCoefficient(B)/LeadingCoefficient(D));
	//"right factor", g;
return g;
end function;



// The algorithm AllSubfields form the paper HKN 2013:
// Fields are represented by r-partition-vetors
// and intersections are computed using the new method ( vH, Sz 2015).


AllSubf:=function(Part)
	// Input  : Set of partitions Part.
	// Output : All the intersections of the partitions in Part.

	r:=#Part[1][1];
	K:=[*[i : i in [1..r]], 0*];
	ListSubfields:=[K];
	e:=Vector(Integers(),[(Part[i][1] eq K[1]) select 1 else 0: i in [1..r]]);

	// The function P^\infty
	PInf:=function(a,P)
		if P[a] eq a then 
			return a; 
		else 
			return $$(P[a],P); 
		end if;
	end function;

	// The intersection of two partitions:
	Intersect:=function(v,w,r)
		P:=w[1];
		for a:=1 to r do
			b:=v[1][a];
			A:=PInf(a,P);
			B:=PInf(b,P);
			if B gt A then 
				P[B]:=A;
			elif A gt B then 
				P[A]:=B; 
			end if;
		end for;
		// Normalizing the partition:
		for i:=1 to r do P[i]:=P[P[i]]; end for;
	if P eq v[1] then return v; else
	return [*P,0*]; end if;
	end function;

	NextSubfield:=procedure(~ListSubfields, M, e, s)
		for i:=s+1 to r do
			if e[i] eq 0 then
				N:=Intersect(Part[i],M,r);
				ee:=Vector(Integers(), [N[1][j] eq 1 select 1 else 0 : j in [1..r]]);
				for k:=1 to i-1 do
					if ee[k] gt e[k] then
						continue i;				
					end if;
				end for;
				Append(~ListSubfields, N);
				$$(~ListSubfields, N, ee, i);
			end if;
		end for;
	end procedure;

	NextSubfield(~ListSubfields, K, e, 0);	

return ListSubfields;
end function;