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