> # Excerpts from a MUG posting > # > # Brendan McKay writes: > # > Suppose you have an algebraic number expressed as an > # > formula involving surds. Possibly deeply nested > # > surds. Is there an easy way to find the minimal > # > integer polynomial it satisfies? > > # Date: Tue, 26 Mar 1996 11:45:29 -0800 > # From: monagan@cecm.sfu.ca (Michael Monagan) > # To: maple-list@daisy.uwaterloo.ca > # Subject: minimal polynomial for algebraic number > # > # Yes, there is -- well it's not necessarilly computationally easy, > # but it is conceptionally easy. It requires factorization of resultants > # over Q, then some way of testing whether an expression involving > # surds is NON-zero in order to identify the minimal polynomial > # from the other factors. > # I've written a program to do this (below) -- Brendan, it > # may get bogged down in indentifying the minimal polynomial; > # I really need complex interval arithmetic to do this efficiently. > # Also, in Release 4 `radnormal/radbasis` will do this for you. > # Just look at the code for factor to see how it is used. > # Mike > > minpoly( sqrt(2)+sqrt(3)+sqrt(6) ); > # Error, (in minpoly) 2nd argument required: must be the variable > minpoly( sqrt(2)+sqrt(3)+sqrt(6), z ); > # > # 4 2 4 2 > # minpoly = [1, [[z - 22 z + 48 z - 23, 1], [z - 22 z - 48 z - 23, 1]]] > # > # -6 > # [537.1923951, .1*10 ] > # > # [537.1923951] > # > # 4 2 > # z - 22 z - 48 z - 23 > # > (4-2*3^(1/2))^(1/2); > # > # 1/2 1/2 > # (4 - 2 3 ) > # > minpoly(",z); > # > # 2 2 > # minpoly = [1, [[z + 2 z - 2, 1], [z - 2 z - 2, 1]]] > # > # -8 > # [.2*10 , 2.928203230] > # > # [2.928203230] > # > # 2 > # z + 2 z - 2 > # > # > # > minpoly := proc(a,x:name,K) local A,m,z,rof,p,r; > > if nargs=1 then ERROR(`2nd argument required: must be the variable`) fi; > if nargs=2 then RETURN( minpoly(a,x,{}) ) fi; > if type(K,radical) or type(K,RootOf) then RETURN( minpoly(a,x,{K}) ) fi; > > if hastype(a,radical) then > A := convert([a,x,K],RootOf); > m := minpoly(op(A)); > m := convert(m,radical); > m := readlib(factors)(m); > print(minpoly=m); > m := selectFactor(a,x,map(x -> x[1],m[2])); > RETURN(collect(m,x)) > fi; > > rof := [op(indets(a,RootOf) minus K)]; > if nops(rof)=0 then > if not has(a,x) then RETURN(x-a) fi; > # make it square free > r := evala(Gcd(a,diff(a,x))); > if not has(r,x) then RETURN(a) fi; > evala(Divide(a,r,'m')); > RETURN(m) > fi; > > rof := sort( rof, proc(x,y) evalb(length(x)>length(y)) end ); > rof := op(1,rof); > m := op(1,rof); > m := frontend( subs, [_Z=alpha,m], [{`+`,`*`,`=`}, {}] ); > p := subs(rof=alpha,a); > if not has(p,x) then p := x-p fi; > r := resultant(p,m,alpha); > r := r/lcoeff(r,x); > r := evala(Expand(r)); > RETURN(minpoly(r,x,K)); > > if has(a,x) then RETURN(r) fi; > #f := readlib(factors)(r); > #for f in f[2] do > # print(a,f[1]); > # if evala( subs(x=a,f[1]) ) = 0 then RETURN( minpoly(f[1],x) ) fi; > #od; > #FAIL > > end: > > selectFactor := proc(a,x,mp) local k,e,m,i,s; > if true then > Digits := 10; > do > e := map(abs@evalf,subs(x=a,mp)); > print(e); > m := min( op(e) ); > member(m,e,'k'); > if m=0 then RETURN(mp[k]) fi; > e := subsop(k=NULL,e); > print(e); > e := map( proc(x,m) if x/m > 1000 then NULL else x fi end, e, m ); > if nops(e) = 0 then RETURN(mp[k]) fi; > Digits := 2*Digits > od; > fi; > > if true then > for i to nops(mm) do > if readlib(radnormal)(mm[i])=0 then RETURN(m[i]) fi; > od; > RETURN(FAIL) > fi; > Digits := 5; > s := map((f,x,a) -> evalr(Signum(subs(x=a,f))), m,x,a); > while nops(select(has,s,FAIL))>1 do > Digits := 2*Digits; > s := map((f,x,a) -> evalr(Signum(subs(x=a,f))), m,x,a); > od; > for i do > if has(s[i],FAIL) then RETURN(m[i]) fi; > od; > ERROR() > end: