BeginPackage["Sturm`", "Polynom`", "Euclid`"]
(* Copyright: Erich Kaltofen, 4/27/1990.
   Permission to use provided the copyright notice is not removed. *)

Sturm`STrace::usage =
"Type \"Sturm`STrace[s__]:=Print[s]\" to trace this package."

SturmSequence::usage =
"SturmSequence[f_, x_] returns the quotient representation of the
 Sturm sequence for f, i.e.,\n
 \t\t {-1, Q[[t+1]], U[[t]], Q[[t]],..., U[[3]], Q[[3]]}\n.
 such that\n
 \t\t F[[i]] = Q[[i+2]] F[[i+1]] + U[[i+2]] F[[i+2]],\n
 where F[[1]]=f, F[[2]]=D[f, x], and U[[i+2]] are negative units,
 for all 1 <= i <= t-1; F[[t+1]]==0."

SignSequence::usage =
"SignSequence[ss_, x_, a_] returns the sign sequence
 corresponding to the Sturm sequence ss (polynomials in x) at a."

SignVariations::usage =
"SignVariations[signs_] returns the number of sign variations in the
 sign sequence signs."

CountRealRoots::usage=
"CountRealRoots[f_, x_, a_, b_] returns the number of real roots
 of f in x which are in the closed intervall a<=x<=b."

RealRoots::usage=
"RealRoots[f_, x_, precision_] returns a list {r,...} of rational
 numbers such that each real root of f (rational polynomial in x) is
 within precision of a number in the list, meaning
 0 <= r-root < precision. Several r can be equal, but correpond to
 different distinct roots."

Clear["Sturm`IntervalRoots"]
IntervalRoots::usage=
"IntervalRoots[ss_, x_, a_, b_, va_, vb_, precision_]"

Begin["`private`"]

SturmSequence[f_, x_]:=
(* Compute a Sturm Sequence for f in x. *)
Block[{localqr, localnu, fp, cf},
   (* used quotient/remainder function *)
   localqr[a_, b_]:= (* Map[Expand, PolynomialQuotientRemainder[a, b, x] ];
                        Impossible under Version 2.0! *)
                     List[Expand[PolynomialQuotient[a, b, x] ],
                          Expand[PolynomialRemainder[a, b, x] ] ];
   (* used normalization/unit function *)
   localnu[a_] := Block[{nu},
                        nu = Polynom`MonicLC[a, x];
                        (* always normalize such that unit is negative *)
                        If[nu[[2]] > 0,
                           nu = List[Expand[-nu[[1]] ], -nu[[2]] ] 
                          ];
                        Return[nu]  ];                        
   fp := Expand[D[f, x]];
   STrace[StringForm["f = ``, D[f, ``]=``", f, x, fp] ];
   cf = Euclid`ContFractGeneric[f, fp, localqr, localnu];
   (* STrace[StringForm["cf = ``", cf] ]; *)
   Return[Drop[Reverse[cf], 1] ] 
] (* end SturmSequence *)

SignSequence[ss_, x_, a_]:=
(* Compute the sign sequence for the Sturm sequence ss
   in the variable x at the point a. *)
Block[{X, Y, seq, U, Q, temp, signs},
   X = 1; Y = 0; (* multiplier of 0 *)
   seq = ss; signs = {1};
   While[seq != {},
            U = seq[[1]]; Q = (seq[[2]] /. x->a);
            temp = X; X = Q X + U Y; Y = temp;
            seq = Drop[seq, 2]; PrependTo[signs, Sign[X] ];
        ];
   Return[signs]
] (* end SignSequence *)

SignVariations[signs_]:=
Block[{ct = 0, seq, s},
   (* Remove 0s *)
   seq = Select[signs, Function[#!=0] ];
   While[seq != {},
            s = seq[[1]]; seq = Rest[seq];
            If[seq != {}, If[s != seq[[1]], ct++] ]
        ];
   Return[ct]
] (* end SignVariations *)

CountRealRoots[f_, x_, a_, b_]:=
Block[{ss, va, vb},
  ss = SturmSequence[f, x];
  STrace[StringForm["ss = ``", ss] ];
  va = SignVariations[ SignSequence[ ss, x, a ] ];
  STrace[StringForm["va = ``", va] ];
  vb = SignVariations[ SignSequence[ ss, x, b ] ];
  STrace[StringForm["vb = ``", vb] ];
  (* Include possible root at a *)
  If[(f/.x->a)==0, STrace[StringForm["Root at a"] ]; va=va+1];
  Return[va - vb]
]

RealRoots[f_, x_, precision_]:=
Block[{ss, a, b, va, vb},
  ss = SturmSequence[f, x];
  STrace[StringForm["ss = ``", ss] ];
  b = Polynom`RootBound[f, x]; a = -b;
  va = SignVariations[ SignSequence[ ss, x, a ] ];
  vb = SignVariations[ SignSequence[ ss, x, b ] ];
  STrace[ StringForm["a = ``, va =``, b = ``, vb =``", a, va, b, vb] ];
  If[va - vb < 1, Return[{}] ];
  Return[IntervalRoots[ss, x, a, b, va, vb, precision] ]
] (* end RealRoots *)

IntervalRoots[ss_, x_, a_, b_, va_, vb_, precision_]:=
(* Returns a list {r, ...} of roots, a < r <= b, such that all
   real roots given of f (given by the Sturm sequence ss in x)
   lie within precision;
   va and vb are the sign variations of ss at a and b, resp.,
   with va - vb > 0. *)
Block[{mid, vmid, lint = {}, rint ={}},
   STrace[ StringForm["Interval searched: a =``, b =``, va - vb =``",
                     a, b, va - vb] ];
   If[b - a < precision,
      (* b is an approximation for va - vb roots *)
      Do[PrependTo[lint, b], {va - vb}];
      Return[ lint ] ];
   mid = (a+b)/2;
   vmid = SignVariations[ SignSequence[ ss, x, mid ] ];
   If[va - vmid > 0,
      lint = IntervalRoots[ss, x, a, mid, va, vmid, precision] ];
   If[vmid - vb > 0,
      rint = IntervalRoots[ss, x, mid, b, vmid, vb, precision] ];
   Return[Join[lint, rint] ]
] (* end IntervalRoots *)

End[] (* `private` *)
EndPackage[] (* Sturm` *)
Print["Package Sturm` loaded"]