001 /*
002 * $Id: ComplexRootTest.java 3563 2011-03-12 12:15:40Z kredel $
003 */
004
005 package edu.jas.root;
006
007
008 import java.util.ArrayList;
009 import java.util.Collections;
010 import java.util.List;
011
012 import junit.framework.Test;
013 import junit.framework.TestCase;
014 import junit.framework.TestSuite;
015
016 import org.apache.log4j.BasicConfigurator;
017 //import org.apache.log4j.Logger;
018
019 import edu.jas.arith.BigDecimal;
020 import edu.jas.arith.BigRational;
021 import edu.jas.kern.ComputerThreads;
022 import edu.jas.poly.Complex;
023 import edu.jas.poly.ComplexRing;
024 import edu.jas.poly.GenPolynomial;
025 import edu.jas.poly.GenPolynomialRing;
026 import edu.jas.poly.PolyUtil;
027 import edu.jas.poly.TermOrder;
028 import edu.jas.structure.Power;
029 import edu.jas.ufd.Squarefree;
030 import edu.jas.ufd.SquarefreeFactory;
031
032
033 /**
034 * RootUtil tests with JUnit.
035 * @author Heinz Kredel.
036 */
037
038 public class ComplexRootTest extends TestCase {
039
040
041 /**
042 * main.
043 */
044 public static void main(String[] args) {
045 //BasicConfigurator.configure();
046 junit.textui.TestRunner.run(suite());
047 ComputerThreads.terminate();
048 }
049
050
051 /**
052 * Constructs a <CODE>ComplexRootTest</CODE> object.
053 * @param name String.
054 */
055 public ComplexRootTest(String name) {
056 super(name);
057 }
058
059
060 /**
061 */
062 public static Test suite() {
063 TestSuite suite = new TestSuite(ComplexRootTest.class);
064 return suite;
065 }
066
067
068 TermOrder to = new TermOrder(TermOrder.INVLEX);
069
070
071 GenPolynomialRing<Complex<BigRational>> dfac;
072
073
074 ComplexRing<BigRational> cfac;
075
076
077 BigRational eps;
078
079
080 Complex<BigRational> ceps;
081
082
083 GenPolynomial<Complex<BigRational>> a;
084
085
086 GenPolynomial<Complex<BigRational>> b;
087
088
089 GenPolynomial<Complex<BigRational>> c;
090
091
092 GenPolynomial<Complex<BigRational>> d;
093
094
095 GenPolynomial<Complex<BigRational>> e;
096
097
098 int rl = 1;
099
100
101 int kl = 3;
102
103
104 int ll = 3;
105
106
107 int el = 5;
108
109
110 float q = 0.7f;
111
112
113 @Override
114 protected void setUp() {
115 a = b = c = d = e = null;
116 cfac = new ComplexRing<BigRational>(new BigRational(1));
117 String[] vars = new String[] { "x" };
118 dfac = new GenPolynomialRing<Complex<BigRational>>(cfac, rl, to, vars);
119 eps = Power.positivePower(new BigRational(1L, 10L), BigDecimal.DEFAULT_PRECISION);
120 ceps = new Complex<BigRational>(cfac,eps);
121 }
122
123
124 @Override
125 protected void tearDown() {
126 a = b = c = d = e = null;
127 dfac = null;
128 cfac = null;
129 eps = null;
130 }
131
132
133 /**
134 * Test root bound.
135 *
136 */
137 public void testRootBound() {
138 //a = dfac.random(kl, ll, el, q);
139 a = dfac.univariate(0, 2L).sum(dfac.getONE()); // x^2 + 1
140 //System.out.println("a = " + a);
141
142 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
143
144 Complex<BigRational> Mb = cr.rootBound(a);
145 BigRational M = Mb.getRe();
146
147 //System.out.println("M = " + M);
148 assertTrue("M >= 1 ", M.compareTo(BigRational.ONE) >= 0);
149
150 //a = a.monic();
151 a = a.multiply(dfac.fromInteger(5));
152 //System.out.println("a = " + a);
153 M = cr.rootBound(a).getRe();
154
155 //System.out.println("M = " + M);
156 assertTrue("M >= 1 ", M.compareTo(BigRational.ONE) >= 0);
157 }
158
159
160 /**
161 * Test Cauchy index.
162 *
163 */
164 public void testCauchyIndex() {
165 a = dfac.random(kl, ll, el, q);
166 b = dfac.random(kl, ll, el, q);
167 //a = dfac.univariate(0,2L).sum(dfac.getONE()); // x^2 + 1
168 //System.out.println("a = " + a);
169 //System.out.println("b = " + b);
170
171 BigRational l = new BigRational(0);
172 BigRational r = new BigRational(1);
173 GenPolynomialRing<BigRational> fac = new GenPolynomialRing<BigRational>(l, dfac);
174
175 ComplexRootsSturm<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
176
177 GenPolynomial<BigRational> f = PolyUtil.<BigRational> realPartFromComplex(fac, a);
178 GenPolynomial<BigRational> g = PolyUtil.<BigRational> imaginaryPartFromComplex(fac, b);
179 //System.out.println("re(a) = " + f);
180 //System.out.println("im(b) = " + g);
181
182 long ci = cr.indexOfCauchy(l, r, g, f);
183 //System.out.println("ci = " + ci);
184
185 assertTrue("ci >= 0 ", ci >= -a.degree(0));
186 }
187
188
189 /**
190 * Test Routh.
191 *
192 */
193 public void testRouth() {
194 ComplexRootsSturm<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
195
196 //a = dfac.random(kl, ll, el, q);
197 //b = dfac.random(kl, ll, el, q);
198 Complex<BigRational> I = cfac.getIMAG();
199 GenPolynomial<Complex<BigRational>> X = dfac.univariate(0);
200 //System.out.println("I = " + I);
201 //System.out.println("X = " + X);
202
203 //a = dfac.univariate(0,2L).sum( dfac.getONE().multiply(I) ); // x^2 + i
204
205 //b = X.subtract( dfac.getONE().multiply( I ) ); // x - i
206 b = X.subtract(dfac.getONE().multiply(I.negate())); // x + i
207 c = X.subtract(dfac.getONE().multiply(I.multiply(cfac.fromInteger(3)))); // x - 3i
208 d = X.subtract(dfac.getONE().multiply(I.multiply(cfac.fromInteger(4)))); // x - 4i
209 e = X.subtract(dfac.getONE().multiply(I.multiply(cfac.fromInteger(5)))); // x - 5i
210
211 a = b.multiply(c).multiply(d).multiply(e);
212 //System.out.println("a = " + a.toScript());
213 //System.out.println("i = " + cfac.getIMAG());
214
215 Complex<BigRational> Mb = cr.rootBound(a);
216 BigRational M = Mb.getRe();
217 //System.out.println("M = " + M);
218
219 BigRational minf = M.negate(); // - infinity
220 BigRational pinf = M; // + infinity
221 GenPolynomialRing<BigRational> fac = new GenPolynomialRing<BigRational>(pinf, dfac);
222
223 GenPolynomial<BigRational> f = PolyUtil.<BigRational> realPartFromComplex(fac, a);
224 GenPolynomial<BigRational> g = PolyUtil.<BigRational> imaginaryPartFromComplex(fac, a);
225 //System.out.println("re(a) = " + f.toScript());
226 //System.out.println("im(a) = " + g.toScript());
227
228 long[] ri = cr.indexOfRouth(minf, pinf, f, g);
229 //System.out.println("ri = [" + ri[0] + ", " + ri[1] + " ]");
230 long deg = ri[0] + ri[1];
231 assertTrue("sum(ri) == deg(a) ", deg >= a.degree(0));
232 }
233
234
235 /**
236 * Test winding number.
237 *
238 */
239 @SuppressWarnings("unchecked")
240 public void testWindingNumber() {
241 ComplexRootsSturm<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
242 Complex<BigRational> I = cfac.getIMAG();
243
244 a = dfac.univariate(0, 2L).sum(cfac.fromInteger(1)); // x^2 + 1
245 //a = dfac.random(kl, ll, el, q);
246 //a = dfac.univariate(0,2L).subtract(cfac.getONE()); // x^2 - 1
247 //a = dfac.univariate(0,2L).subtract(I); // x^2 - I
248 //a = dfac.univariate(0,1L); // x
249 //System.out.println("a = " + a);
250
251 Complex<BigRational> Mb = cr.rootBound(a);
252 BigRational M = Mb.getRe();
253 //System.out.println("M = " + M);
254 BigRational eps = new BigRational(1, 1000);
255 //System.out.println("eps = " + eps);
256
257 Complex<BigRational>[] corner = new Complex[4];
258
259 corner[0] = new Complex<BigRational>(cfac, M.negate(), M); // nw
260 corner[1] = new Complex<BigRational>(cfac, M.negate(), M.negate()); // sw
261 corner[2] = new Complex<BigRational>(cfac, M, M.negate()); // se
262 corner[3] = new Complex<BigRational>(cfac, M, M); // ne
263
264 long v = 0;
265 try {
266 v = cr.windingNumber(new Rectangle<BigRational>(corner), a);
267 } catch (InvalidBoundaryException e) {
268 fail("" + e);
269 }
270 //System.out.println("winding number = " + v);
271 assertTrue("wind(rect,a) == 2 ", v == 2);
272
273 //if ( true ) return;
274
275 corner[0] = new Complex<BigRational>(cfac, M.negate(), M); // nw
276 corner[1] = new Complex<BigRational>(cfac, M.negate(), eps); // sw
277 corner[2] = new Complex<BigRational>(cfac, M, eps); // se
278 corner[3] = new Complex<BigRational>(cfac, M, M); // ne
279
280 try {
281 v = cr.windingNumber(new Rectangle<BigRational>(corner), a);
282 } catch (InvalidBoundaryException e) {
283 fail("" + e);
284 }
285 //System.out.println("winding number = " + v);
286 assertTrue("wind(rect,a) == 1 ", v == 1);
287
288 corner[0] = new Complex<BigRational>(cfac, eps.negate(), eps); // nw
289 corner[1] = new Complex<BigRational>(cfac, eps.negate(), eps.negate()); // sw
290 corner[2] = new Complex<BigRational>(cfac, eps, eps.negate()); // se
291 corner[3] = new Complex<BigRational>(cfac, eps, eps); // ne
292
293 try {
294 v = cr.windingNumber(new Rectangle<BigRational>(corner), a);
295 } catch (InvalidBoundaryException e) {
296 fail("" + e);
297 }
298 //System.out.println("winding number = " + v);
299 assertTrue("wind(rect,a) == 0 ", v == 0);
300 }
301
302
303 /**
304 * Test complex roots, sqrt(-1).
305 *
306 */
307 @SuppressWarnings("unchecked")
308 public void testComplexRootsImag() {
309 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
310 Complex<BigRational> I = cfac.getIMAG();
311
312 a = dfac.univariate(0, 2L).sum(cfac.fromInteger(1)); // x^2 + 1
313 //a = dfac.univariate(0,2L).subtract(cfac.getONE()); // x^2 - 1
314 //a = dfac.univariate(0,2L).subtract(I); // x^2 - I
315 //a = dfac.univariate(0,1L); // x
316 //System.out.println("a = " + a);
317
318 Complex<BigRational> Mb = cr.rootBound(a);
319 BigRational M = Mb.getRe();
320 //System.out.println("M = " + M);
321
322 Complex<BigRational>[] corner = new Complex[4];
323
324 corner[0] = new Complex<BigRational>(cfac, M.negate(), M); // nw
325 corner[1] = new Complex<BigRational>(cfac, M.negate(), M.negate()); // sw
326 corner[2] = new Complex<BigRational>(cfac, M, M.negate()); // se
327 corner[3] = new Complex<BigRational>(cfac, M, M); // ne
328
329 Rectangle<BigRational> rect = new Rectangle<BigRational>(corner);
330
331 List<Rectangle<BigRational>> roots = null;
332 try {
333 roots = cr.complexRoots(rect, a);
334 } catch (InvalidBoundaryException e) {
335 fail("" + e);
336 }
337 //System.out.println("roots = " + roots);
338 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0));
339 }
340
341
342 /**
343 * Test complex roots.
344 */
345 @SuppressWarnings("unchecked")
346 public void testComplexRootsRand() {
347 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
348 Complex<BigRational> I = cfac.getIMAG();
349
350 a = dfac.random(kl, ll, el, q);
351 Squarefree<Complex<BigRational>> engine = SquarefreeFactory
352 .<Complex<BigRational>> getImplementation(cfac);
353 a = engine.squarefreePart(a);
354
355 //a = dfac.univariate(0,2L).subtract(cfac.getONE()); // x^2 - 1
356 //a = dfac.univariate(0,2L).sum(cfac.fromInteger(1)); // x^2 + 1
357 //a = dfac.univariate(0,2L).subtract(I); // x^2 - I
358 //a = dfac.univariate(0,1L); // x
359 //System.out.println("a = " + a);
360
361 Complex<BigRational> Mb = cr.rootBound(a);
362 BigRational M = Mb.getRe();
363 //System.out.println("M = " + M);
364
365 Complex<BigRational>[] corner = new Complex[4];
366
367 corner[0] = new Complex<BigRational>(cfac, M.negate(), M); // nw
368 corner[1] = new Complex<BigRational>(cfac, M.negate(), M.negate()); // sw
369 corner[2] = new Complex<BigRational>(cfac, M, M.negate()); // se
370 corner[3] = new Complex<BigRational>(cfac, M, M); // ne
371
372 Rectangle<BigRational> rect = new Rectangle<BigRational>(corner);
373
374
375 List<Rectangle<BigRational>> roots = null;
376 try {
377 roots = cr.complexRoots(rect, a);
378 } catch (InvalidBoundaryException e) {
379 fail("" + e);
380 }
381 //System.out.println("a = " + a);
382 //System.out.println("roots = " + roots);
383 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0));
384 }
385
386
387 /**
388 * Test complex roots.
389 */
390 public void testComplexRoots() {
391 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
392
393 a = dfac.random(kl, ll, el + 1, q);
394 //System.out.println("a = " + a);
395
396 List<Rectangle<BigRational>> roots = cr.complexRoots(a);
397 //System.out.println("a = " + a);
398 //System.out.println("roots = " + roots);
399 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0));
400 }
401
402
403 /**
404 * Test complex root refinement.
405 */
406 public void testComplexRootRefinement() {
407 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
408
409 a = dfac.random(kl, ll, el - 1, q);
410 //a = dfac.parse("( (x-1)^3 )");
411 Squarefree<Complex<BigRational>> engine = SquarefreeFactory
412 .<Complex<BigRational>> getImplementation(cfac);
413 //System.out.println("a = " + a);
414 a = engine.squarefreePart(a);
415 //System.out.println("a = " + a);
416
417 List<Rectangle<BigRational>> roots = cr.complexRoots(a);
418 //System.out.println("a = " + a);
419 //System.out.println("roots = " + roots);
420 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0));
421
422 BigRational len = new BigRational(1, 1000);
423 //System.out.println("len = " + len);
424
425 for (Rectangle<BigRational> root : roots) {
426 try {
427 Rectangle<BigRational> refine = cr.complexRootRefinement(root, a, len);
428 //System.out.println("refine = " + refine);
429 } catch (InvalidBoundaryException e) {
430 fail("" + e);
431 }
432 }
433 }
434
435
436 /**
437 * Test complex root refinement full.
438 */
439 public void testComplexRootRefinementFull() {
440 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
441
442 a = dfac.random(kl, ll, el - 1, q);
443 //a = dfac.parse("( (x-1)^3 )");
444 //a = dfac.parse("( x^4-2 )");
445 //System.out.println("a = " + a);
446
447 BigRational len = new BigRational(1, 1000);
448 //System.out.println("len = " + len);
449
450 List<Rectangle<BigRational>> refine = cr.complexRoots(a, len);
451 //System.out.println("refine = " + refine);
452 assertTrue("#roots == deg(a) ", refine.size() == a.degree(0));
453 }
454
455
456 /**
457 * Test winding number with wrong precondition.
458 */
459 @SuppressWarnings("unchecked")
460 public void testWindingNumberWrong() {
461 ComplexRootsSturm<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
462 Complex<BigRational> I = cfac.getIMAG();
463
464 a = dfac.univariate(0, 2L).sum(cfac.fromInteger(1)); // x^2 + 1
465 //a = dfac.random(kl, ll, el, q);
466 //a = dfac.univariate(0,2L).subtract(cfac.getONE()); // x^2 - 1
467 //a = dfac.univariate(0,2L).subtract(I); // x^2 - I
468 //a = dfac.univariate(0,1L); // x
469 //System.out.println("a = " + a);
470
471 Complex<BigRational> Mb = cfac.fromInteger(1); //.divide(cfac.fromInteger(2)); //cr.rootBound(a);
472 BigRational M = Mb.getRe();
473 //System.out.println("M = " + M);
474 BigRational eps = new BigRational(1, 1000);
475 //System.out.println("eps = " + eps);
476 BigRational zero = new BigRational();
477 BigRational one = new BigRational(1);
478
479 Complex<BigRational>[] corner = new Complex[4];
480
481 corner[0] = new Complex<BigRational>(cfac, M.negate(), M); // nw
482 corner[1] = new Complex<BigRational>(cfac, M.negate(), zero); // sw
483 corner[2] = new Complex<BigRational>(cfac, M, zero); // se
484 corner[3] = new Complex<BigRational>(cfac, M, M); // ne
485
486 Rectangle<BigRational> rect = new Rectangle<BigRational>(corner);
487 //System.out.println("rect = " + rect);
488
489 try {
490 long v = cr.windingNumber(rect, a);
491 System.out.println("winding number = " + v);
492 fail("wind(rect,a) must throw an exception");
493 } catch (InvalidBoundaryException e) {
494 }
495 }
496
497
498 /**
499 * Test complex root approximation.
500 */
501 public void testComplexRootApproximation() {
502 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
503
504 //a = dfac.random(kl, ll, el-1, q);
505 //a = dfac.parse("( (x-1)*(x-2)*(x-3)*(x - { 0i1 })*(x-5) )*( x^4-2 )");
506 //a = dfac.parse("( (x-1)*(x-2)*(x-3)*( x^4-2 ) )");
507 //a = dfac.parse("( (x-2)*( x^4-2 ) )");
508 a = dfac.parse("( ( x^4-2 ) )");
509 b = dfac.parse("( (x-1)*(x-2)*(x-3) )");
510 c = dfac.parse("( x^4-2 )");
511 d = dfac.parse("( (x - { 0i1 })*(x-5) )");
512 //a = c; // b; //.multiply(c); //.multiply(d);
513 //System.out.println("a = " + a);
514 //System.out.println("b = " + b);
515 //System.out.println("c = " + c);
516 //System.out.println("d = " + d);
517 //a = b.multiply(c).multiply(d);
518 //System.out.println("a = " + a);
519 Squarefree<Complex<BigRational>> engine = SquarefreeFactory
520 .<Complex<BigRational>> getImplementation(cfac);
521 a = engine.squarefreePart(a);
522 //System.out.println("a = " + a);
523
524 eps = eps.multiply(new BigRational(1000000));
525 //System.out.println("eps = " + eps);
526 BigDecimal eps1 = new BigDecimal(eps);
527 BigDecimal eps2 = eps1.multiply(new BigDecimal("10"));
528 //System.out.println("eps1 = " + eps1);
529 //System.out.println("eps2 = " + eps2);
530
531 List<Rectangle<BigRational>> roots = cr.complexRoots(a);
532 //System.out.println("a = " + a);
533 //System.out.println("roots = " + roots);
534 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0));
535
536 for (Rectangle<BigRational> root : roots) {
537 try {
538 Complex<BigDecimal> cd = cr.approximateRoot(root, a, eps);
539 //System.out.println("cd = " + cd);
540 } catch (NoConvergenceException e) {
541 //fail("" + e);
542 }
543 }
544 }
545
546
547 /**
548 * Test complex root approximation full algorithm.
549 */
550 public void testComplexRootApproximationFull() {
551 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
552
553 //a = dfac.random(kl, ll, el-1, q);
554 //a = dfac.parse("( (x-1)*(x-2)*(x-3)*(x - { 0i1 })*(x-5) )*( x^4-2 )");
555 //a = dfac.parse("( (x-1)*(x-2)*(x-3)*( x^4-2 ) )");
556 a = dfac.parse("( (x-2)*( x^4-2 ) )");
557 //a = dfac.parse("( ( x^4-2 ) )");
558 b = dfac.parse("( (x-1)*(x-2)*(x-3) )");
559 c = dfac.parse("( x^4-2 )");
560 d = dfac.parse("( (x - { 0i1 })*(x-5) )");
561 //a = c; // b; //.multiply(c); //.multiply(d);
562 //System.out.println("a = " + a);
563 //System.out.println("b = " + b);
564 //System.out.println("c = " + c);
565 //System.out.println("d = " + d);
566 //a = b.multiply(c).multiply(d);
567 //System.out.println("a = " + a);
568
569 eps = eps.multiply(new BigRational(1000000));
570 //System.out.println("eps = " + eps);
571 BigDecimal eps1 = new BigDecimal(eps);
572 BigDecimal eps2 = eps1.multiply(new BigDecimal("10"));
573 //System.out.println("eps1 = " + eps1);
574 //System.out.println("eps2 = " + eps2);
575
576 List<Complex<BigDecimal>> roots = cr.approximateRoots(a, eps);
577 //System.out.println("a = " + a);
578 //System.out.println("roots = " + roots);
579 //now always true:
580 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0));
581 }
582
583
584 /**
585 * Test complex root approximation full algorithm with Wilkinson
586 * polynomials. p = (x-i0)*(x-i1)*(x-i2)*(x-i3*...*(x-in)
587 */
588 public void testComplexRootApproximationWilkinsonFull() {
589 final int N = 4;
590 d = dfac.getONE();
591 e = dfac.univariate(0);
592
593 BigDecimal br = new BigDecimal();
594 ComplexRing<BigDecimal> cf = new ComplexRing<BigDecimal>(br);
595 Complex<BigDecimal> I = cf.getIMAG();
596 Complex<BigDecimal> cc = null;
597 Complex<BigRational> Ir = cfac.getIMAG();
598
599 List<Complex<BigDecimal>> Rn = new ArrayList<Complex<BigDecimal>>(N);
600 a = d;
601 for (int i = 0; i < N; i++) {
602 cc = cf.fromInteger(i).multiply(I);
603 Rn.add(cc);
604 c = dfac.fromInteger(i).multiply(Ir);
605 b = e.subtract(c);
606 a = a.multiply(b);
607 }
608 //System.out.println("a = " + a);
609 Collections.reverse(Rn);
610 //System.out.println("Rn = " + Rn);
611
612 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
613
614 eps = eps.multiply(new BigRational(100000));
615 //System.out.println("eps = " + eps);
616 BigDecimal eps1 = new BigDecimal(eps);
617 BigDecimal eps2 = eps1.multiply(new BigDecimal("10"));
618 //System.out.println("eps1 = " + eps1);
619 //System.out.println("eps2 = " + eps2);
620
621 List<Complex<BigDecimal>> roots = cr.approximateRoots(a, eps);
622 //System.out.println("a = " + a);
623 //System.out.println("roots = " + roots);
624 //now always true:
625 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0));
626
627 int i = 0;
628 for (Complex<BigDecimal> dd : roots) {
629 Complex<BigDecimal> di = Rn.get(i++);
630 //System.out.print("di = " + di + ", ");
631 //System.out.println("dd = " + dd);
632 assertTrue("|dd - di| < eps ", dd.subtract(di).norm().getRe().compareTo(eps2) <= 0);
633 }
634 }
635
636
637 /**
638 * Test complex root approximation full algorithm with Wilkinson
639 * polynomials, inverse roots. p = (x-1/i1)*(x-1/i2)*(x-1/i3*...*(x-1/in)
640 */
641 public void testComplexRootApproximationWilkinsonInverseFull() {
642 final int N = 5;
643 d = dfac.getONE();
644 e = dfac.univariate(0);
645
646 BigDecimal br = new BigDecimal();
647 ComplexRing<BigDecimal> cf = new ComplexRing<BigDecimal>(br);
648 Complex<BigDecimal> I = cf.getIMAG();
649 Complex<BigDecimal> cc = null;
650 Complex<BigRational> Ir = cfac.getIMAG();
651
652 List<Complex<BigDecimal>> Rn = new ArrayList<Complex<BigDecimal>>(N);
653 a = d;
654 for (int i = 1; i < N; i++) {
655 cc = cf.fromInteger(i).multiply(I);
656 cc = cc.inverse();
657 Rn.add(cc);
658 c = dfac.fromInteger(i).multiply(Ir);
659 c = d.divide(c);
660 b = e.subtract(c);
661 a = a.multiply(b);
662 }
663 //System.out.println("a = " + a);
664 //Collections.sort(Rn);
665 //System.out.println("Rn = " + Rn);
666
667 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
668
669 eps = eps.multiply(new BigRational(100000));
670 //System.out.println("eps = " + eps);
671 BigDecimal eps1 = new BigDecimal(eps);
672 BigDecimal eps2 = eps1.multiply(new BigDecimal("10"));
673 //System.out.println("eps1 = " + eps1);
674 //System.out.println("eps2 = " + eps2);
675
676 List<Complex<BigDecimal>> roots = cr.approximateRoots(a, eps);
677 //System.out.println("a = " + a);
678 //System.out.println("roots = " + roots);
679 //now always true:
680 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0));
681 //Collections.sort(roots);
682 //System.out.println("roots = " + roots);
683
684 for (Complex<BigDecimal> dd : roots) {
685 //System.out.println("dd = " + dd);
686 boolean t = false;
687 for (Complex<BigDecimal> di : Rn) {
688 //System.out.println("di = " + di);
689 t = dd.subtract(di).norm().getRe().compareTo(eps2) <= 0;
690 if (t) {
691 break;
692 }
693 }
694 if (!t) {
695 //assertTrue("|dd - di| < eps ", dd.subtract(di).norm().getRe().compareTo(eps2) <= 0);
696 fail("|dd - di| < eps ");
697 }
698 }
699 }
700
701
702 /**
703 * Test complex root invariant rectangle.
704 */
705 public void testComplexRootInvariant() {
706 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
707
708 a = dfac.random(kl, ll, el - 1, q);
709 b = dfac.random(kl, ll, 2, q);
710 //a = dfac.parse("( (x-1)^3 )");
711 //a = dfac.parse("( x^4-2 )");
712 if ( a.degree() == 0 ) {
713 return;
714 }
715 Squarefree<Complex<BigRational>> engine = SquarefreeFactory
716 .<Complex<BigRational>> getImplementation(cfac);
717 a = engine.squarefreePart(a);
718 b = engine.squarefreePart(b);
719 //System.out.println("a = " + a);
720 //System.out.println("b = " + b);
721
722 List<Rectangle<BigRational>> roots = cr.complexRoots(a);
723 //System.out.println("roots = " + roots);
724 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0));
725
726 Rectangle<BigRational> rect = roots.get(0);
727 //System.out.println("rect = " + rect);
728
729 try {
730 Rectangle<BigRational> ref = cr.invariantRectangle(rect,a,b);
731 //System.out.println("ref = " + ref);
732 } catch (InvalidBoundaryException e) {
733 e.printStackTrace();
734 }
735 }
736
737
738 /**
739 * Test complex root invariant magnitude rectangle.
740 */
741 public void testComplexRootInvariantMagnitude() {
742 ComplexRootsAbstract<BigRational> cr = new ComplexRootsSturm<BigRational>(cfac);
743
744 a = dfac.random(kl, ll, el - 1, q);
745 b = dfac.random(kl, ll, 3, q);
746 //a = dfac.parse("( x^2 + 1 )");
747 //b = dfac.parse("( x - 0i1 )");
748 if ( a.degree() == 0 ) {
749 return;
750 }
751 Squarefree<Complex<BigRational>> engine = SquarefreeFactory
752 .<Complex<BigRational>> getImplementation(cfac);
753 a = engine.squarefreePart(a);
754 b = engine.squarefreePart(b);
755 //System.out.println("a = " + a);
756 //System.out.println("b = " + b);
757
758 List<Rectangle<BigRational>> roots = cr.complexRoots(a);
759 //System.out.println("roots = " + roots);
760 assertTrue("#roots == deg(a) ", roots.size() == a.degree(0));
761
762 Rectangle<BigRational> rect = roots.get(0);
763 //System.out.println("rect = " + rect);
764
765 try {
766 Rectangle<BigRational> ref = cr.invariantMagnitudeRectangle(rect,a,b,eps);
767 //System.out.println("ref = " + ref);
768 Complex<BigRational> mag = cr.complexRectangleMagnitude(ref,a,b);
769 //System.out.println("mag = " + mag);
770 Complex<BigRational> cmag = cr.complexMagnitude(ref,a,b,eps);
771 //System.out.println("cmag = " + cmag);
772 assertEquals("mag == cmag: " + cmag, mag, cmag);
773 BigRational rmag = cmag.getRe();
774 //System.out.println("rmag = " + new BigDecimal(cmag.getRe()) + " i " + new BigDecimal(cmag.getIm()));
775 } catch (InvalidBoundaryException e) {
776 e.printStackTrace();
777 }
778 }
779
780 }