(*:Copyright: Copyright 1988-1999 Wolfram Research, Inc. *)

(*:Mathematica Version: 3.0 *)

(*:Package Version: 1.0 *)

(*:Name: Algebra`RootIsolation` *)

(*:Title: Polynomial root isolation*) 

(*:Author: Adam Strzebonski *)

(*:Keywords:
	polynomial roots, root isolation, polynomial, complex roots
*)

(*:Requirements: None. *)

(*:Sources:
	An Efficient Algorithm for Infallible Polynomial Complex Roots 
	Isolation, G.E.Collins & W.Krandick, Proceedings of ISSAC'92.
	An implementation of Vincent's theorem, A.G.Akritas, Numerische
	Mathematic 36, 1980, pp.53-62.
*)

(*:Summary:
This package defines functions for counting and isolating real and complex
roots of polynomials.
*)

(*:Discussion:
CountRoots[f,{a,b}] computes the number of roots (multiplicities counted) of 
a polynomial f in the interval (a,b) (for complex numbers a,b interval (a,b) 
is the open rectangle (or open line segment or point) of which a is the 
lower-left vertex and b is the upper-right vertex).  The function uses
the winding number computation for rectangles and Descartes' rule of signs
for intervals.    
RealRootIntervals finds "isolating intervals" for all real roots of a given
rational univariate polynomial or a list of such polynomials. 
ComplexRootIntervals finds "isolating intervals" for all roots of a given
rational univariate polynomial or a list of such polynomials. 
If a root is real, the isolating interval is an openreal interval, if a root
is not real, the isolating interval is an open  rectangle, disjoint from the 
real axis.
For a given algebraic number a, and a given integer prec, ContractInterval 
finds an isolating interval for a which specifies a up to the decimal 
precision prec.
*)

BeginPackage["Algebra`RootIsolation`"]

ComplexRootIntervals::usage = "ComplexRootIntervals[pp] computes isolating intervals for all roots of pp. pp is a rational univariate polynomial or a list of such polynomials. If a root is real, the isolating interval is an open real interval, if a root is not real, the isolating interval is an open rectangle, disjoint from the real axis."

RealRootIntervals::usage = "RealRootIntervals[pp] computes isolating intervals for all real roots of pp. pp is a rational univariate polynomial or a list of such polynomials."
	  
ContractInterval::usage = "ContractInterval[a, prec] returns an interval containing the algebraic number a which specifies a up to the decimal precision prec." 	  
	  	  
CountRoots::usage = "CountRoots[f, {x,a,b}] computes the number of roots (multiplicities counted) of a univariate polynomial f[x] in the interval (a,b) (for complex numbers a,b interval (a,b) is the open rectangle (or open line segment or point) of which a is the lower-left vertex and b is the upper-right vertex)."

Begin["Algebra`RootIsolation`Private`"]

Unprotect[CountRoots]
Unprotect[ComplexRootIntervals]
Unprotect[ContractInterval]
Unprotect[RealRootIntervals]

RootReduce[0];
algebraic[a__] := RootsDump`algebraic[a]
algnumQ[a__] := RootsDump`algnumQ[a]
codeContract[a__] := Algebra`Algebraics`Private`codeContract[a]
countsqf[f_,z_,int_] := Algebra`Algebraics`Private`countsqf[f,z,int]
gaussianRationalQ[a__] := RootsDump`gaussianRationalQ[a]
hshrink[a__] := Algebra`Algebraics`Private`hshrink[a]
redsf[a__] := Algebra`Algebraics`Private`redsf[a]

fsfl[f_] := Select[FactorSquareFreeList[f],!NumberQ[#[[1]]]&]

gooddataQ[f_,{x_,int1_,int2_}] := (NumericQ[f] && (f!=0) || 
         (Variables[f]==={x}) && PolynomialQ[f,x] &&
	 And@@(gaussianRationalQ[#]&/@CoefficientList[f,x])) &&
	 (gaussianRationalQ[int1] || int1===-Infinity && FreeQ[int2,Complex]) && 
	 (gaussianRationalQ[int2] || int2===+Infinity && FreeQ[int1,Complex])

gooddataQ[others__] := False

goodpolyQ[f_] := rationalQ[f] && (f!=0) || (Length[Variables[f]]==1) &&
         And@@(rationalQ[#]&/@CoefficientList[f,Variables[f][[1]]])
	 
goodpolyQ[others__] := False	 

rationalQ[a_] := IntegerQ[a] || (Head[a]===Rational)
 
fl[f_] := If[Expand[f]===1,{{1,1}},FactorList[f]]

polydataQ[ff_] := If[Head[ff]===List, And@@(goodpolyQ[#]&/@ff), goodpolyQ[ff]] 
		                                   	 
polydataQ[others__] := False

ComplexRootIntervals[polys_]/;polydataQ[polys] :=
    Module[{l,ll,n,nn={},fs,i,rts},
	    If[Head[polys]===List,
	       l=Length[polys];
	       Do[fs=Select[Transpose[fl[polys[[i]]]][[1]],
	                    !NumberQ[#]&];
	          nn=Union[nn,Table[{fs[[j]],i},{j,Length[fs]}]],
		  {i,l}];
	       If[nn=={}, Return[{{},{}}]];
	       nn=Transpose[nn]; fs=Union[nn[[1]]]; 
	       l=Length[fs]; ll={};
	       Do[n=Position[nn[[1]],fs[[i]],1];
	          If[Length[n]==1,
		    ll=Append[ll,nn[[2]][[n[[1]][[1]]]]],
		    ll=Append[ll,Table[nn[[2]][[n[[j]][[1]]]],{j,Length[n]}]]],
		  {i,l}];
	       rts=System`Private`ComplexRoots[fs]; 
	       nn={}; l=Length[rts[[2]]];
	       Do[nn=Append[nn,ll[[rts[[2]][[i]]]]], {i,l}];
	       {rts[[1]],nn},
	       fs=Select[Transpose[fl[polys]][[1]],!NumberQ[#]&];
	       System`Private`ComplexRoots[fs][[1]]]]

ContractInterval[a_,prec_Integer]/;algnumQ[a] :=
	If[gaussianRationalQ[a],
           {a,a},
	   codeContract[algebraic[a],prec]]

CountRoots[f_,{x_,a_,b_}]/;gooddataQ[f,{x,a,b}] :=
          Module[{factors=fsfl[f],i,n=0},
	         If[NumericQ[f], Return[0]];
	         Do[n=n+factors[[i]][[2]] countsqf[factors[[i]][[1]],x,{a,b}],
		    {i,Length[factors]}];
		 n]

CountRoots[f_,int_]/;
          gooddataQ[SetPrecision[f,Infinity],SetPrecision[int,Infinity]] :=
          CountRoots[SetPrecision[f,Infinity],SetPrecision[int,Infinity]]
          
RealRootIntervals[polys_]/;polydataQ[polys] :=
     Module[{l,ll,n,nn={},fs,i,rts},
	    If[Head[polys]===List, 
	       l=Length[polys]; 
	       Do[fs=Select[Transpose[fl[polys[[i]]]][[1]],
	                    !NumberQ[#]&];
	          nn=Union[nn,Table[{fs[[j]],i},{j,Length[fs]}]],
		  {i,l}]; 
	       If[nn=={}, Return[{{},{}}]];
	       nn=Transpose[nn]; fs=Union[nn[[1]]];  
	       l=Length[fs]; ll={};
	       Do[n=Position[nn[[1]],fs[[i]],1];
	          If[Length[n]==1,
		    ll=Append[ll,nn[[2]][[n[[1]][[1]]]]],
		    ll=Append[ll,Table[nn[[2]][[n[[j]][[1]]]],{j,Length[n]}]]],
		  {i,l}]; 
	       rts=System`Private`RealRoots[fs]; 
	       nn={}; l=Length[rts[[2]]];
	       Do[nn=Append[nn,ll[[rts[[2]][[i]]]]], {i,l}];
	       {rts[[1]],nn},
	       fs=redsf[polys];
	       System`Private`RealRoots[fs][[1]]]]

End[]  (* Algebra`RootIsolation`Private` *)

SetAttributes[#, {Protected, ReadProtected}]&/@
{CountRoots, ComplexRootIntervals, ContractInterval, RealRootIntervals}

EndPackage[]  (* Algebra`RootIsolation` *)

