## We load first the Maple routine : minpoly (recent version of Monagan and
## Fee). Precision used : 64 decimal digits.
## That program uses advanced techniques of linear algebra and this is
## why we read the library package (lattice ).
## 
##
##


readlib(lattice):
norma:=proc(a,n) coeff(minpoly(a,n,x),x,0);end:
rtp:=proc(A)  global 
_x;Minpoly(evala(convert(A,RootOf)),_x);end:#radical_to_poly

minp:=proc() global d; local p;
	if nargs>2 then d:=args[3] else d:=30 fi;
Digits:=d;print(`Precision`,d);
p:=factor(minpoly(args[1],args[2]));
	print(`Error`,subs(_X=args[1],p));
p;
end:

Minpoly := proc(a,x:name,K) local A,m,rof,p,alpha,r,i,s;

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:
       
minpoly := proc(r, n:posint, x, w)
    local i, j, r1, A, L;
    option `Copyright 1991 by the University of Waterloo`;
    if nargs=2 then RETURN( minpoly(r,n,_X,10^(Digits-2)) ) fi;
    if nargs=3 and type(x,name) then RETURN( 
minpoly(r,n,x,10^(Digits-2))) fi;
    if nargs=3 and type(x,numeric) then RETURN( minpoly(r,n,_X,x) ) fi;
    r1 := evalf(r);
    if not type(r1,complex(numeric)) then ERROR(`wrong type of arguments`) fi;
    r1 := convert(r1,rational,exact);
    if type(r1,numeric) then
       A := [seq([seq(delta(i, j),j=0..n), w*Re(r1^i)], i=0..n)];
    else # complex numeric case
       A := [seq([seq(delta(i, j),j=0..n),w*Re(r1^i),w*Im(r1^i)], i=0..n)];
    fi;
    # Important: round the inputs to Digits precision
    A := convert(evalf(A),rational,exact);
    userinfo(2,minpoly,`Input matrix to LLL is`,print(A));
    L := lattice(A);
    userinfo(2,minpoly,`Lattice computed is`,print(L));
    sum('L[1][i+1]*x^i','i'=0..n)
end:

lindep := proc(v:list,x:name) local v1,w,A,n,L,i,j;
    v1 := convert( evalf(v), rational, exact );
    w := 10^(Digits-2);
    n := nops(v)-1;
    if type(v1,list(numeric)) then
       A := [seq([seq(delta(i, j),j=0..n), w*v1[i+1]],i=0..n)];
    elif type(v1,list(complex(numeric))) then
       A := [seq([seq(delta(i, j),j=0..n),w*Re(v1[i+1]),w*Im(v1[i+1])],
                 i=0..n)];
    else ERROR(`bad input`)
    fi;
    # Important: round the inputs to Digits precision
    A := convert(evalf(A),rational,exact);
    userinfo(2,lindep,`Input matrix is`,print(A));
    L := lattice(A,integer);
    userinfo(2,lindep,`Lattice computed is`,print(L));
    v1 := L[1][1..n+1];
    if nargs=2 then add(v1[i+1]*x[i], i=0..n);
    else v1
    fi
end:

Digits:=64:
f:=proc(n) local k;evalf(sum(1/16**k/(8*k+n),k=0..infinity)):end:
F:=proc(n) local k;Sum(1/16**k/(8*k+n),k=0..infinity):end:
liste:=[f(1),f(2),f(3),f(4),f(5),f(6),f(7),f(8),Pi]:
lindep(");vector_found:=["];
` Pi`=sum(F(i)*vector_found[i],i=1..8);