home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / ALGEBRA.PAK / COUNTROO.M next >
Encoding:
Text File  |  1992-07-29  |  3.7 KB  |  137 lines
  1.  
  2. (* Copyright 1988 Wolfram Research, Inc. *)
  3.  
  4. (*:Version: Mathematica 2.0 *)
  5.  
  6. (*:Name: Algebra`CountRoots` *)
  7.  
  8. (*:Title: Sturm's Method for Counting Real Roots in an Interval *) 
  9.  
  10. (*:Author:
  11.     Dan Grayson, modified September, 1990 by Bruce Sawhill,
  12.         modified October, 1990 by Hon Wah Tam
  13.         modified December, 1990 by E C Martin
  14. *)
  15.  
  16. (*:Keywords:
  17.     Sturm's method, polynomial, real roots, endpoints
  18. *)
  19.  
  20. (*:Requirements: none. *)
  21.  
  22. (*:Warnings: Not working for complex polynomials. *)
  23.  
  24. (*:Sources:
  25.     The Theory of Equations,
  26.     by W.S. Burnside and A.W. Panton, S. Chand & Co., 1972.
  27. *)
  28.  
  29. (*:Summary:
  30. For a given real polynomial f[x] and a given interval a <= x <= b,
  31. CountRoots computes the number of roots of f[x] within the given
  32. interval, using Sturm's method.
  33. *)
  34.  
  35. BeginPackage["Algebra`CountRoots`"]
  36.  
  37. CountRoots::usage = "CountRoots[f,{x,a,b}] uses Sturm's method
  38.     to compute the number of zeros of the real polynomial f 
  39.     on the interval x->[a,b].  The endpoints may be 
  40.     infinite.  Duplicated roots are only counted once."
  41.  
  42. Begin["Algebra`CountRoots`Private`"]
  43.  
  44. CountRoots::npoly = "`1` is not a polynomial in `2`."
  45. CountRoots::llim = "Lower limit `1` is neither a real number nor -Infinity."
  46. CountRoots::ulim = "Upper limit `1` is neither a real number nor Infinity."
  47. CountRoots::intv = "Lower limit `1` is not less than upper limit `2`."
  48.  
  49. CountRoots[ x___] := 
  50.     With[{result = CountRoots0[x]}, result /; result =!= $Failed] 
  51.  
  52.  
  53. OmitZeroes[m_] := Select[m,#!=0&]
  54.  
  55. CountChanges[m_] := 
  56.     Module[{i,len=Length[m],num=0} ,
  57.         Do[ If[m[[i]] != m[[i+1]], num++], {i,len-1} ]; num
  58.     ]
  59.  
  60. CountSignChanges[m_] := CountChanges[OmitZeroes[Sign[m]]]
  61.  
  62. TopCoeff[f_,x_] := Coefficient[f,x,Exponent[f,x]]
  63.  
  64. SignedTopCoeff[f_,x_] := 
  65.     Module[{d=Exponent[f,x]},
  66.          Coefficient[f,x,d] (-1)^d 
  67.      ]
  68.  
  69. SturmSequence[f_,g_,x_] := Module[ { seq = {f}, p=f, q=g} ,
  70.     While[ If[q!=0,True,False,True] ,
  71.         seq = Append[seq,q];
  72.         {p,q} = {q,Expand[-PolynomialRemainder[p,q,x]]};
  73.         ];
  74.         seq
  75.     ]
  76.  
  77. SturmSequence[f_,x_] := SturmSequence[f,D[f,x],x]
  78.  
  79. CountRoots0[f_,{x_,a_,b_}] :=
  80.     Which[ !PolynomialQ[f,x],
  81.         Message[CountRoots::npoly,f,x]; $Failed,
  82.            ( (!NumberQ[a] || SameQ[Head[a],Complex]) &&
  83.             !SameQ[a, -Infinity] ),
  84.          Message[CountRoots::llim,a]; $Failed,  
  85.            ( (!NumberQ[b] || SameQ[Head[b],Complex]) &&
  86.             !SameQ[b, Infinity] ),
  87.         Message[CountRoots::ulim,b]; $Failed,
  88.            True, countroots[f,{x,a,b}] ]
  89.  
  90. CountRoots0[ _] := $Failed
  91.  
  92. countroots[f_,{x_,a_,b_}] := 
  93.     Module[{seq = SturmSequence[f,x] } ,
  94.             If[ !TrueQ[ a <= b ],
  95.                 Message[CountRoots::intv,a,b]; $Failed,
  96.                 CountSignChanges[ seq/.x->a ] 
  97.                 -CountSignChanges[ seq/.x->b ]+
  98.                 If[(f/.x->a)==0,1,0] ]
  99.     ] /; NumberQ[a] && NumberQ[b]
  100.     
  101. countroots[f_,{x_,-Infinity,b_}] := 
  102.     Module[{seq = SturmSequence[f,x]} ,
  103.         If[ !TrueQ[ -Infinity < b ],
  104.             Message[CountRoots::intv,-Infinity,b]; $Failed,
  105.             CountSignChanges[ SignedTopCoeff[#,x]& /@ seq ] 
  106.             -CountSignChanges[ seq/.x->b ] ]
  107.     ] /; NumberQ[b]
  108.     
  109. countroots[f_,{x_,a_,Infinity}] := 
  110.     Module[{seq = SturmSequence[f,x]} ,
  111.         If[ !TrueQ[ a < Infinity ],
  112.             Message[CountRoots::intv,a,Infinity]; $Failed,
  113.             CountSignChanges[ seq/.x->a ] 
  114.             -CountSignChanges[ TopCoeff[#,x]& /@ seq ]+
  115.             If[(f/.x->a)==0,1,0] ]
  116.     ] /; NumberQ[a]
  117.     
  118. countroots[f_,{x_,-Infinity,Infinity}] := 
  119.     Module[{seq = SturmSequence[f,x]} ,
  120.         CountSignChanges[ SignedTopCoeff[#,x]& /@ seq ] 
  121.         -CountSignChanges[ TopCoeff[#,x]& /@ seq ]
  122.     ]
  123.     
  124. End[]  (* Algebra`CountRoots`Private` *)
  125.  
  126. Protect[CountRoots]
  127.  
  128. EndPackage[]  (* Algebra`CountRoots` *)
  129.  
  130. (*:Limitations: Real intervals only. *)
  131.  
  132. (*:Examples:
  133. CountRoots[ x^3 + 4 x^2 - x + 7,{x,-10,10}]
  134.  
  135. *)
  136.  
  137.