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