001/* 002 * $Id: RealRootsSturm.java 4961 2014-10-17 18:59:39Z 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 @SuppressWarnings("cast") 078 @Override 079 public List<Interval<C>> realRoots(GenPolynomial<C> f) { 080 List<Interval<C>> R = new ArrayList<Interval<C>>(); 081 if (f == null || f.isConstant()) { 082 return R; 083 } 084 if (f.isZERO()) { 085 C z = f.ring.coFac.getZERO(); 086 R.add(new Interval<C>(z)); 087 return R; 088 } 089 if (f.degree(0) == 1L) { 090 C z = f.monic().trailingBaseCoefficient().negate(); 091 R.add(new Interval<C>(z)); 092 return R; 093 } 094 GenPolynomial<C> F = f; 095 C M = realRootBound(F); // M != 0, since >= 2 096 Interval<C> iv = new Interval<C>(M.negate(), M); 097 //System.out.println("iv = " + iv); 098 List<GenPolynomial<C>> S = sturmSequence(F); 099 //System.out.println("S = " + S); 100 //System.out.println("f_S = " + S.get(0)); 101 List<Interval<C>> Rp = realRoots(iv, S); 102 if (logger.isInfoEnabled() && !(((Object) f.ring.coFac) instanceof BigRational)) { 103 //logger.info("realRoots bound: " + iv); 104 logger.info("realRoots: " + Rp); 105 } 106 R.addAll(Rp); 107 return R; 108 } 109 110 111 /** 112 * Isolating intervals for the real roots. 113 * @param iv interval with f(left) * f(right) != 0. 114 * @param S sturm sequence for f and I. 115 * @return a list of isolating intervals for the real roots of f in I. 116 */ 117 public List<Interval<C>> realRoots(Interval<C> iv, List<GenPolynomial<C>> S) { 118 List<Interval<C>> R = new ArrayList<Interval<C>>(); 119 GenPolynomial<C> f = S.get(0); // squarefree part 120 if (f.isZERO()) { 121 C z = f.leadingBaseCoefficient(); 122 if (!iv.contains(z)) { 123 throw new IllegalArgumentException("root not in interval: f = " + f + ", iv = " + iv 124 + ", z = " + z); 125 } 126 Interval<C> iv1 = new Interval<C>(z); 127 R.add(iv1); 128 return R; 129 } 130 if (f.isConstant()) { 131 return R; 132 //throw new IllegalArgumentException("f has no root: f = " + f + ", iv = " + iv); 133 } 134 if (f.degree(0) == 1L) { 135 C z = f.monic().trailingBaseCoefficient().negate(); 136 if (!iv.contains(z)) { 137 return R; 138 //throw new IllegalArgumentException("root not in interval: f = " + f + ", iv = " + iv + ", z = " + z); 139 } 140 Interval<C> iv1 = new Interval<C>(z); 141 R.add(iv1); 142 return R; 143 } 144 //System.out.println("iv = " + iv); 145 // check sign variations at interval bounds 146 long v = realRootCount(iv, S); 147 //System.out.println("v = " + v); 148 if (v == 0) { 149 return R; 150 } 151 if (v == 1) { 152 R.add(iv); 153 return R; 154 } 155 // now v > 1 156 // bi-sect interval, such that f(c) != 0 157 C c = bisectionPoint(iv, f); 158 //System.out.println("c = " + c); 159 // recursion on both sub-intervals 160 Interval<C> iv1 = new Interval<C>(iv.left, c); 161 Interval<C> iv2 = new Interval<C>(c, iv.right); 162 List<Interval<C>> R1 = realRoots(iv1, S); 163 //System.out.println("R1 = " + R1); 164 if (debug) { 165 logger.info("R1 = " + R1); 166 } 167 List<Interval<C>> R2 = realRoots(iv2, S); 168 //System.out.println("R2 = " + R2); 169 if (debug) { 170 logger.info("R2 = " + R2); 171 } 172 173 // refine isolating intervals if adjacent 174 if (R1.isEmpty()) { 175 R.addAll(R2); 176 return R; 177 } 178 if (R2.isEmpty()) { 179 R.addAll(R1); 180 return R; 181 } 182 iv1 = R1.get(R1.size() - 1); // last 183 iv2 = R2.get(0); // first 184 if (iv1.right.compareTo(iv2.left) < 0) { 185 R.addAll(R1); 186 R.addAll(R2); 187 return R; 188 } 189 // now iv1.right == iv2.left 190 //System.out.println("iv1 = " + iv1); 191 //System.out.println("iv2 = " + iv2); 192 R1.remove(iv1); 193 R2.remove(iv2); 194 while (iv1.right.equals(iv2.left)) { 195 C d1 = bisectionPoint(iv1, f); 196 C d2 = bisectionPoint(iv2, f); 197 Interval<C> iv11 = new Interval<C>(iv1.left, d1); 198 Interval<C> iv12 = new Interval<C>(d1, iv1.right); 199 Interval<C> iv21 = new Interval<C>(iv2.left, d2); 200 Interval<C> iv22 = new Interval<C>(d2, iv2.right); 201 202 boolean b11 = signChange(iv11, f); 203 boolean b12 = signChange(iv12, f); // TODO check unnecessary 204 //boolean b21 = signChange(iv21, f); // TODO check unused or unnecessary 205 boolean b22 = signChange(iv22, f); 206 if (b11) { 207 iv1 = iv11; 208 if (b22) { 209 iv2 = iv22; 210 } else { 211 iv2 = iv21; 212 } 213 break; // done, refine 214 } 215 if (b22) { 216 iv2 = iv22; 217 if (b12) { 218 iv1 = iv12; 219 } else { 220 iv1 = iv11; 221 } 222 break; // done, refine 223 } 224 iv1 = iv12; 225 iv2 = iv21; 226 //System.out.println("iv1 = " + iv1); 227 //System.out.println("iv2 = " + iv2); 228 } 229 R.addAll(R1); 230 R.add(iv1); 231 R.add(iv2); 232 R.addAll(R2); 233 return R; 234 } 235 236 237 /** 238 * Number of real roots in interval. 239 * @param iv interval with f(left) * f(right) != 0. 240 * @param S sturm sequence for f and I. 241 * @return number of real roots of f in I. 242 */ 243 public long realRootCount(Interval<C> iv, List<GenPolynomial<C>> S) { 244 // check sign variations at interval bounds 245 GenPolynomial<C> f = S.get(0); // squarefree part 246 //System.out.println("iv = " + iv); 247 RingFactory<C> cfac = f.ring.coFac; 248 List<C> l = PolyUtil.<C> evaluateMain(cfac, S, iv.left); 249 List<C> r = PolyUtil.<C> evaluateMain(cfac, S, iv.right); 250 long v = RootUtil.<C> signVar(l) - RootUtil.<C> signVar(r); 251 //System.out.println("v = " + v); 252 if (v < 0L) { 253 v = -v; 254 } 255 return v; 256 } 257 258 259 /** 260 * Number of real roots in interval. 261 * @param iv interval with f(left) * f(right) != 0. 262 * @param f univariate polynomial. 263 * @return number of real roots of f in I. 264 */ 265 @Override 266 public long realRootCount(Interval<C> iv, GenPolynomial<C> f) { 267 if (f == null || f.isConstant()) { // ? 268 return 0L; 269 } 270 if (f.isZERO()) { 271 C z = f.leadingBaseCoefficient(); 272 if (!iv.contains(z)) { 273 return 0L; 274 } 275 return 1L; 276 } 277 List<GenPolynomial<C>> S = sturmSequence(f); 278 return realRootCount(iv, S); 279 } 280 281 282 /** 283 * Invariant interval for algebraic number sign. 284 * @param iv root isolating interval for f, with f(left) * f(right) < 0. 285 * @param f univariate polynomial, non-zero. 286 * @param g univariate polynomial, gcd(f,g) == 1. 287 * @return v with v a new interval contained in iv such that g(w) != 0 for w 288 * in v. 289 */ 290 @Override 291 public Interval<C> invariantSignInterval(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g) { 292 Interval<C> v = iv; 293 if (g == null || g.isZERO()) { 294 //throw new IllegalArgumentException("g == 0"); 295 return v; 296 } 297 if (g.isConstant()) { 298 return v; 299 } 300 if (f == null || f.isZERO()) { // ? || f.isConstant() 301 throw new IllegalArgumentException("f == 0"); 302 //return v; 303 } 304 List<GenPolynomial<C>> Sg = sturmSequence(g.monic()); 305 Interval<C> ivp = invariantSignInterval(iv, f, Sg); 306 return ivp; 307 } 308 309 310 /** 311 * Invariant interval for algebraic number sign. 312 * @param iv root isolating interval for f, with f(left) * f(right) < 0. 313 * @param f univariate polynomial, non-zero. 314 * @param Sg Sturm sequence for g, a univariate polynomial with gcd(f,g) == 315 * 1. 316 * @return v with v a new interval contained in iv such that g(w) != 0 for w 317 * in v. 318 */ 319 public Interval<C> invariantSignInterval(Interval<C> iv, GenPolynomial<C> f, List<GenPolynomial<C>> Sg) { 320 Interval<C> v = iv; 321 GenPolynomial<C> g = Sg.get(0); 322 if (g == null || g.isZERO()) { 323 return v; 324 } 325 if (g.isConstant()) { 326 return v; 327 } 328 if (f == null || f.isZERO()) { // ? || f.isConstant() 329 return v; 330 } 331 RingFactory<C> cfac = f.ring.coFac; 332 C two = cfac.fromInteger(2); 333 334 while (true) { 335 long n = realRootCount(v, Sg); 336 logger.debug("n = " + n); 337 if (n == 0) { 338 return v; 339 } 340 C c = v.left.sum(v.right); 341 c = c.divide(two); 342 Interval<C> im = new Interval<C>(c, v.right); 343 if (signChange(im, f)) { 344 v = im; 345 } else { 346 v = new Interval<C>(v.left, c); 347 } 348 } 349 // return v; 350 } 351 352}