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