Source code

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