001 /*
002 * $Id: RealRootsSturm.java 3677 2011-07-04 20:29:49Z kredel $
003 */
004
005 package edu.jas.root;
006
007
008 import java.util.ArrayList;
009 import java.util.List;
010
011 import org.apache.log4j.Logger;
012
013 import edu.jas.arith.Rational;
014 import edu.jas.poly.GenPolynomial;
015 import edu.jas.poly.PolyUtil;
016 import edu.jas.structure.RingElem;
017 import edu.jas.structure.RingFactory;
018
019
020 /**
021 * Real root isolation using Sturm sequences.
022 * @param <C> coefficient type.
023 * @author Heinz Kredel
024 */
025 public class RealRootsSturm<C extends RingElem<C> & Rational> extends RealRootAbstract<C> {
026
027
028 private static final Logger logger = Logger.getLogger(RealRootsSturm.class);
029
030
031 private static boolean debug = logger.isDebugEnabled();
032
033
034 /**
035 * Sturm sequence.
036 * @param f univariate polynomial.
037 * @return a Sturm sequence for f.
038 */
039 public List<GenPolynomial<C>> sturmSequence(GenPolynomial<C> f) {
040 List<GenPolynomial<C>> S = new ArrayList<GenPolynomial<C>>();
041 if (f == null || f.isZERO()) {
042 return S;
043 }
044 if (f.isConstant()) {
045 S.add(f.monic());
046 return S;
047 }
048 GenPolynomial<C> F = f;
049 S.add(F);
050 GenPolynomial<C> G = PolyUtil.<C> baseDeriviative(f);
051 while (!G.isZERO()) {
052 GenPolynomial<C> r = F.remainder(G);
053 F = G;
054 G = r.negate();
055 S.add(F/*.monic()*/);
056 }
057 //System.out.println("F = " + F);
058 if (F.isConstant()) {
059 return S;
060 }
061 // make squarefree
062 List<GenPolynomial<C>> Sp = new ArrayList<GenPolynomial<C>>(S.size());
063 for (GenPolynomial<C> p : S) {
064 p = p.divide(F);
065 Sp.add(p);
066 }
067 return Sp;
068 }
069
070
071 /**
072 * Isolating intervals for the real roots.
073 * @param f univariate polynomial.
074 * @return a list of isolating intervals for the real roots of f.
075 */
076 @Override
077 public List<Interval<C>> realRoots(GenPolynomial<C> f) {
078 List<Interval<C>> R = new ArrayList<Interval<C>>();
079 if (f == null || f.isConstant()) {
080 return R;
081 }
082 if (f.isZERO()) {
083 C z = f.ring.coFac.getZERO();
084 R.add(new Interval<C>(z));
085 return R;
086 }
087 if ( f.degree(0) == 1L ) {
088 C z = f.monic().trailingBaseCoefficient().negate();
089 R.add(new Interval<C>(z));
090 return R;
091 }
092 GenPolynomial<C> F = f;
093 C M = realRootBound(F); // M != 0, since >= 2
094 Interval<C> iv = new Interval<C>(M.negate(), M);
095 //System.out.println("iv = " + iv);
096 List<GenPolynomial<C>> S = sturmSequence(F);
097 //System.out.println("S = " + S);
098 //System.out.println("f_S = " + S.get(0));
099 List<Interval<C>> Rp = realRoots(iv, S);
100 R.addAll(Rp);
101 return R;
102 }
103
104
105 /**
106 * Isolating intervals for the real roots.
107 * @param iv interval with f(left) * f(right) != 0.
108 * @param S sturm sequence for f and I.
109 * @return a list of isolating intervals for the real roots of f in I.
110 */
111 public List<Interval<C>> realRoots(Interval<C> iv, List<GenPolynomial<C>> S) {
112 List<Interval<C>> R = new ArrayList<Interval<C>>();
113 GenPolynomial<C> f = S.get(0); // squarefree part
114 if ( f.isZERO() ) {
115 C z = f.leadingBaseCoefficient();
116 if ( !iv.contains(z) ) {
117 throw new IllegalArgumentException("root not in interval: f = " + f + ", iv = " + iv + ", z = " + z);
118 }
119 Interval<C> iv1 = new Interval<C>(z);
120 R.add(iv1);
121 return R;
122 }
123 if ( f.isConstant() ) {
124 throw new IllegalArgumentException("f has no root: f = " + f + ", iv = " + iv);
125 }
126 if ( f.degree(0) == 1L ) {
127 C z = f.monic().trailingBaseCoefficient().negate();
128 if ( !iv.contains(z) ) {
129 throw new IllegalArgumentException("root not in interval: f = " + f + ", iv = " + iv + ", z = " + z);
130 }
131 Interval<C> iv1 = new Interval<C>(z);
132 R.add(iv1);
133 return R;
134 }
135 //System.out.println("iv = " + iv);
136 // check sign variations at interval bounds
137 long v = realRootCount(iv, S);
138 //System.out.println("v = " + v);
139 if (v == 0) {
140 return R;
141 }
142 if (v == 1) {
143 R.add(iv);
144 return R;
145 }
146 // now v > 1
147 // bi-sect interval, such that f(c) != 0
148 C c = bisectionPoint(iv, f);
149 //System.out.println("c = " + c);
150 // recursion on both sub-intervals
151 Interval<C> iv1 = new Interval<C>(iv.left, c);
152 Interval<C> iv2 = new Interval<C>(c, iv.right);
153 List<Interval<C>> R1 = realRoots(iv1, S);
154 //System.out.println("R1 = " + R1);
155 if (debug) {
156 logger.info("R1 = " + R1);
157 }
158 List<Interval<C>> R2 = realRoots(iv2, S);
159 //System.out.println("R2 = " + R2);
160 if (debug) {
161 logger.info("R2 = " + R2);
162 }
163
164 // refine isolating intervals if adjacent
165 if (R1.isEmpty()) {
166 R.addAll(R2);
167 return R;
168 }
169 if (R2.isEmpty()) {
170 R.addAll(R1);
171 return R;
172 }
173 iv1 = R1.get(R1.size() - 1); // last
174 iv2 = R2.get(0); // first
175 if (iv1.right.compareTo(iv2.left) < 0) {
176 R.addAll(R1);
177 R.addAll(R2);
178 return R;
179 }
180 // now iv1.right == iv2.left
181 //System.out.println("iv1 = " + iv1);
182 //System.out.println("iv2 = " + iv2);
183 R1.remove(iv1);
184 R2.remove(iv2);
185 while (iv1.right.equals(iv2.left)) {
186 C d1 = bisectionPoint(iv1, f);
187 C d2 = bisectionPoint(iv2, f);
188 Interval<C> iv11 = new Interval<C>(iv1.left, d1);
189 Interval<C> iv12 = new Interval<C>(d1, iv1.right);
190 Interval<C> iv21 = new Interval<C>(iv2.left, d2);
191 Interval<C> iv22 = new Interval<C>(d2, iv2.right);
192
193 boolean b11 = signChange(iv11, f);
194 boolean b12 = signChange(iv12, f);
195 boolean b21 = signChange(iv21, f);
196 boolean b22 = signChange(iv22, f);
197 if (b11) {
198 iv1 = iv11;
199 if (b22) {
200 iv2 = iv22;
201 } else {
202 iv2 = iv21;
203 }
204 break; // done, refine
205 }
206 if (b22) {
207 iv2 = iv22;
208 if (b12) {
209 iv1 = iv12;
210 } else {
211 iv1 = iv11;
212 }
213 break; // done, refine
214 }
215 iv1 = iv12;
216 iv2 = iv21;
217 //System.out.println("iv1 = " + iv1);
218 //System.out.println("iv2 = " + iv2);
219 }
220 R.addAll(R1);
221 R.add(iv1);
222 R.add(iv2);
223 R.addAll(R2);
224 return R;
225 }
226
227
228 /**
229 * Number of real roots in interval.
230 * @param iv interval with f(left) * f(right) != 0.
231 * @param S sturm sequence for f and I.
232 * @return number of real roots of f in I.
233 */
234 public long realRootCount(Interval<C> iv, List<GenPolynomial<C>> S) {
235 // check sign variations at interval bounds
236 GenPolynomial<C> f = S.get(0); // squarefree part
237 //System.out.println("iv = " + iv);
238 RingFactory<C> cfac = f.ring.coFac;
239 List<C> l = PolyUtil.<C> evaluateMain(cfac, S, iv.left);
240 List<C> r = PolyUtil.<C> evaluateMain(cfac, S, iv.right);
241 long v = RootUtil.<C> signVar(l) - RootUtil.<C> signVar(r);
242 //System.out.println("v = " + v);
243 if (v < 0L) {
244 v = -v;
245 }
246 return v;
247 }
248
249
250 /**
251 * Number of real roots in interval.
252 * @param iv interval with f(left) * f(right) != 0.
253 * @param f univariate polynomial.
254 * @return number of real roots of f in I.
255 */
256 @Override
257 public long realRootCount(Interval<C> iv, GenPolynomial<C> f) {
258 if (f == null || f.isZERO() || f.isConstant()) {
259 return 0L;
260 }
261 List<GenPolynomial<C>> S = sturmSequence(f);
262 return realRootCount(iv, S);
263 }
264
265
266 /**
267 * Invariant interval for algebraic number sign.
268 * @param iv root isolating interval for f, with f(left) * f(right) < 0.
269 * @param f univariate polynomial, non-zero.
270 * @param g univariate polynomial, gcd(f,g) == 1.
271 * @return v with v a new interval contained in iv such that g(w) != 0 for w in v.
272 */
273 @Override
274 public Interval<C> invariantSignInterval(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g) {
275 Interval<C> v = iv;
276 if (g == null || g.isZERO()) {
277 return v;
278 }
279 if (g.isConstant()) {
280 return v;
281 }
282 if (f == null || f.isZERO() || f.isConstant()) { // ?
283 return v;
284 }
285 List<GenPolynomial<C>> Sg = sturmSequence(g.monic());
286 Interval<C> ivp = invariantSignInterval(iv, f, Sg);
287 return ivp;
288 }
289
290
291 /**
292 * Invariant interval for algebraic number sign.
293 * @param iv root isolating interval for f, with f(left) * f(right) < 0.
294 * @param f univariate polynomial, non-zero.
295 * @param Sg Sturm sequence for g, a univariate polynomial with gcd(f,g) == 1.
296 * @return v with v a new interval contained in iv such that g(w) != 0 for w in v.
297 */
298 public Interval<C> invariantSignInterval(Interval<C> iv, GenPolynomial<C> f, List<GenPolynomial<C>> Sg) {
299 Interval<C> v = iv;
300 GenPolynomial<C> g = Sg.get(0);
301 if (g == null || g.isZERO()) {
302 return v;
303 }
304 if (g.isConstant()) {
305 return v;
306 }
307 if (f == null || f.isZERO() || f.isConstant()) { // ?
308 return v;
309 }
310 RingFactory<C> cfac = f.ring.coFac;
311 C two = cfac.fromInteger(2);
312
313 while (true) {
314 long n = realRootCount(v, Sg);
315 logger.debug("n = " + n);
316 if (n == 0) {
317 return v;
318 }
319 C c = v.left.sum(v.right);
320 c = c.divide(two);
321 Interval<C> im = new Interval<C>(c, v.right);
322 if (signChange(im, f)) {
323 v = im;
324 } else {
325 v = new Interval<C>(v.left, c);
326 }
327 }
328 // return v;
329 }
330
331 }