iso1p := proc(f,g,x) # input two polynomials f(x) and g(x), and this program outputs all the isomorphisms from Q[x]/(g) -> Q[x]/(f). # This happens by working with one prime at a time up to a maximum of k primes. # In this implementation, k is the number of primes and 5 is where to start picking the primes. # This implementation uses on;y one trace at a time, which should be faster when there are many isomorphisms. local p, L, i, j, l, count, A, M, n, m, b, e, B, C, q, h, W, Basis, prime, counter, Base, roots, iso, num, X, time1, time2, MaxP, k, w, fprime, isos, root, fp, gn, time3, galois, pp, D, stoped; #options trace; time1:= time(); lprint("Pre-Processing program starting at:", time1); # Find the bound for which vectors we can throw out. #Step 1 counter:=0; count:=0; b:= bound(f,g,x); `mod` := mods; gn:= lcoeff(g,x); lprint("the bound on the square length is:", b); n := degree(f); p[1] := 4; e := n+1; q:=0; # Step 1(e), note: they will be called Basis[i] when they have been reduced mod p^a, this is to make it faster # Make new basis 1/f', # alpha/ f', alpha^2/f', ... fprime:= subs(RootOf(f)=x, evala(1/subs(x=RootOf(f),diff(f,x)))); for i from 1 to n do Base[i] := [rem( expand(fprime* x^(i-1)) , f,x)]; od; #MaxP := round(15+sqrt(n)); #k:= MaxP; # Step 4 #for q from 1 to k do do q:= q+1; # Step 4(b) #Finds the suitable primes and splits the factors into lists by degree prime[q], p[q], roots[q], M[q] := FindSuitablePrime(f,g,x,p[q],b, e); p[q+1]:=prime[q]; # If the fields appear to be Galois we want to move to the second method if M[q] = 0 then C := Matrix(n+1,n+1,shape=identity); galois:=1; break; fi; lprint("the prime being used is:", prime[q]); lprint("p^a is:", p[q]); # Reduce the coefficients (mod p^a) of Basis to remove fractions for i from 1 to n do Basis[i]:= Base[i] mod p[q]; od; # Step 4(c) #Converts all of the lists of factors into lists with the 1..n-1 powers of the #sum of the roots # Repeating for the correct variable name m := nops(M[q]); for j from 1 to m do M[q][j] := Newton(M[q][j], n, p[q]); od; # Now we construct the matrix to run LLL on # Step 4(d) for j from 1 to m-1 do if q = 1 then A := Matrix(n+1+1, n+1+1); # Constructing matrix C (this is the identity) for i from 1 to n+1 do A[i,i] := 1; od; #count := 1; # Constructing marix T (traces under Basis_i) #for j from 1 to nops(M[q]) do for l from 1 to n do # Here we weight the traces of \alpha with the coefficients of 1/f'(\alpha) W[q][j][1][l] := add(coeff(op(Basis[l]),x,i-1)*M[q][j][1][i], i=1..n) mod p[q]; if _EnvMethod = 1 then A[ l , n + 1+1] := (10^10)*W[q][j][1][l] mod p[q]; else A[ l , n + 1+1] := W[q][j][1][l] mod p[q]; fi; od; # Constructing matrix P A[(n + 1+1), (n + 1+1)] := p[q]; if _EnvMethod = 1 then A[n+1,n+1+1] := (10^10)*gn*M[q][j][2][2] mod p[q]; else A[n+1,n+1+1] := gn*M[q][j][2][2] mod p[q]; fi; # count := count +1; #od; else #count := 1; A:= Matrix(e+1, n+1+1); # Constructing matrix C (coefficients) A[1..e,1..n+1]:=C; #for j from 1 to nops(M[q]) do # Here we weight the traces of \alpha with the coefficients of 1/f'(\alpha) # Constructing marix T (traces under Basis_i) for l from 1 to n do W[q][j][1][l] := add(coeff(op(Basis[l]),x,i-1)*M[q][j][1][i], i=1..n) mod p[q]; od; # Constructing marix CT (traces under Basis_i weighted with coefficients matrix) for l from 1 to e do A[l, n + 1 + 1] := ( add(A[l, i]*W[q][j][1][i], i=1..n) + A[l,n+1]*gn*M[q][j][2][2] ) mod p[q]; od; # Constructing matrix P A[(e + 1), (n + 1+1)] := p[q]; #count := count +1; #od; fi; #lprint("this is A", A); # If the traces are already zero modulo that prime then we want to go to the next prime # If for too many primes this is zero then we want to quit # Note, all the traces being zero is given by: [seq(seq(A[i,j],j=n+2..n+m+1),i=1..e)] mod prime[q] =[seq(0, i=1..e*m)] # Step 4(e) if [seq(A[i,n+2],i=1..e)] mod p[q] =[seq(0, i=1..e)] then C:= Matrix(e,n+1); C:= A[1..e,1..n+1]; count := count+1; lprint("This prime has zero traces"); if count < 10 then next; else lprint("Pre-Processing has now finished"); stoped :=1; k:=q; break; fi; else counter:= counter+1; fi; #if _EnvMethod = 1 then # This is to see whether putting the p^a rows first in LLL makes it faster # This appears to be slower #lprint("Inputting the matrix in LLL with the p^a rows first"); # # D := Matrix(e, n+m+1); # # D[1..m, 1..n+m+1] := A[e-m+1..e,1..n+m+1]; # # D[m+1..e,1..n+m+1] := A[1..e-m,1..n+m+1]; # # A[1..e,1..n+m+1] := D[1..e,1..n+m+1]; # # D:= 'D'; # #fi; # This version below is LLL with the vectors above the bound b, being thrown out. # Step 4(f) lprint("this is A:",A); lprint("calling LLL at time:", time()); L:= LLLInteger(A,'cut_away',b); lprint("number of vectors remaining after LLL:", nops(L), time()); # Here we split the result from LLL into two parts < C B> were C are the coefficients and B is the result # of the traces # Step 4(g) e:= nops(L); if e = 0 then lprint(" No isomorphism here at:", time()); return []; fi; B:= 'B'; C:= 'C': for i from 1 to e do B[i] := [seq(L[i][j], j=n+2..n+1+1)]; C[i] := [seq(L[i][j], j=1..n+1)]; od; B:= Matrix([seq(B[i],i=1..e)]); C:= Matrix([seq(C[i],i=1..e)]); # If the traces are not zero then we continue to multiply them by 10^n and run LLL again # Step 4(h) while not(verify(B, Matrix(e,1), 'Matrix')) do A:= Matrix(e, n+1+1); A[1..e, n+2..n+1+1] := (10^n)*B; A[1..e,1..n+1] := C; lprint("calling LLL after multiplying by 10^n at time:", time()); L:= LLLInteger(A,'cut_away',b); lprint("done with LLL 10^n at time:", time()); lprint("number of vectors remaining:", nops(L)); e:= nops(L); # Step 4(i) if e = 0 then time2:= time(); lprint("There are no isomorphisms"); lprint("program done early at:", time2); lprint("program took a total of:", (time2-time1)); lprint("number of vectors added:", counter); return e, L; fi; B:='B'; C:='C'; for i from 1 to e do B[i] := [seq(L[i][j], j=n+2..n+1+1)]; C[i] := [seq(L[i][j], j=1..n+1)]; od; B:= Matrix([seq(B[i],i=1..e)]); C:= Matrix([seq(C[i],i=1..e)]); od; if e = 1 then time3:= time(); # Test to make sure the isomorphism is really an isomorphism h:= subs( x=RootOf(f), add( (C[1,i]/(-C[1,n+1])) * x^(i-1) ,i=1..n) ); fprime:=subs(x=RootOf(f),fprime); iso:= evala( Expand(h * fprime*gn)); fp:= diff(f,x); # if evala( Expand( subs(x=iso,g) )) =0 then if VerifyIso(f,g,iso,n,x) = 1 then time2:= time(); lprint("There is one Verified isomorphism"); lprint("program done early at:", time2); lprint("program took a total of:", (time2-time1)); lprint("number of vectors added:", counter); lprint("Time taken to verify isomorphism:", (time2-time3)); return {subs(RootOf(f)=x, h/fp)}, lprint("where x=", RootOf(f)); else time2:= time(); lprint("There are no isomorphisms"); lprint("program done early at:", time2); lprint("program took a total of:", (time2-time1)); lprint("number of vectors added:", counter); return 0, {}; fi; fi; od; if stoped =1 then break; fi; od; lprint("going into recursive method at:", time()); lprint("This implementation assumes that there are linear roots"); # Here we will find the isomorphism by checking all possibilities, only if none is found above, this is slow isos := {}; # Steps 2 and 3 roots:= [seq(roots[i], i=1..k)]; # If there are no linear roots then we need to get some if [seq(roots[i],i=1..k)] = [seq(0,i=1..k)] then prime[1], p[1], roots[1], M[1] := FindSuitablePrime(f,g,x,prime[1]-1, b, e,'need_linear'); fi; # Step 4 for q from 1 to k do if roots[q] = 0 or roots[q][2] =[] then next else # Step 4(a) # We try all the combinations and see which ones work and if this doesn't give a definitive answer we Hensel Lift more utill it give a definitive answer for j from 1 to nops(roots[q][2]) do #prime[q], p[q], roots[q], M[q] := FindSuitablePrime(f,g,x,prime[q]-1, b, e, 'galois'); num, X[q][j]:= BruteIso(b, n, C, p[q], [seq(Base[i], i=1..n)], roots[q][2][j], roots[q][3][1], g); # Step 4(b), note: this comes fist because it must be programmed that way while num >1 do #b:= 5*b; e:= 3*e; prime[q], pp[q], root[q], M[q] := FindSuitablePrime(f,g,x,prime[q]-1, b, e, 'galois'); num, X[q][j]:= BruteIso(b, n, X[q][j], pp[q] , [seq(Base[i], i=1..n)], root[q][2][j], root[q][3][1], g); od; e:= n+1; if num = 1 then # Test to make sure the isomorphism is really an isomorphism h:= subs( x=RootOf(f), add( (X[q][j][1,i]/(-X[q][j][1,n+1])) * x^(i-1) ,i=1..n) ); fprime:=subs(x=RootOf(f),fprime); iso:= evala( Expand(h * fprime*gn)); fp:= diff(f,x); if VerifyIso(f,g,iso,n,x) = 1 then lprint("This is another Verified isomorphism"); isos:= isos union {subs(RootOf(f)=x, h/fp)}; else lprint("No verified isomorphism here", time()) fi; fi; od; # Output time2:= time(); lprint("program done at:", time2); lprint("program took a total of:", (time2-time1)); lprint("number of vectors added:", counter); lprint("Number of isomorphisms:", nops(isos)); return isos, lprint("where x=", RootOf(f)); fi; od; end: VerifyIso:= proc(f, g, iso, n, x) local B, i, fprime; #options trace; fprime:=subs(x=RootOf(f), diff(f,x)); for i from 1 to n+1 do B[i] := expandB(fprime, iso, i-1); if denom(B[i]) <> 1 then lprint("This is not an isomorphism, because it is not an Alg Int."); return 0; fi; od; # evala(Expand(add( coeff(g,x,i-1) * B[i] , i=1..n+1))); # evala(Expand( subs(x=iso, g))); if evala(Expand(add( coeff(g,x,i-1) * B[i]/ lcoeff(g,x)^(2*i-2), i=1..n+1))) = 0 then lprint("There is one Verified isomorphism"); return 1 else lprint("This is not an isomorphism"); return 0; fi; end: expandB := proc(fprime,A,i) option remember; if i=0 then return fprime; elif i=1 then return evala(Expand(fprime*A)); else evala(Expand(A* expandB(fprime,A,i-1) )); fi; end: Newton := proc(M, deg, p) local i, n, m, j, k, l, N; #options trace; # Input a list of the form [degree of factors, [factors of f],[factors of g]] and degree of min poly # Output a list of the form [[seq of ith power of roots of f of given degree],[seq of 0,1st powers of roots for g]] # also reduces mod p^a to have smaller elements # This is used to find the traces [[seq( add(PowSum(j,x,i,p), j = M[2]) mod p ,i=0..deg-1)], [seq( add(PowSum(j,x,i,p), j = M[3]) mod p ,i=0..1)]] end: PowSum := proc(f,x,i, p) # Input: f in K[x], i in N. # Output: the sum of the i'th powers of the roots of f. modulo p^a #options trace; option remember: local j,n; n := degree(f,x); if i=0 then n mod p elif i=1 then -coeff(f,x,n-1)/lcoeff(f,x) mod p else add( -coeff(f,x,n-j)*PowSum(f,x,i-j, p)/lcoeff(f,x) ,j=1..i-1) - i*coeff(f,x,n-i)/lcoeff(f,x) mod p fi; end: bound := proc(f, g, x) local n,G,F,i; #options trace; n:= degree(f); #Step 2 G := n*rootbound(g,x)*lcoeff(g,x); F:= add(i^2,i=coeffs(f)); #Step 3 (n^4)*(G^2)*F; end: FindSuitablePrime := proc(F, G, x, bp, bound, e) # Given polynomials F, G; Find a suitable prime > bp # Then lifts the factorization to the desired height so that p^a > bound^n # Output: p,p^a, linear roots, [L1, L2, ... ] # where each Li is of the form # [deg, [factors of F], [factors of G]] option remember; #options trace; local counter, p, F1, G1, F2, G2, i, degs, deg,n, a, s, linear, numf, numg, dummy1, dummy2, linf, ling, p_power, Fn, Gn, marked, j; # Step 1 marked:=0; counter := 0; p := bp; numf :=0; numg := 0; Fn := lcoeff(F,x); Gn := lcoeff(G,x); do # Step 2(a)(b)(c) p := nextprime(p); # Avoid p | discrim(F,x) if has(Gcd(diff(F,x),F) mod p, x) then next fi; if has(Gcd(diff(G,x),G) mod p, x) then next fi; if gcd(p, Fn) <> 1 then next fi; if gcd(p, Gn) <> 1 then next fi; # Step 2(d) F1:= expand(F/Fn) mod p; F1 := DistDeg(F1 ,x) mod p; #F1[1][1] := expand(Fn*F1[1][1]) mod p; # if we are running this because we need linear roots if nargs>6 and args[7]='need_linear' then if F1[1][2] = 1 then lprint("Found Linear Roots"); else next; fi; fi; if nops(F1)=1 then if nargs <7 then # Step 2(d) # DistDeg factorization is trivial, increase p: counter := counter + 1; if counter > 25 and F1[1][2] = 1 then lprint("They appear to be Galois"); marked:=1; else next; fi; else marked:=1; fi; fi; G1 := expand(G/Gn) mod p; G1 := DistDeg( G1 ,x) mod p; #G1[1][1] := expand(Gn*G1[1][1]) mod p; #Step 2(f) degs := sort([seq(degree(i[1],x), i= F1)]); if degs <> sort([seq(degree(i[1],x), i= G1)]) then error "not isomorphic", p fi; # note: will be repeated for the correct variable name F2[0] := nops(F1); G2[0] := nops(G1); # Find the places in the factorization where the linear roots are for i from 1 to F2[0] do if F1[i][2] = 1 then linf := Fact(F1[i][1],x,p); lprint( "The linear factors", linf); numf:= nops(linf[2]); F1 := [seq([linf[2][j][1],1], j=1..numf), op({op(F1)} minus {[F1[i][1],1]}) ]; fi; if G1[i][2] =1 then ling := Fact(G1[i][1],x,p); numg:= nops(ling[2]); G1 := [seq([ling[2][j][1],1], j=1..numg), op({op(G1)} minus {[G1[i][1],1]}) ]; fi; od; F2[0] := nops(F1); G2[0] := nops(G1); for i from 1 to F2[0] do # Convert the factors to modp1-notation, see: ?modp1 F2[i] := modp1('ConvertIn'(F1[i,1],x),p); od; for i from 1 to G2[0] do # Convert the factors to modp1-notation, see: ?modp1 G2[i] := modp1('ConvertIn'(G1[i,1],x),p); od; # This is the code that Maple uses for Hensel lifting, because of Quadratic Hensel Lifting it is lifted between p^a and p^2a: # Step 2(h) #lprint( "a is :", ceil(evalf(bound^(e/(2))))); #ceil( (bound^(e/10)) * (2^(1*e/4)) ); #ceil(evalf(bound^(e/(2)))) `factor/lift`(F, F2, x, p, dummy1, p_power, dummy2, `lift to`, ceil(evalf( (bound^(e/10)) * (2^(1*e/4)) ) )); dummy1 := 'dummy1'; dummy2:= 'dummy2'; p_power := 'p_power'; #seq(F2[i], i=0..nops(F1)); `factor/lift`(G, G2, x, p, dummy1, p_power, dummy2, `lift to`, ceil(evalf( (bound^(e/10)) * (2^(1*e/4)) ) )); # Step 2(g), make sure it is Hensel Lifted to the correct level (i.e. correct for quadratic Hensel Lifting) a := min(ceil(evalf( (bound^(e/10)) * (2^(1*e/4)) )), ln(p_power)/ln(p)); F2[1] := modp1('ConvertOut'(F2[1],x),p_power); F2[1] := expand(Fn*F2[1]); F2[1] := modp1('ConvertIn'(F2[1],x),p_power); F2:= [seq(modp1('ConvertOut'(F2[i],x),p_power) mod p^a, i=1..nops(F1))]; G2[1] := modp1('ConvertOut'(G2[1],x),p_power); G2[1] := expand(Gn*G2[1]); G2[1] := modp1('ConvertIn'(G2[1],x),p_power); G2:= [seq(modp1('ConvertOut'(G2[i],x),p_power) mod p^a, i=1..nops(G1))]; # Step 2(i) # Find the linear roots after Hensel Lifting if numf <> 0 then linf:= [seq(-coeff(F2[i],x,0)/lcoeff(F2[i],x), i=1..numf)]; ling := [seq(-coeff(G2[i],x,0)/lcoeff(G2[i],x), i=1..numg)]; linear:= [p^a, linf, ling]; lprint("There are linear roots", linear); else linear:=0; fi; if marked = 0 then # A list whose elements look like this: # p, p^a, linear , [deg, [factors of F], [factors of G]] return p, p^a, linear, [seq( [deg, [seq(`if`(degree(i,x)=deg, i, NULL), i=F2)], [seq(`if`(degree(i,x)=deg, i, NULL), i=G2)]], deg = degs)] elif marked = 1 then return p, p^a, linear, 0 fi; od; end: Fact := proc(f,x,p) # This is exactly maple's factoring, but with remember so it gives the factors in the same order each time. local L; option remember; L:= Factors(f) mod p; return L end: LLLInteger := proc(AA) local i,j,n,m,A,B,l,t,dd,k,r,d,c,d0,d1,d2,q,flag,kmax, y; t := AA; n,m := op(1,t); t := convert(t,listlist); for i to n do A[i] := t[i] od: B := Matrix(n,n); d := Array(0..n,[1,add(A[1][l]^2,l=1..m)]); userinfo(1,'LLL',`beginning reduction`); kmax := 1; for y in [3/10, 1/2, 3/4] do # y is the number in the LLL paper k := 2; flag := true; while kkmax then for i to k do t := add(A[i][l]*A[k][l],l=1..m); for l to i-1 do t := iquo(t*d[l]-B[l,i]*B[l,k],d[l-1]) od; if i=k then d[k] := t else B[i,k] := t fi; od; if d[k]=0 then userinfo(1,'LLL',`removing a dependent vector`); A[k] := A[n]; n := n-1; next fi; kmax := kmax+1; userinfo(2,'LLL',`reached column ` || kmax); fi; # Reduce column k. for i from k-1 by -1 to 1 while flag do r := mods(B[i,k],d[i]); q := iquo(B[i,k]-r,d[i]); if q=0 then next end if; A[k] := A[k]-q*A[i]; for j to i-1 do B[j,k] := B[j,k]-q*B[j,i] od; B[i,k] := r; od; l := k-1; d0 := d[k]; d1 := d[l]; d2 := d[k-2]; c := B[l,k]; dd := iquo(d0*d2+c^2,d1); if denom(y)*dd < numer(y)*d1 then userinfo(3,'LLL',`switching at row ` || k); A[k],A[l] := A[l],A[k]; for i to k-2 do B[i,k],B[i,l] := B[i,l],B[i,k] od; for j from k+1 to kmax do t := B[l,j]; B[l,j] := iquo(B[k,j]*d2+c*t,d1); B[k,j] := iquo(t*dd-c*B[l,j],d2) od; d[l] := dd; if k>2 then k := k-1; flag := false fi; else k := k+1; flag := true; fi; od; # The 'cut_away' option will be used to speed up `factor/knapsack` while nargs>2 and args[2]='cut_away' and n>0 and d[n] > args[3]*d[n-1] do n := n-1; kmax := n od; od; [seq(A[i],i=1..n)] end: BruteIso := proc(b, n, D, p, Base, rootf, rootg, g) # this is what we will use to construct the martix to run LLL on for the brute force part # Output either 0, [] or 1, M or e, C local l, i, W, A, M, B, C, j, e, Basis, L; #options trace; C:= D; e := linalg[rowdim](C); # Reduce the coefficients (mod p^a) of Basis to remove fractions for i from 1 to n do Basis[i]:= [op(Base[i]) mod p]; od; A:= Matrix(e+1, n+1+1); A[1..e,1..n+1]:=C; for l from 1 to n do # Here we weight the traces of \alpha with the coefficients of 1/f'(\alpha) W[l] := subs( x= rootf, op(Basis[l]) ) mod p; od; for l from 1 to e do A[l, n + 1 + 1] := ( add(A[l, i]*W[i], i=1..n) + A[l,n+1]*rootg*lcoeff(g,x) ) mod p; od; A[e + 1, n + 1 +1] := p; A; # This version below is LLL with the vectors above the bound b, being # thrown out. lprint("calling LLL at time:", time()); lprint("calling LLL with:", A); L:= LLLInteger(A,'cut_away',b); lprint("done LLL at time:", time()); lprint("number of vectors remaining:", nops(L)); # Here we split the result from LLL into two parts < C B> were C are the coefficients and B is the result # of the traces e:= nops(L); if e = 0 then lprint(" No isomorphism here at:", time()); return 0, []; fi; B:= 'B'; C:= 'C': for i from 1 to e do B[i] := [L[i][n+2]]; C[i] := [seq(L[i][j], j=1..n+1)]; od; B:= Matrix([seq(B[i],i=1..e)]); C:= Matrix([seq(C[i],i=1..e)]); # If the traces are not zero then we continue to multiply them by 10^n and run LLL again while not(verify(B, Matrix(e,1), 'Matrix')) do A:= Matrix(e, n+1+1); A[1..e,1..n+1] := C; A[1..e, n+2..n+2] := (10^n)*B; lprint("calling LLL after multiplying by 10^n at time:", time()); L:= LLLInteger(A,'cut_away',b); lprint("done with LLL 10^n at time:", time()); e:= nops(L); if e = 0 then lprint(" No isomorphism here at:", time()); return 0, []; fi; B:='B'; C:='C'; for i from 1 to e do B[i] := [L[i][n+2]]; C[i] := [seq(L[i][j], j=1..n+1)]; od; B:= Matrix([seq(B[i],i=1..e)]); C:= Matrix([seq(C[i],i=1..e)]); od; if e = 1 then lprint("Found potential isomorphism at:", time()); return 1, L; fi; lprint("Needs more Hensel Lifting at:", time()); return e, C; end: NumberOfAutomorphisms := proc(f,x) local v,i; v := evala(Factors(f, RootOf(f,x) )); add(`if`(degree(i[1],x)=1, 1, 0), i=v[2]) end: testing := proc(n,p) # where n is the number of primes and p is where the primes #start at local i, j, f, g, M, N; global Errors; options trace; for i from 4 to 8 do Errors[i] := {}; for j from 1 to nops(L[i]) do f := L[i][j]: #g := NewMinpoly(f,x): M[j] :=traperror(iso1p(f,f,x,n,p)); if M[j] = lasterror then if M[j] = "They appear to be Galois" then Errors[i] := Errors[i] union {f}; else error M[j] fi; else [f, M[j], NumberOfAutomorphisms(f,x)] fi; od; od; end: f1, f2 := -2571-4700*x+17712*x^5-31634*x^4-16827*x^3-20032*x^2+10508*x^6+22675*x^7-59352*x^8-89565*x^9- 97744*x^10-41533*x^11-32024*x^12+3140*x^13-1622*x^14+5011*x^15-1843*x^16+2760*x^17-1808*x^18+ 1002*x^19-624*x^20+273*x^21-102*x^22+40*x^23-10*x^24+x^25, 1997-4764*x-26224*x^5+25016*x^4- 12989*x^3+12358*x^2+70324*x^6-40059*x^7+91800*x^8-19093*x^9+87762*x^10-44799*x^11+31952*x^12- 6608*x^13+2914*x^14+6821*x^15+3457*x^16+3416*x^17+2228*x^18+1082*x^19+634*x^20+273*x^21+102*x^22 +40*x^23+10*x^24+x^25: # f1 = -2571-4700*y+17712*y^5-31634*y^4-16827*y^3-20032*y^2+10508*y^6+22675*y^7-59352*y^8-89565*y^9- 97744*y^10-41533*y^11- 32024*y^12+3140*y^13-1622*y^14+5011*y^15-1843*y^16+2760*y^17-1808*y^18+1002*y^19-624*y^20+273*y^21-102*y^22+40*y^23-10*y^24+y^25; # f2 = 1997-4764*y-26224*y^5+25016*y^4-12989*y^3+12358*y^2+70324*y^6-40059*y^7+91800*y^8-19093*y^9+87762*y^10-44799*y^11+31952*y^12-6608*y^13+2914*y^14+6821*y^15+3457*y^16+3416*y^17+2228*y^18+1082*y^19+634*y^20+273*y^21+102*y^22+40*y^23+10*y^24+y^25; g,h := -42885793067858008650*x-161285249796825311495*x^5-429207332210687640181*x^4+ 264147194777000152867*x^3+6032198632961699729*x^2+1332292171242725153638*x^6-\ 2802452375464033976482*x^7+4191186502263628451150*x^8-4343360968383689825174*x ^9+2762117617850418790424*x^10+528408698276686662696*x^11-\ 4525061881234027826296*x^12+6920777688226436002712*x^13-6463590834754261173232 *x^14+4352260622811059705104*x^15-2131245874043847615600*x^16+ 815826517614521346000*x^17-384197583430502150000*x^18+292282923494920350000*x^ 19-249634002590534050000*x^20+180537842164766250000*x^21-102418940816662500000 *x^22+44254465332187500000*x^23-12927273797812500000*x^24+2174026154062500000* x^25+13774402803823804220, -42885793067858008650*x-161285249796825311495*x^5+429207332210687640181*x^4+ 264147194777000152867*x^3-6032198632961699729*x^2-1332292171242725153638*x^6-\ 2802452375464033976482*x^7-4191186502263628451150*x^8-4343360968383689825174*x ^9-2762117617850418790424*x^10+528408698276686662696*x^11+ 4525061881234027826296*x^12+6920777688226436002712*x^13+6463590834754261173232 *x^14+4352260622811059705104*x^15+2131245874043847615600*x^16+ 815826517614521346000*x^17+384197583430502150000*x^18+292282923494920350000*x^ 19+249634002590534050000*x^20+180537842164766250000*x^21+102418940816662500000 *x^22+44254465332187500000*x^23+12927273797812500000*x^24+2174026154062500000* x^25-13774402803823804220: f := x^14+2*x^13-5*x^12-184*x^11-314*x^10+474*x^9+1760*x^8+1504*x^7-400*x^6-1478*x^5-818*x^4+73*x^3+260*x^2+121*x+23 ; # In order to test the isomorphism # L:=factor(f1,RootOf(f2)): # map(degree, [op(L)],x); #pick the one with degree 1 in x # [op(L)][-1]: # iso:= coeff(%,x,0)/coeff(%,x,1): # k:=evala(eval(diff(f2,x), x=RootOf(f2))): # subs(RootOf(f2)=y,evala(Expand(iso*k)));