/****************************************************************** 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:=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; 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:=FunctionField(k); P:=PolynomialRing(F); P2:=PolynomialRing(k); H2x:=homP2|P2.1>; PP:=PolynomialRing(Integers()); Hx:=homP|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:=FunctionField(K); P:=PolynomialRing(F); P2:=PolynomialRing(K); H2x:=homP2|P2.1>; Hx:=homP|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:=ExtensionField; FFx:=PolynomialRing(FF); halpha:=homFF|alpha>; halpha2:=homPx|Px.1>; else FF:=K; FFx:=PolynomialRing(FF); halpha:=homFF|FF!(PP.1-ideal_poly)>; halpha2:=homFFx| 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:=homFFx|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:=FunctionField(k); P:=PolynomialRing(F); P2:=PolynomialRing(k); H2x:=homP2|P2.1>; Hx:=homP|x>; if Degree(ideal_poly) ne 1 then FF:=ExtensionField; FFx:=PolynomialRing(FF); halpha:=homFF|alpha>; halpha2:=homPx|Px.1>; else FF:=k; FFx:=PolynomialRing(FF); halpha:=homFF|FF!(PP.1-ideal_poly)>; halpha2:=homFFx| 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:=FunctionField(k); P:=PolynomialRing(F); P2:=PolynomialRing(k); H2x:=homP2|P2.1>; PP:=PolynomialRing(Integers()); Hx:=homP|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:=FieldOfFractions(PolynomialRing(k)); P:=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:=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:=FieldOfFractions(PolynomialRing(k)); //f:=F!T; P:=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:=FieldOfFractions(PolynomialRing(k)); //f:=F!T; P:=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:=FieldOfFractions(PolynomialRing(k)); //f:=F!T; P:=PolynomialRing(F); K:=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;