001/* 002 * $Id$ 003 */ 004 005package edu.jas.root; 006 007 008import java.util.ArrayList; 009// import java.util.Arrays; 010import java.util.List; 011 012import org.apache.logging.log4j.LogManager; 013import org.apache.logging.log4j.Logger; 014 015import edu.jas.arith.BigRational; 016import edu.jas.arith.Rational; 017import edu.jas.poly.Complex; 018import edu.jas.poly.ComplexRing; 019import edu.jas.poly.GenPolynomial; 020import edu.jas.poly.PolyUtil; 021import edu.jas.structure.RingElem; 022import edu.jas.structure.RingFactory; 023 024 025/** 026 * Complex roots implemented by Sturm sequences. Algorithms use exact method 027 * derived from Wilf's numeric Routh-Hurwitz method. 028 * @param <C> coefficient type. 029 * @author Heinz Kredel 030 */ 031public class ComplexRootsSturm<C extends RingElem<C> & Rational> extends ComplexRootsAbstract<C> { 032 033 034 private static final Logger logger = LogManager.getLogger(ComplexRootsSturm.class); 035 036 037 private static final boolean debug = logger.isDebugEnabled(); 038 039 040 /** 041 * Constructor. 042 * @param cf coefficient factory. 043 */ 044 public ComplexRootsSturm(RingFactory<Complex<C>> cf) { 045 super(cf); 046 //ufd = GCDFactory.<Complex<C>> getImplementation(cf); 047 } 048 049 050 /** 051 * Cauchy index of rational function f/g on interval. 052 * @param a interval bound for I = [a,b]. 053 * @param b interval bound for I = [a,b]. 054 * @param f univariate polynomial. 055 * @param g univariate polynomial. 056 * @return winding number of f/g in I. 057 */ 058 public long indexOfCauchy(C a, C b, GenPolynomial<C> f, GenPolynomial<C> g) { 059 List<GenPolynomial<C>> S = sturmSequence(g, f); 060 //System.out.println("S = " + S); 061 if (debug) { 062 logger.info("sturmSeq = {}", S); 063 } 064 RingFactory<C> cfac = f.ring.coFac; 065 List<C> l = PolyUtil.<C> evaluateMain(cfac, S, a); 066 List<C> r = PolyUtil.<C> evaluateMain(cfac, S, b); 067 long v = RootUtil.<C> signVar(l) - RootUtil.<C> signVar(r); 068 //System.out.println("v = " + v); 069 // if (v < 0L) { 070 // v = -v; 071 // } 072 return v; 073 } 074 075 076 /** 077 * Routh index of complex function f + i g on interval. 078 * @param a interval bound for I = [a,b]. 079 * @param b interval bound for I = [a,b]. 080 * @param f univariate polynomial. 081 * @param g univariate polynomial != 0. 082 * @return index number of f + i g. 083 */ 084 public long[] indexOfRouth(C a, C b, GenPolynomial<C> f, GenPolynomial<C> g) { 085 List<GenPolynomial<C>> S = sturmSequence(f, g); 086 //System.out.println("S = " + S); 087 RingFactory<C> cfac = f.ring.coFac; 088 List<C> l = PolyUtil.<C> evaluateMain(cfac, S, a); 089 List<C> r = PolyUtil.<C> evaluateMain(cfac, S, b); 090 long v = RootUtil.<C> signVar(l) - RootUtil.<C> signVar(r); 091 //System.out.println("v = " + v); 092 093 long d = f.degree(0); 094 if (d < g.degree(0)) { 095 d = g.degree(0); 096 } 097 //System.out.println("d = " + d); 098 long ui = (d - v) / 2; 099 long li = (d + v) / 2; 100 //System.out.println("upper = " + ui); 101 //System.out.println("lower = " + li); 102 return new long[] { ui, li }; 103 } 104 105 106 /** 107 * Sturm sequence. 108 * @param f univariate polynomial. 109 * @param g univariate polynomial. 110 * @return a Sturm sequence for f and g. 111 */ 112 public List<GenPolynomial<C>> sturmSequence(GenPolynomial<C> f, GenPolynomial<C> g) { 113 List<GenPolynomial<C>> S = new ArrayList<GenPolynomial<C>>(); 114 if (f == null || f.isZERO()) { 115 return S; 116 } 117 if (f.isConstant()) { 118 S.add(f.monic()); 119 return S; 120 } 121 GenPolynomial<C> F = f; 122 S.add(F); 123 GenPolynomial<C> G = g; //PolyUtil.<C> baseDerivative(f); 124 while (!G.isZERO()) { 125 GenPolynomial<C> r = F.remainder(G); 126 F = G; 127 G = r.negate(); 128 S.add(F/*.monic()*/); 129 } 130 //System.out.println("F = " + F); 131 if (F.isConstant()) { 132 return S; 133 } 134 // make squarefree 135 List<GenPolynomial<C>> Sp = new ArrayList<GenPolynomial<C>>(S.size()); 136 for (GenPolynomial<C> p : S) { 137 p = p.divide(F); 138 Sp.add(p); 139 } 140 return Sp; 141 } 142 143 144 /** 145 * Complex root count of complex polynomial on rectangle. 146 * @param rect rectangle. 147 * @param a univariate complex polynomial. 148 * @return root count of a in rectangle. 149 */ 150 @Override 151 public long complexRootCount(Rectangle<C> rect, GenPolynomial<Complex<C>> a) 152 throws InvalidBoundaryException { 153 C rl = rect.lengthReal(); 154 C il = rect.lengthImag(); 155 //System.out.println("complexRootCount: rl = " + rl + ", il = " + il); 156 // only linear polynomials have zero length intervals 157 if (rl.isZERO() && il.isZERO()) { 158 Complex<C> e = PolyUtil.<Complex<C>> evaluateMain(a.ring.coFac, a, rect.getSW()); 159 if (e.isZERO()) { 160 return 1; 161 } 162 return 0; 163 } 164 if (rl.isZERO() || il.isZERO()) { 165 //RingFactory<C> cf = (RingFactory<C>) rl.factory(); 166 //GenPolynomialRing<C> rfac = new GenPolynomialRing<C>(cf,a.ring); 167 //cf = (RingFactory<C>) il.factory(); 168 //GenPolynomialRing<C> ifac = new GenPolynomialRing<C>(cf,a.ring); 169 //GenPolynomial<C> rp = PolyUtil.<C> realPartFromComplex(rfac, a); 170 //GenPolynomial<C> ip = PolyUtil.<C> imaginaryPartFromComplex(ifac, a); 171 //RealRoots<C> rr = new RealRootsSturm<C>(); 172 if (rl.isZERO()) { 173 //logger.info("lengthReal == 0: {}", rect); 174 //Complex<C> r = rect.getSW(); 175 //r = new Complex<C>(r.ring,r.getRe()/*,0*/); 176 //Complex<C> e = PolyUtil.<Complex<C>> evaluateMain(a.ring.coFac, a, r); 177 //logger.info("a(re(rect)): {}", e); 178 //if ( !e.getRe().isZERO() ) { 179 // return 0; 180 //} 181 //C ev = PolyUtil.<C> evaluateMain(rp.ring.coFac, rp, rl); 182 //logger.info("re(a)(re(rect)): {}", ev); 183 //Interval<C> iv = new Interval<C>(rect.getSW().getIm(),rect.getNE().getIm()); 184 //logger.info("iv: {}", iv); 185 //long ic = rr.realRootCount(iv,ip); 186 //logger.info("ic: {}", ic); 187 188 Complex<C> sw = rect.getSW(); 189 Complex<C> ne = rect.getNE(); 190 C delta = sw.ring.ring.getONE(); //parse("1"); // works since linear polynomial 191 Complex<C> cd = new Complex<C>(sw.ring, delta/*, 0*/); 192 sw = sw.subtract(cd); 193 ne = ne.sum(cd); 194 rect = rect.exchangeSW(sw); 195 rect = rect.exchangeNE(ne); 196 logger.info("new rectangle: {}", rect.toScript()); 197 } 198 if (il.isZERO()) { 199 //logger.info("lengthImag == 0: {}", rect); 200 //Interval<C> rv = new Interval<C>(rect.getSW().getRe(),rect.getNE().getRe()); 201 //logger.info("rv: {}", rv); 202 //long rc = rr.realRootCount(rv,rp); 203 //logger.info("rc: {}", rc); 204 205 Complex<C> sw = rect.getSW(); 206 Complex<C> ne = rect.getNE(); 207 C delta = sw.ring.ring.getONE(); //parse("1"); // works since linear polynomial 208 Complex<C> cd = new Complex<C>(sw.ring, sw.ring.ring.getZERO(), delta); 209 sw = sw.subtract(cd); 210 ne = ne.sum(cd); 211 rect = rect.exchangeSW(sw); 212 rect = rect.exchangeNE(ne); 213 logger.info("new rectangle: {}", rect.toScript()); 214 } 215 } 216 long wn = windingNumber(rect, a); 217 //System.out.println("complexRootCount: wn = " + wn); 218 return wn; 219 } 220 221 222 /** 223 * Winding number of complex function A on rectangle. 224 * @param rect rectangle. 225 * @param A univariate complex polynomial. 226 * @return winding number of A around rect. 227 */ 228 public long windingNumber(Rectangle<C> rect, GenPolynomial<Complex<C>> A) 229 throws InvalidBoundaryException { 230 Boundary<C> bound = new Boundary<C>(rect, A); // throws InvalidBoundaryException 231 ComplexRing<C> cr = (ComplexRing<C>) A.ring.coFac; 232 RingFactory<C> cf = cr.ring; 233 C zero = cf.getZERO(); 234 C one = cf.getONE(); 235 long ix = 0L; 236 for (int i = 0; i < 4; i++) { 237 long ci = indexOfCauchy(zero, one, bound.getRealPart(i), bound.getImagPart(i)); 238 //System.out.println("ci[" + i + "," + (i + 1) + "] = " + ci); 239 ix += ci; 240 } 241 if (ix % 2L != 0) { 242 throw new InvalidBoundaryException("odd winding number " + ix); 243 } 244 return ix / 2L; 245 } 246 247 248 /** 249 * List of complex roots of complex polynomial a on rectangle. 250 * @param rect rectangle. 251 * @param a univariate squarefree complex polynomial. 252 * @return list of complex roots. 253 */ 254 @SuppressWarnings({ "cast", "unchecked" }) 255 @Override 256 public List<Rectangle<C>> complexRoots(Rectangle<C> rect, GenPolynomial<Complex<C>> a) 257 throws InvalidBoundaryException { 258 ComplexRing<C> cr = (ComplexRing<C>) a.ring.coFac; 259 List<Rectangle<C>> roots = new ArrayList<Rectangle<C>>(); 260 if (a.isConstant() || a.isZERO()) { 261 return roots; 262 } 263 //System.out.println("rect = " + rect); 264 long n = windingNumber(rect, a); 265 if (n < 0) { // can this happen? 266 throw new RuntimeException("negative winding number " + n); 267 //System.out.println("negative winding number " + n); 268 //return roots; 269 } 270 if (n == 0) { 271 return roots; 272 } 273 if (n == 1) { 274 //not ok: rect = excludeZero(rect, a); 275 roots.add(rect); 276 return roots; 277 } 278 Complex<C> eps = cr.fromInteger(1); 279 eps = eps.divide(cr.fromInteger(1000)); // 1/1000 280 //System.out.println("eps = " + eps); 281 //System.out.println("rect = " + rect); 282 // construct new center 283 Complex<C> delta = rect.corners[3].subtract(rect.corners[1]); 284 delta = delta.divide(cr.fromInteger(2)); 285 //System.out.println("delta = " + delta); 286 boolean work = true; 287 while (work) { 288 Complex<C> center = rect.corners[1].sum(delta); 289 //System.out.println("center = " + toDecimal(center)); 290 if (debug) { 291 logger.info("new center = {}", center); 292 } 293 try { 294 Complex<C>[] cp = (Complex<C>[]) copyOfComplex(rect.corners, 4); 295 // (Complex<C>[]) new Complex[4]; cp[0] = rect.corners[0]; 296 // cp[0] fix 297 cp[1] = new Complex<C>(cr, cp[1].getRe(), center.getIm()); 298 cp[2] = center; 299 cp[3] = new Complex<C>(cr, center.getRe(), cp[3].getIm()); 300 Rectangle<C> nw = new Rectangle<C>(cp); 301 //System.out.println("nw = " + nw); 302 List<Rectangle<C>> nwr = complexRoots(nw, a); 303 //System.out.println("#nwr = " + nwr.size()); 304 roots.addAll(nwr); 305 if (roots.size() == a.degree(0)) { 306 work = false; 307 break; 308 } 309 310 cp = (Complex<C>[]) copyOfComplex(rect.corners, 4); 311 cp[0] = new Complex<C>(cr, cp[0].getRe(), center.getIm()); 312 // cp[1] fix 313 cp[2] = new Complex<C>(cr, center.getRe(), cp[2].getIm()); 314 cp[3] = center; 315 Rectangle<C> sw = new Rectangle<C>(cp); 316 //System.out.println("sw = " + sw); 317 List<Rectangle<C>> swr = complexRoots(sw, a); 318 //System.out.println("#swr = " + swr.size()); 319 roots.addAll(swr); 320 if (roots.size() == a.degree(0)) { 321 work = false; 322 break; 323 } 324 325 cp = (Complex<C>[]) copyOfComplex(rect.corners, 4); 326 cp[0] = center; 327 cp[1] = new Complex<C>(cr, center.getRe(), cp[1].getIm()); 328 // cp[2] fix 329 cp[3] = new Complex<C>(cr, cp[3].getRe(), center.getIm()); 330 Rectangle<C> se = new Rectangle<C>(cp); 331 //System.out.println("se = " + se); 332 List<Rectangle<C>> ser = complexRoots(se, a); 333 //System.out.println("#ser = " + ser.size()); 334 roots.addAll(ser); 335 if (roots.size() == a.degree(0)) { 336 work = false; 337 break; 338 } 339 340 cp = (Complex<C>[]) copyOfComplex(rect.corners, 4); 341 cp[0] = new Complex<C>(cr, center.getRe(), cp[0].getIm()); 342 cp[1] = center; 343 cp[2] = new Complex<C>(cr, cp[2].getRe(), center.getIm()); 344 // cp[3] fix 345 Rectangle<C> ne = new Rectangle<C>(cp); 346 //System.out.println("ne = " + ne); 347 List<Rectangle<C>> ner = complexRoots(ne, a); 348 //System.out.println("#ner = " + ner.size()); 349 roots.addAll(ner); 350 work = false; 351 } catch (InvalidBoundaryException e) { 352 // repeat with new center 353 delta = delta.sum(delta.multiply(eps)); // distort 354 //System.out.println("new delta = " + toDecimal(delta)); 355 eps = eps.sum(eps.multiply(cr.getIMAG())); 356 } 357 } 358 return roots; 359 } 360 361 362 /** 363 * Invariant rectangle for algebraic number. 364 * @param rect root isolating rectangle for f which contains exactly one 365 * root. 366 * @param f univariate polynomial, non-zero. 367 * @param g univariate polynomial, gcd(f,g) == 1. 368 * @return v a new rectangle contained in rect such that g(w) != 0 for w in 369 * v. 370 */ 371 @Override 372 public Rectangle<C> invariantRectangle(Rectangle<C> rect, GenPolynomial<Complex<C>> f, 373 GenPolynomial<Complex<C>> g) throws InvalidBoundaryException { 374 Rectangle<C> v = rect; 375 if (g == null || g.isZERO()) { 376 return v; 377 } 378 if (g.isConstant()) { 379 return v; 380 } 381 if (f == null || f.isZERO() || f.isConstant()) { // ? 382 return v; 383 } 384 BigRational len = v.rationalLength(); 385 BigRational half = new BigRational(1, 2); 386 while (true) { 387 long n = windingNumber(v, g); 388 //System.out.println("n = " + n); 389 if (n < 0) { // can this happen? 390 throw new RuntimeException("negative winding number " + n); 391 } 392 if (n == 0) { 393 return v; 394 } 395 len = len.multiply(half); 396 Rectangle<C> v1 = v; 397 v = complexRootRefinement(v, f, len); 398 if (v.equals(v1)) { 399 //System.out.println("len = " + len); 400 if (!f.gcd(g).isONE()) { 401 System.out.println("f.gcd(g) = " + f.gcd(g)); 402 throw new RuntimeException("no convergence " + v); 403 } 404 //break; // no convergence 405 } 406 } 407 //return v; 408 } 409 410 411 /** 412 * Exclude zero. If an axis intersects with the rectangle, it is shrunk to 413 * exclude the axis. Not used. 414 * @param rect root isolating rectangle for f which contains exactly one 415 * root. 416 * @return a new rectangle r such that re(r) < 0 or (re)r > 0 and 417 * im(r) < 0 or (im)r > 0. 418 */ 419 public Rectangle<C> excludeZero(Rectangle<C> rect, GenPolynomial<Complex<C>> f) 420 throws InvalidBoundaryException { 421 if (f == null || f.isZERO()) { 422 return rect; 423 } 424 System.out.println("\nexcludeZero: rect = " + rect + ", f = " + f); 425 Complex<C> zero = f.ring.coFac.getZERO(); 426 ComplexRing<C> cr = zero.ring; 427 Complex<C> sw = rect.getSW(); 428 Complex<C> ne = rect.getNE(); 429 Interval<C> ir = new Interval<C>(sw.getRe(), ne.getRe()); 430 Interval<C> ii = new Interval<C>(sw.getIm(), ne.getIm()); 431 System.out.println("intervals, ir = " + ir + ", ii = " + ii); 432 if (!(ir.contains(zero.getRe()) || ii.contains(zero.getIm()))) { 433 // !rect.contains(zero) not correct 434 return rect; 435 } 436 //System.out.println("contains: ir = " + ir.contains(zero.getRe()) + ", ii = " + ii.contains(zero.getIm()) ); 437 Rectangle<C> rn = rect; 438 // shrink real part 439 if (ir.contains(zero.getRe())) { 440 Complex<C> sw0 = new Complex<C>(cr, zero.getRe(), sw.getIm()); 441 Complex<C> ne0 = new Complex<C>(cr, zero.getRe(), ne.getIm()); 442 Rectangle<C> rl = new Rectangle<C>(sw, ne0); 443 Rectangle<C> rr = new Rectangle<C>(sw0, ne); 444 System.out.println("rectangle, rl = " + rl + ", rr = " + rr); 445 if (complexRootCount(rr, f) == 1) { 446 rn = rr; 447 } else { // complexRootCount(rl,f) == 1 448 rn = rl; 449 } 450 System.out.println("rectangle, real = " + rn); 451 } 452 // shrink imaginary part 453 sw = rn.getSW(); 454 ne = rn.getNE(); 455 ii = new Interval<C>(sw.getIm(), ne.getIm()); 456 System.out.println("interval, ii = " + ii); 457 if (ii.contains(zero.getIm())) { 458 Complex<C> sw1 = new Complex<C>(cr, sw.getRe(), zero.getIm()); 459 Complex<C> ne1 = new Complex<C>(cr, ne.getRe(), zero.getIm()); 460 Rectangle<C> iu = new Rectangle<C>(sw1, ne); 461 Rectangle<C> il = new Rectangle<C>(sw, ne1); 462 System.out.println("rectangle, il = " + il + ", iu = " + iu); 463 if (complexRootCount(il, f) == 1) { 464 rn = il;; 465 } else { // complexRootCount(iu,f) == 1 466 rn = iu; 467 } 468 System.out.println("rectangle, imag = " + rn); 469 } 470 return rn; 471 } 472 473}