001/* 002 * $Id: ComplexRootsAbstract.java 5872 2018-07-20 16:01:46Z kredel $ 003 */ 004 005package edu.jas.root; 006 007 008import java.util.ArrayList; 009import java.util.List; 010import java.util.Map; 011import java.util.SortedMap; 012 013import org.apache.logging.log4j.Logger; 014import org.apache.logging.log4j.LogManager; 015 016import edu.jas.arith.BigDecimal; 017import edu.jas.arith.BigRational; 018import edu.jas.arith.Rational; 019import edu.jas.poly.Complex; 020import edu.jas.poly.ComplexRing; 021import edu.jas.poly.GenPolynomial; 022import edu.jas.poly.GenPolynomialRing; 023import edu.jas.poly.PolyUtil; 024import edu.jas.structure.RingElem; 025import edu.jas.structure.RingFactory; 026import edu.jas.structure.UnaryFunctor; 027import edu.jas.ufd.Squarefree; 028import edu.jas.ufd.SquarefreeFactory; 029 030 031/** 032 * Complex roots abstract class. 033 * @param <C> coefficient type. 034 * @author Heinz Kredel 035 */ 036public abstract class ComplexRootsAbstract<C extends RingElem<C> & Rational> implements ComplexRoots<C> { 037 038 039 private static final Logger logger = LogManager.getLogger(ComplexRootsAbstract.class); 040 041 042 private static final boolean debug = logger.isDebugEnabled(); 043 044 045 /** 046 * Engine for square free decomposition. 047 */ 048 public final Squarefree<Complex<C>> engine; 049 050 051 /** 052 * Constructor. 053 * @param cf coefficient factory. 054 */ 055 public ComplexRootsAbstract(RingFactory<Complex<C>> cf) { 056 if (!(cf instanceof ComplexRing)) { 057 throw new IllegalArgumentException("cf not supported coefficients " + cf); 058 } 059 engine = SquarefreeFactory.<Complex<C>> getImplementation(cf); 060 } 061 062 063 /** 064 * Root bound. With f(-M + i M) * f(-M - i M) * f(M - i M) * f(M + i M) != 065 * 0. 066 * @param f univariate polynomial. 067 * @return M such that root(f) is contained in the rectangle spanned by M. 068 */ 069 public Complex<C> rootBound(GenPolynomial<Complex<C>> f) { 070 if (f == null) { 071 return null; 072 } 073 RingFactory<Complex<C>> cfac = f.ring.coFac; 074 Complex<C> M = cfac.getONE(); 075 if (f.isZERO() || f.isConstant()) { 076 return M; 077 } 078 Complex<C> a = f.leadingBaseCoefficient().norm(); 079 for (Complex<C> c : f.getMap().values()) { 080 Complex<C> d = c.norm().divide(a); 081 if (M.compareTo(d) < 0) { 082 M = d; 083 } 084 } 085 M = M.sum(cfac.getONE()); 086 //System.out.println("M = " + M); 087 return M; 088 } 089 090 091 /** 092 * Magnitude bound. 093 * @param rect rectangle. 094 * @param f univariate polynomial. 095 * @return B such that |f(c)| < B for c in rect. 096 */ 097 public C magnitudeBound(Rectangle<C> rect, GenPolynomial<Complex<C>> f) { 098 if (f == null) { 099 return null; 100 } 101 if (f.isZERO()) { 102 return f.ring.coFac.getONE().getRe(); 103 } 104 //System.out.println("f = " + f); 105 if (f.isConstant()) { 106 Complex<C> c = f.leadingBaseCoefficient(); 107 return c.norm().getRe(); 108 } 109 GenPolynomial<Complex<C>> fa = f.map(new UnaryFunctor<Complex<C>, Complex<C>>() { 110 111 112 public Complex<C> eval(Complex<C> a) { 113 return a.norm(); 114 } 115 }); 116 //System.out.println("fa = " + fa); 117 Complex<C> Mc = rect.getNW().norm(); 118 C M = Mc.getRe(); 119 //System.out.println("M = " + M); 120 Complex<C> M1c = rect.getSW().norm(); 121 C M1 = M1c.getRe(); 122 if (M.compareTo(M1) < 0) { 123 M = M1; 124 Mc = M1c; 125 } 126 M1c = rect.getSE().norm(); 127 M1 = M1c.getRe(); 128 if (M.compareTo(M1) < 0) { 129 M = M1; 130 Mc = M1c; 131 } 132 M1c = rect.getNE().norm(); 133 M1 = M1c.getRe(); 134 if (M.compareTo(M1) < 0) { 135 //M = M1; 136 Mc = M1c; 137 } 138 //System.out.println("M = " + M); 139 Complex<C> B = PolyUtil.<Complex<C>> evaluateMain(f.ring.coFac, fa, Mc); 140 //System.out.println("B = " + B); 141 return B.getRe(); 142 } 143 144 145 /** 146 * Complex root count of complex polynomial on rectangle. 147 * @param rect rectangle. 148 * @param a univariate complex polynomial. 149 * @return root count of a in rectangle. 150 */ 151 public abstract long complexRootCount(Rectangle<C> rect, GenPolynomial<Complex<C>> a) 152 throws InvalidBoundaryException; 153 154 155 /** 156 * List of complex roots of complex polynomial a on rectangle. 157 * @param rect rectangle. 158 * @param a univariate squarefree complex polynomial. 159 * @return list of complex roots. 160 */ 161 public abstract List<Rectangle<C>> complexRoots(Rectangle<C> rect, GenPolynomial<Complex<C>> a) 162 throws InvalidBoundaryException; 163 164 165 /** 166 * List of complex roots of complex polynomial. 167 * @param a univariate complex polynomial. 168 * @return list of complex roots. 169 */ 170 @SuppressWarnings({"cast","unchecked"}) 171 public List<Rectangle<C>> complexRoots(GenPolynomial<Complex<C>> a) { 172 List<Rectangle<C>> roots = new ArrayList<Rectangle<C>>(); 173 if (a.isConstant() || a.isZERO()) { 174 return roots; 175 } 176 ComplexRing<C> cr = (ComplexRing<C>) a.ring.coFac; 177 SortedMap<GenPolynomial<Complex<C>>, Long> sa = engine.squarefreeFactors(a); 178 for (Map.Entry<GenPolynomial<Complex<C>>, Long> me : sa.entrySet()) { 179 GenPolynomial<Complex<C>> p = me.getKey(); 180 Complex<C> Mb = rootBound(p); 181 C M = Mb.getRe(); 182 C M1 = M.sum(M.factory().fromInteger(1)); // asymmetric to origin 183 //System.out.println("M = " + M); 184 if (debug) { 185 logger.info("rootBound = " + M); 186 } 187 Complex<C>[] corner = (Complex<C>[]) new Complex[4]; 188 corner[0] = new Complex<C>(cr, M1.negate(), M); // nw 189 corner[1] = new Complex<C>(cr, M1.negate(), M1.negate()); // sw 190 corner[2] = new Complex<C>(cr, M, M1.negate()); // se 191 corner[3] = new Complex<C>(cr, M, M); // ne 192 Rectangle<C> rect = new Rectangle<C>(corner); 193 try { 194 List<Rectangle<C>> rs = complexRoots(rect, p); 195 long e = me.getValue(); // sa.get(p); 196 for (int i = 0; i < e; i++) { // add with multiplicity 197 roots.addAll(rs); 198 } 199 } catch (InvalidBoundaryException e) { 200 //logger.error("invalid boundary for p = " + p); 201 throw new RuntimeException("this should never happen " + e); 202 } 203 } 204 return roots; 205 } 206 207 208 /** 209 * Complex root refinement of complex polynomial a on rectangle. 210 * @param rect rectangle containing exactly one complex root. 211 * @param a univariate squarefree complex polynomial. 212 * @param len rational length for refinement. 213 * @return refined complex root. 214 */ 215 @SuppressWarnings({"cast","unchecked"}) 216 public Rectangle<C> complexRootRefinement(Rectangle<C> rect, GenPolynomial<Complex<C>> a, BigRational len) 217 throws InvalidBoundaryException { 218 ComplexRing<C> cr = (ComplexRing<C>) a.ring.coFac; 219 Rectangle<C> root = rect; 220 long w; 221 if (debug) { 222 w = complexRootCount(root, a); 223 if (w != 1) { 224 System.out.println("#root = " + w); 225 System.out.println("root = " + root); 226 throw new ArithmeticException("no initial isolating rectangle " + rect); 227 } 228 } 229 Complex<C> eps = cr.fromInteger(1); 230 eps = eps.divide(cr.fromInteger(1000)); // 1/1000 231 BigRational length = len.multiply(len); 232 Complex<C> delta = null; 233 boolean work = true; 234 while (work) { 235 try { 236 while (root.rationalLength().compareTo(length) > 0) { 237 //System.out.println("root = " + root + ", len = " + new BigDecimal(root.rationalLength())); 238 if (delta == null) { 239 delta = root.corners[3].subtract(root.corners[1]); 240 delta = delta.divide(cr.fromInteger(2)); 241 //System.out.println("delta = " + toDecimal(delta)); 242 } 243 Complex<C> center = root.corners[1].sum(delta); 244 //System.out.println("refine center = " + toDecimal(center)); 245 if (debug) { 246 logger.info("new center = " + center); 247 } 248 249 Complex<C>[] cp = (Complex<C>[]) copyOfComplex(root.corners, 4); 250 // cp[0] fix 251 cp[1] = new Complex<C>(cr, cp[1].getRe(), center.getIm()); 252 cp[2] = center; 253 cp[3] = new Complex<C>(cr, center.getRe(), cp[3].getIm()); 254 Rectangle<C> nw = new Rectangle<C>(cp); 255 w = complexRootCount(nw, a); 256 if (w == 1) { 257 root = nw; 258 delta = null; 259 continue; 260 } 261 262 cp = (Complex<C>[]) copyOfComplex(root.corners, 4); 263 cp[0] = new Complex<C>(cr, cp[0].getRe(), center.getIm()); 264 // cp[1] fix 265 cp[2] = new Complex<C>(cr, center.getRe(), cp[2].getIm()); 266 cp[3] = center; 267 Rectangle<C> sw = new Rectangle<C>(cp); 268 w = complexRootCount(sw, a); 269 //System.out.println("#swr = " + w); 270 if (w == 1) { 271 root = sw; 272 delta = null; 273 continue; 274 } 275 276 cp = (Complex<C>[]) copyOfComplex(root.corners, 4); 277 cp[0] = center; 278 cp[1] = new Complex<C>(cr, center.getRe(), cp[1].getIm()); 279 // cp[2] fix 280 cp[3] = new Complex<C>(cr, cp[3].getRe(), center.getIm()); 281 Rectangle<C> se = new Rectangle<C>(cp); 282 w = complexRootCount(se, a); 283 //System.out.println("#ser = " + w); 284 if (w == 1) { 285 root = se; 286 delta = null; 287 continue; 288 } 289 290 cp = (Complex<C>[]) copyOfComplex(root.corners, 4); 291 cp[0] = new Complex<C>(cr, center.getRe(), cp[0].getIm()); 292 cp[1] = center; 293 cp[2] = new Complex<C>(cr, cp[2].getRe(), center.getIm()); 294 // cp[3] fix 295 Rectangle<C> ne = new Rectangle<C>(cp); 296 w = complexRootCount(ne, a); 297 //System.out.println("#ner = " + w); 298 if (w == 1) { 299 root = ne; 300 delta = null; 301 continue; 302 } 303 if (true) { 304 w = complexRootCount(root, a); 305 System.out.println("#root = " + w); 306 System.out.println("root = " + root); 307 } 308 throw new ArithmeticException("no isolating rectangle " + rect); 309 } 310 work = false; 311 } catch (InvalidBoundaryException e) { 312 // repeat with new center 313 delta = delta.sum(delta.multiply(eps)); // distort 314 //System.out.println("new refine delta = " + toDecimal(delta)); 315 eps = eps.sum(eps.multiply(cr.getIMAG())); 316 } 317 } 318 return root; 319 } 320 321 322 /** 323 * List of complex roots of complex polynomial. 324 * @param a univariate complex polynomial. 325 * @param len rational length for refinement. 326 * @return list of complex roots to desired precision. 327 */ 328 @SuppressWarnings({"cast","unchecked"}) 329 public List<Rectangle<C>> complexRoots(GenPolynomial<Complex<C>> a, BigRational len) { 330 ComplexRing<C> cr = (ComplexRing<C>) a.ring.coFac; 331 SortedMap<GenPolynomial<Complex<C>>, Long> sa = engine.squarefreeFactors(a); 332 List<Rectangle<C>> roots = new ArrayList<Rectangle<C>>(); 333 for (Map.Entry<GenPolynomial<Complex<C>>, Long> me : sa.entrySet()) { 334 GenPolynomial<Complex<C>> p = me.getKey(); 335 Complex<C> Mb = rootBound(p); 336 C M = Mb.getRe(); 337 C M1 = M.sum(M.factory().fromInteger(1)); // asymmetric to origin 338 if (debug) { 339 logger.info("rootBound = " + M); 340 } 341 Complex<C>[] corner = (Complex<C>[]) new Complex[4]; 342 corner[0] = new Complex<C>(cr, M1.negate(), M); // nw 343 corner[1] = new Complex<C>(cr, M1.negate(), M1.negate()); // sw 344 corner[2] = new Complex<C>(cr, M, M1.negate()); // se 345 corner[3] = new Complex<C>(cr, M, M); // ne 346 Rectangle<C> rect = new Rectangle<C>(corner); 347 try { 348 List<Rectangle<C>> rs = complexRoots(rect, p); 349 List<Rectangle<C>> rf = new ArrayList<Rectangle<C>>(rs.size()); 350 for (Rectangle<C> r : rs) { 351 Rectangle<C> rr = complexRootRefinement(r, p, len); 352 rf.add(rr); 353 } 354 long e = me.getValue(); // sa.get(p); 355 for (int i = 0; i < e; i++) { // add with multiplicity 356 roots.addAll(rf); 357 } 358 } catch (InvalidBoundaryException e) { 359 throw new RuntimeException("this should never happen " + e); 360 } 361 } 362 return roots; 363 } 364 365 366 /** 367 * Invariant rectangle for algebraic number. 368 * @param rect root isolating rectangle for f which contains exactly one 369 * root. 370 * @param f univariate polynomial, non-zero. 371 * @param g univariate polynomial, gcd(f,g) == 1. 372 * @return v with v a new rectangle contained in iv such that g(w) != 0 for 373 * w in v. 374 */ 375 public abstract Rectangle<C> invariantRectangle(Rectangle<C> rect, GenPolynomial<Complex<C>> f, 376 GenPolynomial<Complex<C>> g) throws InvalidBoundaryException; 377 378 379 /** 380 * Get decimal approximation. 381 * @param a complex number. 382 * @return decimal(a). 383 */ 384 public String toDecimal(Complex<C> a) { 385 C r = a.getRe(); 386 String s = r.toString(); 387 BigRational rs = new BigRational(s); 388 BigDecimal rd = new BigDecimal(rs); 389 C i = a.getIm(); 390 s = i.toString(); 391 BigRational is = new BigRational(s); 392 BigDecimal id = new BigDecimal(is); 393 //System.out.println("rd = " + rd); 394 //System.out.println("id = " + id); 395 return rd.toString() + " i " + id.toString(); 396 } 397 398 399 /** 400 * Approximate complex root. 401 * @param rt root isolating rectangle. 402 * @param f univariate polynomial, non-zero. 403 * @param eps requested interval length. 404 * @return a decimal approximation d such that |d-v| < eps, for f(v) = 0, 405 * v in rt. 406 */ 407 public Complex<BigDecimal> approximateRoot(Rectangle<C> rt, GenPolynomial<Complex<C>> f, BigRational eps) 408 throws NoConvergenceException { 409 if (rt == null) { 410 throw new IllegalArgumentException("null interval not allowed"); 411 } 412 Complex<BigDecimal> d = rt.getDecimalCenter(); 413 //System.out.println("d = " + d); 414 if (f == null || f.isZERO() || f.isConstant() || eps == null) { 415 return d; 416 } 417 if (rt.rationalLength().compareTo(eps) < 0) { 418 return d; 419 } 420 ComplexRing<BigDecimal> cr = d.ring; 421 Complex<C> sw = rt.getSW(); 422 BigDecimal swr = new BigDecimal(sw.getRe().getRational()); 423 BigDecimal swi = new BigDecimal(sw.getIm().getRational()); 424 Complex<BigDecimal> ll = new Complex<BigDecimal>(cr, swr, swi); 425 Complex<C> ne = rt.getNE(); 426 BigDecimal ner = new BigDecimal(ne.getRe().getRational()); 427 BigDecimal nei = new BigDecimal(ne.getIm().getRational()); 428 Complex<BigDecimal> ur = new Complex<BigDecimal>(cr, ner, nei); 429 430 BigDecimal e = new BigDecimal(eps.getRational()); 431 Complex<BigDecimal> q = new Complex<BigDecimal>(cr, new BigDecimal("0.25")); 432 e = e.multiply(d.norm().getRe()); // relative error 433 //System.out.println("e = " + e); 434 435 // polynomials with decimal coefficients 436 GenPolynomialRing<Complex<BigDecimal>> dfac = new GenPolynomialRing<Complex<BigDecimal>>(cr, f.ring); 437 GenPolynomial<Complex<BigDecimal>> df = PolyUtil.<C> complexDecimalFromRational(dfac, f); 438 GenPolynomial<Complex<C>> fp = PolyUtil.<Complex<C>> baseDeriviative(f); 439 GenPolynomial<Complex<BigDecimal>> dfp = PolyUtil.<C> complexDecimalFromRational(dfac, fp); 440 441 // Newton Raphson iteration: x_{n+1} = x_n - f(x_n)/f'(x_n) 442 int i = 0; 443 final int MITER = 50; 444 int dir = -1; 445 while (i++ < MITER) { 446 Complex<BigDecimal> fx = PolyUtil.<Complex<BigDecimal>> evaluateMain(cr, df, d); // f(d) 447 //BigDecimal fs = fx.norm().getRe(); 448 //System.out.println("fs = " + fs); 449 if (fx.isZERO()) { 450 return d; 451 } 452 Complex<BigDecimal> fpx = PolyUtil.<Complex<BigDecimal>> evaluateMain(cr, dfp, d); // f'(d) 453 if (fpx.isZERO()) { 454 throw new NoConvergenceException("zero deriviative should not happen"); 455 } 456 Complex<BigDecimal> x = fx.divide(fpx); 457 Complex<BigDecimal> dx = d.subtract(x); 458 //System.out.println("dx = " + dx); 459 if (d.subtract(dx).norm().getRe().compareTo(e) <= 0) { 460 return dx; 461 } 462 // if ( false ) { // not useful: 463 // Complex<BigDecimal> fxx = PolyUtil.<Complex<BigDecimal>> evaluateMain(cr, df, dx); // f(dx) 464 // //System.out.println("fxx = " + fxx); 465 // BigDecimal fsx = fxx.norm().getRe(); 466 // System.out.println("fsx = " + fsx); 467 // while ( fsx.compareTo( fs ) >= 0 ) { 468 // System.out.println("trying to increase f(d) "); 469 // if ( i++ > MITER ) { // dx > right: dx - right > 0 470 // throw new NoConvergenceException("no convergence after " + i + " steps"); 471 // } 472 // x = x.multiply(q); // x * 1/4 473 // dx = d.subtract(x); 474 // //System.out.println(" x = " + x); 475 // System.out.println("dx = " + dx); 476 // fxx = PolyUtil.<Complex<BigDecimal>> evaluateMain(cr, df, dx); // f(dx) 477 // //System.out.println("fxx = " + fxx); 478 // fsx = fxx.norm().getRe(); 479 // System.out.println("fsx = " + fsx); 480 // } 481 // } 482 // check interval bounds 483 while (dx.getRe().compareTo(ll.getRe()) < 0 || dx.getIm().compareTo(ll.getIm()) < 0 484 || dx.getRe().compareTo(ur.getRe()) > 0 || dx.getIm().compareTo(ur.getIm()) > 0) { 485 // dx < ll: dx - ll < 0 486 // dx > ur: dx - ur > 0 487 if (i++ > MITER) { // dx > right: dx - right > 0 488 throw new NoConvergenceException("no convergence after " + i + " steps"); 489 } 490 if (i > MITER / 2 && dir == 0) { 491 Complex<C> cc = rt.getCenter(); 492 Rectangle<C> nrt = rt.exchangeSE(cc); 493 Complex<BigDecimal> sd = nrt.getDecimalCenter(); 494 d = sd; 495 x = cr.getZERO(); 496 logger.info("trying new SE starting point " + d); 497 i = 0; 498 dir = 1; 499 } 500 if (i > MITER / 2 && dir == 1) { 501 Complex<C> cc = rt.getCenter(); 502 Rectangle<C> nrt = rt.exchangeNW(cc); 503 Complex<BigDecimal> sd = nrt.getDecimalCenter(); 504 d = sd; 505 x = cr.getZERO(); 506 logger.info("trying new NW starting point " + d); 507 i = 0; 508 dir = 2; 509 } 510 if (i > MITER / 2 && dir == 2) { 511 Complex<C> cc = rt.getCenter(); 512 Rectangle<C> nrt = rt.exchangeSW(cc); 513 Complex<BigDecimal> sd = nrt.getDecimalCenter(); 514 d = sd; 515 x = cr.getZERO(); 516 logger.info("trying new SW starting point " + d); 517 i = 0; 518 dir = 3; 519 } 520 if (i > MITER / 2 && dir == 3) { 521 Complex<C> cc = rt.getCenter(); 522 Rectangle<C> nrt = rt.exchangeNE(cc); 523 Complex<BigDecimal> sd = nrt.getDecimalCenter(); 524 d = sd; 525 x = cr.getZERO(); 526 logger.info("trying new NE starting point " + d); 527 i = 0; 528 dir = 4; 529 } 530 if (i > MITER / 2 && (dir == -1 || dir == 4 || dir == 5)) { 531 Complex<C> sr = rt.randomPoint(); 532 BigDecimal srr = new BigDecimal(sr.getRe().getRational()); 533 BigDecimal sri = new BigDecimal(sr.getIm().getRational()); 534 Complex<BigDecimal> sd = new Complex<BigDecimal>(cr, srr, sri); 535 d = sd; 536 x = cr.getZERO(); 537 logger.info("trying new random starting point " + d); 538 if (dir == -1) { 539 i = 0; 540 dir = 0; 541 } else if (dir == 4) { 542 i = 0; 543 dir = 5; 544 } else { 545 //i = 0; 546 dir = 6; // end 547 } 548 } 549 x = x.multiply(q); // x * 1/4 550 dx = d.subtract(x); 551 //System.out.println(" x = " + x); 552 //System.out.println("dx = " + dx); 553 } 554 d = dx; 555 } 556 throw new NoConvergenceException("no convergence after " + i + " steps"); 557 } 558 559 560 /** 561 * List of decimal approximations of complex roots of complex polynomial. 562 * @param a univariate complex polynomial. 563 * @param eps length for refinement. 564 * @return list of complex decimal root approximations to desired precision. 565 */ 566 @SuppressWarnings({"cast","unchecked"}) 567 public List<Complex<BigDecimal>> approximateRoots(GenPolynomial<Complex<C>> a, BigRational eps) { 568 ComplexRing<C> cr = (ComplexRing<C>) a.ring.coFac; 569 SortedMap<GenPolynomial<Complex<C>>, Long> sa = engine.squarefreeFactors(a); 570 List<Complex<BigDecimal>> roots = new ArrayList<Complex<BigDecimal>>(); 571 for (Map.Entry<GenPolynomial<Complex<C>>, Long> me : sa.entrySet()) { 572 GenPolynomial<Complex<C>> p = me.getKey(); 573 List<Complex<BigDecimal>> rf = null; 574 if (p.degree(0) <= 1) { 575 Complex<C> tc = p.trailingBaseCoefficient(); 576 tc = tc.negate(); 577 BigDecimal rr = new BigDecimal(tc.getRe().getRational()); 578 BigDecimal ri = new BigDecimal(tc.getIm().getRational()); 579 ComplexRing<BigDecimal> crf = new ComplexRing<BigDecimal>(rr); 580 Complex<BigDecimal> r = new Complex<BigDecimal>(crf, rr, ri); 581 rf = new ArrayList<Complex<BigDecimal>>(1); 582 rf.add(r); 583 } else { 584 Complex<C> Mb = rootBound(p); 585 C M = Mb.getRe(); 586 C M1 = M.sum(M.factory().fromInteger(1)); // asymmetric to origin 587 if (debug) { 588 logger.info("rootBound = " + M); 589 } 590 Complex<C>[] corner = (Complex<C>[]) new Complex[4]; 591 corner[0] = new Complex<C>(cr, M1.negate(), M); // nw 592 corner[1] = new Complex<C>(cr, M1.negate(), M1.negate()); // sw 593 corner[2] = new Complex<C>(cr, M, M1.negate()); // se 594 corner[3] = new Complex<C>(cr, M, M); // ne 595 Rectangle<C> rect = new Rectangle<C>(corner); 596 List<Rectangle<C>> rs = null; 597 try { 598 rs = complexRoots(rect, p); 599 } catch (InvalidBoundaryException e) { 600 throw new RuntimeException("this should never happen " + e); 601 } 602 rf = new ArrayList<Complex<BigDecimal>>(rs.size()); 603 for (Rectangle<C> r : rs) { 604 Complex<BigDecimal> rr = null; 605 while (rr == null) { 606 try { 607 rr = approximateRoot(r, p, eps); 608 rf.add(rr); 609 } catch (NoConvergenceException e) { 610 // fall back to exact algorithm 611 BigRational len = r.rationalLength(); 612 len = len.multiply(new BigRational(1, 1000)); 613 try { 614 r = complexRootRefinement(r, p, len); 615 logger.info("fall back rootRefinement = " + r); 616 //System.out.println("len = " + len); 617 } catch (InvalidBoundaryException ee) { 618 throw new RuntimeException("this should never happen " + ee); 619 } 620 } 621 } 622 } 623 } 624 long e = me.getValue(); // sa.get(p); 625 for (int i = 0; i < e; i++) { // add with multiplicity 626 roots.addAll(rf); 627 } 628 } 629 return roots; 630 } 631 632 633 /** 634 * Copy the specified array. 635 * @param original array. 636 * @param newLength new array length. 637 * @return copy of this. 638 */ 639 public Complex[] copyOfComplex(Complex[] original, int newLength) { 640 Complex[] copy = new Complex[newLength]; 641 System.arraycopy(original, 0, copy, 0, Math.min(original.length, newLength)); 642 return copy; 643 } 644 645 646 /** 647 * Invariant rectangle for algebraic number magnitude. 648 * @param rect root isolating rectangle for f which contains exactly one 649 * root. 650 * @param f univariate polynomial, non-zero. 651 * @param g univariate polynomial, gcd(f,g) == 1. 652 * @param eps length limit for rectangle length. 653 * @return v with v a new rectangle contained in rect such that |g(a) - 654 * g(b)| < eps for a, b in v in rect. 655 */ 656 public Rectangle<C> invariantMagnitudeRectangle(Rectangle<C> rect, GenPolynomial<Complex<C>> f, 657 GenPolynomial<Complex<C>> g, BigRational eps) throws InvalidBoundaryException { 658 Rectangle<C> v = rect; 659 if (g == null || g.isZERO()) { 660 return v; 661 } 662 if (g.isConstant()) { 663 return v; 664 } 665 if (f == null || f.isZERO() || f.isConstant()) { // ? 666 return v; 667 } 668 GenPolynomial<Complex<C>> gp = PolyUtil.<Complex<C>> baseDeriviative(g); 669 //System.out.println("g = " + g); 670 //System.out.println("gp = " + gp); 671 BigRational B = magnitudeBound(rect, gp).getRational(); 672 //System.out.println("B = " + B + " : " + B.getClass()); 673 674 BigRational len = v.rationalLength(); 675 BigRational half = new BigRational(1, 2); 676 677 BigRational vlen = v.rationalLength(); 678 vlen = vlen.multiply(vlen); 679 //eps = eps.multiply(eps); 680 //System.out.println("v = " + v); 681 //System.out.println("vlen = " + vlen); 682 while (B.multiply(vlen).compareTo(eps) >= 0) { // TODO: test squared 683 len = len.multiply(half); 684 v = complexRootRefinement(v, f, len); 685 //System.out.println("v = " + v); 686 vlen = v.rationalLength(); 687 vlen = vlen.multiply(vlen); 688 //System.out.println("vlen = " + vlen); 689 } 690 //System.out.println("vlen = " + vlen); 691 return v; 692 } 693 694 695 /** 696 * Complex algebraic number magnitude. 697 * @param rect root isolating rectangle for f which contains exactly one 698 * root, with rect such that |g(a) - g(b)| < eps for a, b in 699 * rect. 700 * @param f univariate polynomial, non-zero. 701 * @param g univariate polynomial, gcd(f,g) == 1. 702 * @return g(rect) . 703 */ 704 public Complex<C> complexRectangleMagnitude(Rectangle<C> rect, GenPolynomial<Complex<C>> f, 705 GenPolynomial<Complex<C>> g) { 706 if (g.isZERO() || g.isConstant()) { 707 return g.leadingBaseCoefficient(); 708 } 709 RingFactory<Complex<C>> cfac = f.ring.coFac; 710 //System.out.println("cfac = " + cfac + " : " + cfac.getClass()); 711 Complex<C> c = rect.getCenter(); 712 Complex<C> ev = PolyUtil.<Complex<C>> evaluateMain(cfac, g, c); 713 return ev; 714 } 715 716 717 /** 718 * Complex algebraic number magnitude. 719 * @param rect root isolating rectangle for f which contains exactly one 720 * root, with rect such that |g(a) - g(b)| < eps for a, b in 721 * rect. 722 * @param f univariate polynomial, non-zero. 723 * @param g univariate polynomial, gcd(f,g) == 1. 724 * @param eps length limit for rectangle length. 725 * @return g(rect) . 726 */ 727 public Complex<C> complexMagnitude(Rectangle<C> rect, GenPolynomial<Complex<C>> f, 728 GenPolynomial<Complex<C>> g, BigRational eps) throws InvalidBoundaryException { 729 if (g.isZERO() || g.isConstant()) { 730 return g.leadingBaseCoefficient(); 731 } 732 Rectangle<C> v = invariantMagnitudeRectangle(rect, f, g, eps); 733 //System.out.println("ref = " + ref); 734 return complexRectangleMagnitude(v, f, g); 735 } 736 737}