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