# $Source: /u/maple/research/lib/evala/src/RCS/Subfields,v $ # $Notify: mvanhoei@daisy.uwaterloo.ca macro( RAD={'radical','nonreal'}, g_factors=`DEtools/kovacicsols/factors`, shorten=`polytools/shorten`, sort_poly=`polytools/sort_poly`, degKprod=`evala/Subfields/degKprod` ): `evala/Subfields`:=proc(f::{polynom,set},deg::posint,F::set,x::name) local i,p,w,ans,R,G,d,j; option `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved.\ Author: Mark van Hoeij`; i:=indets(f,'RootOf'); if deg<2 then error "deg should be >1" elif f={} then return f elif nargs=2 then return procname(f,deg,i) elif indets([args],RAD)<>{} then p := radfield(indets([args],RAD)); return eval(subs(p[2],procname(op(eval(subs(p[1],[args])))))) fi; i:=i union F; if nargs=3 then p:=indets(f,'name'); if nops(p)<>1 then x # = error message. fi; return procname(f,deg,i,p[1]) elif F<>i then return procname(f,deg,i,x) fi; Normalizer:=`evala/Normal`; if type(f,set) then if f={} or {0}<>(map(degree,f,x) mod deg) then return {} fi; w:=sort_poly([op(f)],x); ans:=NULL; for R in {procname(w[1],deg,F,x,x)} do if nops(f)=1 and degree(f[1],x)=deg then return {R} fi; G:=F union {R}; d:=0; for p in w do if nops(g_factors(p,x,G))=1 then d:=1; break fi od; if d=0 then for p in [ans] do if nops(g_factors(p,x,G))>1 then # Remove doubles d:=1; break fi od fi; if d=0 then ans:=ans,R fi od; {ans} elif nargs=4 then procname({f},deg,F,x) else if not type(f,polynom(anything,x)) then error f, "is not a polynomial" fi; d:=degree(f,x); if not type(d/deg,posint) then return NULL fi; p:=shorten(f,x); R:=RootOf(p,x); if d=deg then return R fi; G:=F union {R}; p:={op(map2(op,1,g_factors(p,x,G)))}; if deg=2 then p:=p minus {p[1]} fi; # to avoid doubles p:=degKprod(p,x,d/deg); ans:=NULL; for j in p do w:=Normalizer(add(coeff(i,x,degree(i,x)-1),i=j)); if not has(w,R) then w:=sort_poly(select(has,[coeffs(collect(convert(j,`*`) ,x,Normalizer),x)],R),x)[1] fi; i:=sqrfree(evala(Norm(x-w,G,F)),x)[2][1][1]; if degree(i,x)>deg then next elif degree(i,x)=deg then ans:=ans,RootOf(shorten(i,x),x) elif not isprime(deg) then error "failed" else w:=sort_poly(select(has,[coeffs(collect(convert(j,`*`) ,x,Normalizer),x)],R),x); w:=[op(w),seq(randpoly(w,'coeffs'=rand(-d..d),'dense', degree=1),d=[1,2,5,9,99,999,10^9])]; for d in w do i:=sqrfree(evala(Norm(x-d,G,F)),x)[2][1][1]; if degree(i,x)>deg then break elif degree(i,x)=deg then ans:=ans,RootOf(shorten(i,x),x); break fi od; if degree(i,x)