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.arith.Roots;
015import edu.jas.poly.GenPolynomial;
016import edu.jas.poly.GenPolynomialRing;
017import edu.jas.poly.TermOrder;
018import edu.jas.structure.Power;
019
020import junit.framework.Test;
021import junit.framework.TestCase;
022import junit.framework.TestSuite;
023
024
025/**
026 * RealRoot tests with JUnit.
027 * @author Heinz Kredel
028 */
029
030public class RealRootTest extends TestCase {
031
032
033    /**
034     * main.
035     */
036    public static void main(String[] args) {
037        junit.textui.TestRunner.run(suite());
038    }
039
040
041    /**
042     * Constructs a <CODE>RealRootTest</CODE> object.
043     * @param name String.
044     */
045    public RealRootTest(String name) {
046        super(name);
047    }
048
049
050    /**
051     */
052    public static Test suite() {
053        TestSuite suite = new TestSuite(RealRootTest.class);
054        return suite;
055    }
056
057
058    TermOrder to = new TermOrder(TermOrder.INVLEX);
059
060
061    GenPolynomialRing<BigRational> dfac;
062
063
064    BigRational ai, bi, ci, di, ei, eps;
065
066
067    GenPolynomial<BigRational> a, b, c, d, e;
068
069
070    int rl = 1;
071
072
073    int kl = 5;
074
075
076    int ll = 7;
077
078
079    int el = 7;
080
081
082    float q = 0.7f;
083
084
085    @Override
086    protected void setUp() {
087        a = b = c = d = e = null;
088        ai = bi = ci = di = ei = null;
089        String[] vars = new String[] { "x" };
090        dfac = new GenPolynomialRing<BigRational>(new BigRational(1), rl, to, vars);
091        // eps = new BigRational(1L,1000000L*1000000L*1000000L);
092        eps = Power.positivePower(new BigRational(1L, 10L), BigDecimal.DEFAULT_PRECISION);
093    }
094
095
096    @Override
097    protected void tearDown() {
098        a = b = c = d = e = null;
099        ai = bi = ci = di = ei = null;
100        dfac = null;
101        eps = null;
102    }
103
104
105    /**
106     * Test Sturm sequence.
107     */
108    public void testSturmSequence() {
109        a = dfac.random(kl, ll, el, q);
110        //System.out.println("a = " + a);
111
112        RealRootsSturm<BigRational> rrs = new RealRootsSturm<BigRational>();
113
114        List<GenPolynomial<BigRational>> S = rrs.sturmSequence(a);
115        //System.out.println("S = " + S);
116
117        try {
118            b = a.remainder(S.get(0));
119        } catch (Exception e) {
120            fail("not S(0)|f " + e);
121        }
122        assertTrue("a mod S(0) == 0 ", b.isZERO());
123
124        assertTrue("S(-1) == 1 ", S.get(S.size() - 1).isConstant());
125    }
126
127
128    /**
129     * Test root bound.
130     */
131    public void testRootBound() {
132        a = dfac.random(kl, ll, el, q);
133        //System.out.println("a = " + a);
134
135        RealRootsAbstract<BigRational> rr = new RealRootsSturm<BigRational>();
136
137        // used root bound
138        BigRational M = rr.realRootBound(a);
139        //System.out.println("M = " + M);
140        assertTrue("M >= 1 ", M.compareTo(BigRational.ONE) >= 0);
141        Interval<BigRational> v1 = new Interval<BigRational>(M.negate(), M);
142        long r1 = rr.realRootCount(v1, a);
143        //System.out.println("v1 = " + v1 + ", r1 = " + r1);
144
145        a = a.monic();
146        //System.out.println("a = " + a);
147        BigDecimal ar = M.getDecimal();
148        //System.out.println("ar = " + ar);
149        assertTrue("ar >= 1 ", ar.compareTo(BigDecimal.ONE) >= 0);
150
151        // maxNorm root bound
152        BigRational mr = a.maxNorm().getRational().sum(BigRational.ONE);
153        BigDecimal dr = mr.getDecimal();
154        //System.out.println("dr = " + dr);
155        //assertTrue("ar >= maxNorm(a): " + (ar.subtract(dr)), ar.compareTo(dr) >= 0);
156        Interval<BigRational> v2 = new Interval<BigRational>(mr.negate(), mr);
157        long r2 = rr.realRootCount(v2, a);
158        //System.out.println("v2 = " + v2 + ", r2 = " + r2);
159        assertTrue("r1 == r2: " + (r2 - r1), r1 == r2);
160
161        // squareNorm root bound
162        BigRational qr = a.squareNorm().getRational();
163        BigDecimal ir = Roots.sqrt(qr.getDecimal());
164        //qr = Roots.sqrt(qr);
165        //System.out.println("ir = " + ir);
166        //assertTrue("ar >= squareNorm(a): " + (ar.subtract(ir)), ar.compareTo(ir) >= 0);
167        Interval<BigRational> v3 = new Interval<BigRational>(qr.negate(), qr);
168        long r3 = rr.realRootCount(v3, a);
169        //System.out.println("v3 = " + v3 + ", r3 = " + r3);
170        assertTrue("r1 == r3: " + (r3 - r1), r1 == r3);
171
172        // sumNorm root bound
173        BigRational pr = a.sumNorm().getRational();
174        BigDecimal sr = pr.getDecimal();
175        //System.out.println("sr = " + sr);
176        //assertTrue("ar >= squareNorm(a): " + (ar.subtract(sr)), ar.compareTo(sr) >= 0);
177        Interval<BigRational> v4 = new Interval<BigRational>(pr.negate(), pr);
178        long r4 = rr.realRootCount(v4, a);
179        //System.out.println("v4 = " + v4 + ", r4 = " + r4);
180        assertTrue("r1 == r4: " + (r4 - r1), r1 == r4);
181
182        // minimal root bound
183        BigDecimal dri = dr.sum(BigDecimal.ONE).inverse();
184        //System.out.println("dri = " + dri + ", sign(dri) = " + dri.signum());
185        assertTrue("minimal root > 0: " + dri, dri.signum() > 0);
186        BigDecimal mri = rr.realMinimalRootBound(a).getDecimal();
187        //System.out.println("mri = " + mri + ", sign(mri) = " + mri.signum());
188        BigDecimal s = dri.subtract(mri).abs();
189        eps = eps.multiply(BigRational.ONE.fromInteger(10));
190        //System.out.println("s = " + s + ", eps = " + eps.getDecimal());
191        assertTrue("minimal root: " + dri, s.compareTo(eps.getDecimal()) < 0);
192
193        // minimal root separation
194        long n = a.degree();
195        if (n > 0) {
196            BigDecimal sep = sr.sum(BigDecimal.ONE).power(2 * n).multiply(sr.fromInteger(n).power(n + 1))
197                            .inverse();
198            //System.out.println("sep = " + sep + ", sign(sep) = " + sep.signum());
199            assertTrue("separation(a) > 0: " + sep, sep.signum() > 0);
200            BigDecimal sri = rr.realMinimalRootSeparation(a).getDecimal();
201            BigDecimal ss = sep.subtract(sri).abs();
202            assertTrue("minimal separation: " + sep, ss.compareTo(eps.getDecimal()) < 0);
203        }
204    }
205
206
207    /**
208     * Test real root isolation.
209     */
210    public void testRealRootIsolation() {
211        a = dfac.random(kl, ll * 2, el * 2, q);
212        //a = a.multiply( dfac.univariate(0) );
213        //System.out.println("a = " + a);
214
215        RealRoots<BigRational> rr = new RealRootsSturm<BigRational>();
216
217        List<Interval<BigRational>> R = rr.realRoots(a);
218        //System.out.println("R = " + R);
219        //assertTrue("#roots >= 0 ", R.size() >= 0);
220        assertTrue("#roots >= 0 ", R != null);
221    }
222
223
224    /**
225     * Test Thom lemma real root sign sequence.
226     */
227    public void testRealRootSignSequence() {
228        a = dfac.random(kl, ll * 2, el * 2, q);
229        if (a.degree() % 2 == 0) {
230            a = a.multiply(dfac.univariate(0).subtract(dfac.getONE()));
231        }
232        if (a.trailingExpVector().degree() > 0) {
233            a = a.subtract(dfac.getONE()); // exclude root 0
234        }
235        //System.out.println("a = " + a);
236        RealRootsAbstract<BigRational> rr = new RealRootsSturm<BigRational>();
237
238        List<Interval<BigRational>> R = rr.realRoots(a);
239        //System.out.println("R = " + R);
240        //assertTrue("#roots >= 0 ", R.size() >= 0);
241        assertTrue("#roots >= 0 ", R != null);
242
243        int l = R.size();
244        Interval<BigRational> v = R.get(l - 1);
245        Interval<BigRational> u = R.get(0);
246        if (u.left.isZERO() && u.right.isZERO()) {
247            Interval<BigRational> w = v;
248            v = u;
249            u = w;
250        }
251        Interval<BigRational> vm = new Interval<BigRational>(u.left, v.right);
252        //System.out.println("v  = " + v);
253        //System.out.println("u  = " + u);
254        //System.out.println("vm = " + vm);
255        long rc = rr.realRootCount(vm, a);
256        //System.out.println("rc = " + rc);
257        assertTrue("root number: R = " + R + ", a = " + a + ", rc = " + rc, rc == l);
258        long rn = rr.realRootNumber(a, vm);
259        assertTrue("root number == " + rn, rn == l);
260
261        long d = a.degree();
262        List<GenPolynomial<BigRational>> fs = rr.fourierSequence(a);
263        //System.out.println("fs = " + fs);
264        assertTrue("len(fs) == " + (d + 1 - fs.size()), fs.size() == (d + 1));
265
266        //List<Integer> ss = rr.signSequence(a, v);
267        //System.out.println("ss = " + ss);
268        //assertTrue("len(ss) == " + (d-ss.size()), ss.size() == d);
269        for (Interval<BigRational> t : R) {
270            List<Integer> ss = rr.signSequence(a, t);
271            //System.out.println("ss = " + ss);
272            assertTrue("len(ss) == " + (d - ss.size()), ss.size() == d);
273        }
274    }
275
276
277    /**
278     * Test real root isolation Wilkinson polynomials. p =
279     * (x-0)*(x-1)*(x-2)*(x-3)*...*(x-n)
280     */
281    public void testRealRootIsolationWilkinson() {
282        final int N = 10;
283        d = dfac.getONE();
284        e = dfac.univariate(0);
285
286        List<Interval<BigRational>> Rn = new ArrayList<Interval<BigRational>>(N);
287        a = d;
288        for (int i = 0; i < N; i++) {
289            c = dfac.fromInteger(i);
290            Rn.add(new Interval<BigRational>(c.leadingBaseCoefficient()));
291            b = e.subtract(c);
292            a = a.multiply(b);
293        }
294        //System.out.println("a = " + a);
295
296        RealRoots<BigRational> rr = new RealRootsSturm<BigRational>();
297
298        List<Interval<BigRational>> R = rr.realRoots(a);
299        //System.out.println("R = " + R);
300
301        assertTrue("#roots = " + N + " ", R.size() == N);
302
303        eps = eps.multiply(new BigRational("1/10"));
304        //System.out.println("eps = " + eps);
305
306        R = rr.refineIntervals(R, a, eps);
307        //System.out.println("R = " + R);
308        int i = 0;
309        for (Interval<BigRational> v : R) {
310            BigDecimal dd = v.toDecimal();
311            BigDecimal di = Rn.get(i++).toDecimal();
312            //System.out.println("v  = " + dd);
313            //System.out.println("vi = " + di);
314            //System.out.println("|dd - di| < eps: " + dd.compareTo(di));
315            assertTrue("|dd - di| < eps ", dd.compareTo(di) == 0);
316        }
317    }
318
319
320    /**
321     * Test real root isolation Wilkinson polynomials inverse. p =
322     * (x-1)*(x-1/2)*(x-1/3)*...*(x-1/n)
323     */
324    public void testRealRootIsolationWilkinsonInverse() {
325        final int N = 9;
326        d = dfac.getONE();
327        e = dfac.univariate(0);
328
329        List<Interval<BigRational>> Rn = new ArrayList<Interval<BigRational>>(N);
330        a = d;
331        for (int i = 1; i < N; i++) { // use only for i > 0, since reverse
332            c = dfac.fromInteger(i);
333            if (i != 0) {
334                c = d.divide(c);
335            }
336            Rn.add(new Interval<BigRational>(c.leadingBaseCoefficient()));
337            b = e.subtract(c);
338            a = a.multiply(b);
339        }
340        //System.out.println("a = " + a);
341        //System.out.println("Rn = " + Rn);
342        Collections.reverse(Rn);
343        //System.out.println("Rn = " + Rn);
344
345        RealRoots<BigRational> rr = new RealRootsSturm<BigRational>();
346
347        List<Interval<BigRational>> R = rr.realRoots(a);
348        //System.out.println("R = " + R);
349
350        assertTrue("#roots = " + (N - 1) + " ", R.size() == (N - 1));
351
352        eps = eps.multiply(new BigRational("1/100"));
353        //System.out.println("eps = " + eps);
354
355        R = rr.refineIntervals(R, a, eps);
356        //System.out.println("R = " + R);
357        int i = 0;
358        for (Interval<BigRational> v : R) {
359            BigDecimal dd = v.toDecimal(); //.sum(eps1);
360            BigDecimal di = Rn.get(i++).toDecimal();
361            //System.out.println("v  = " + dd);
362            //System.out.println("vi = " + di);
363            //System.out.println("|dd - di| < eps: " + dd.compareTo(di));
364            assertTrue("|dd - di| < eps ", dd.compareTo(di) == 0);
365        }
366    }
367
368
369    /**
370     * Test real algebraic number sign.
371     */
372    public void testRealAlgebraicNumberSign() {
373        d = dfac.fromInteger(2);
374        e = dfac.univariate(0);
375
376        a = e.multiply(e);
377        // a = a.multiply(e).multiply(e).multiply(e);
378        a = a.subtract(d); // x^2 -2
379        //System.out.println("a = " + a);
380
381        RealRoots<BigRational> rr = new RealRootsSturm<BigRational>();
382
383        ai = new BigRational(1);
384        bi = new BigRational(2);
385        Interval<BigRational> iv = new Interval<BigRational>(ai, bi);
386        //System.out.println("iv = " + iv);
387        assertTrue("sign change", rr.signChange(iv, a));
388
389        b = dfac.random(kl, (int) a.degree() + 1, (int) a.degree(), 1.0f);
390        //b = dfac.getZERO();
391        //b = dfac.random(kl,ll,el,q);
392        //b = b.multiply(b);
393        //b = b.abs().negate();
394        //System.out.println("b = " + b);
395        if (b.isZERO()) {
396            int s = rr.realSign(iv, a, b);
397            assertTrue("algebraic sign", s == 0);
398            return;
399        }
400
401        int as = rr.realSign(iv, a, b);
402        //System.out.println("as = " + as);
403        // how to test?
404        int asn = rr.realSign(iv, a, b.negate());
405        //System.out.println("asn = " + asn);
406        assertTrue("algebraic sign", as != asn);
407
408        iv = new Interval<BigRational>(bi.negate(), ai.negate());
409        //System.out.println("iv = " + iv);
410        assertTrue("sign change", rr.signChange(iv, a));
411
412        int as1 = rr.realSign(iv, a, b);
413        //System.out.println("as1 = " + as1);
414        // how to test?
415        int asn1 = rr.realSign(iv, a, b.negate());
416        //System.out.println("asn1 = " + asn1);
417        assertTrue("algebraic sign", as1 != asn1);
418
419        assertTrue("algebraic sign", as * as1 == asn * asn1);
420    }
421
422
423    /**
424     * Test real root isolation and decimal refinement of Wilkinson polynomials.
425     * p = (x-0)*(x-1)*(x-2)*(x-3)*...*(x-n)
426     */
427    public void testRealRootIsolationDecimalWilkinson() {
428        final int N = 10;
429        d = dfac.getONE();
430        e = dfac.univariate(0);
431
432        List<Interval<BigRational>> Rn = new ArrayList<Interval<BigRational>>(N);
433        a = d;
434        for (int i = 0; i < N; i++) {
435            c = dfac.fromInteger(i);
436            Rn.add(new Interval<BigRational>(c.leadingBaseCoefficient()));
437            b = e.subtract(c);
438            a = a.multiply(b);
439        }
440        //System.out.println("a = " + a);
441
442        RealRootsAbstract<BigRational> rr = new RealRootsSturm<BigRational>();
443
444        List<Interval<BigRational>> R = rr.realRoots(a);
445        //System.out.println("R = " + R);
446
447        assertTrue("#roots = " + N + " ", R.size() == N);
448
449        eps = eps.multiply(new BigRational(100000));
450        //System.out.println("eps = " + eps);
451        BigDecimal eps1 = new BigDecimal(eps);
452        BigDecimal eps2 = eps1.multiply(new BigDecimal("100"));
453        //System.out.println("eps1 = " + eps1);
454        //System.out.println("eps2 = " + eps2);
455
456        try {
457            int i = 0;
458            for (Interval<BigRational> v : R) {
459                //System.out.println("v = " + v);
460                BigDecimal dd = rr.approximateRoot(v, a, eps);
461                BigDecimal di = Rn.get(i++).toDecimal();
462                //System.out.println("di = " + di);
463                //System.out.println("dd = " + dd);
464                assertTrue("|dd - di| < eps ", dd.subtract(di).abs().compareTo(eps2) <= 0);
465            }
466        } catch (NoConvergenceException e) {
467            fail(e.toString());
468        }
469    }
470
471
472    /**
473     * Test real root isolation and decimal refinement of Wilkinson polynomials,
474     * inverse roots. p = (x-1)*(x-1/2)*(x-1/3)*...*(x-1/n)
475     */
476    public void testRealRootIsolationDecimalWilkinsonInverse() {
477        final int N = 10;
478        d = dfac.getONE();
479        e = dfac.univariate(0);
480
481        List<Interval<BigRational>> Rn = new ArrayList<Interval<BigRational>>(N);
482        a = d;
483        for (int i = 1; i < N; i++) { // use only for i > 0, since reverse
484            c = dfac.fromInteger(i);
485            if (i != 0) {
486                c = d.divide(c);
487            }
488            Rn.add(new Interval<BigRational>(c.leadingBaseCoefficient()));
489            b = e.subtract(c);
490            a = a.multiply(b);
491        }
492        //System.out.println("a = " + a);
493        //System.out.println("Rn = " + Rn);
494        Collections.reverse(Rn);
495        //System.out.println("Rn = " + Rn);
496
497        RealRootsAbstract<BigRational> rr = new RealRootsSturm<BigRational>();
498
499        List<Interval<BigRational>> R = rr.realRoots(a);
500        //System.out.println("R = " + R);
501
502        assertTrue("#roots = " + (N - 1) + " ", R.size() == (N - 1));
503
504        eps = eps.multiply(new BigRational(1000000));
505        //System.out.println("eps = " + eps);
506        BigDecimal eps1 = new BigDecimal(eps);
507        BigDecimal eps2 = eps1.multiply(new BigDecimal("10"));
508        //System.out.println("eps1 = " + eps1);
509        //System.out.println("eps2 = " + eps2);
510
511        try {
512            int i = 0;
513            for (Interval<BigRational> v : R) {
514                //System.out.println("v = " + v);
515                BigDecimal dd = rr.approximateRoot(v, a, eps);
516                BigDecimal di = Rn.get(i++).toDecimal();
517                //System.out.println("di = " + di);
518                //System.out.println("dd = " + dd);
519                assertTrue("|dd - di| < eps ", dd.subtract(di).abs().compareTo(eps2) <= 0);
520            }
521        } catch (NoConvergenceException e) {
522            fail(e.toString());
523        }
524    }
525
526
527    /**
528     * Test real root isolation and decimal refinement of Wilkinson polynomials,
529     * all roots. p = (x-0)*(x-1)*(x-2)*(x-3)*...*(x-n)
530     */
531    public void testRealRootIsolationDecimalWilkinsonAll() {
532        final int N = 10;
533        d = dfac.getONE();
534        e = dfac.univariate(0);
535
536        List<Interval<BigRational>> Rn = new ArrayList<Interval<BigRational>>(N);
537        a = d;
538        for (int i = 0; i < N; i++) {
539            c = dfac.fromInteger(i);
540            Rn.add(new Interval<BigRational>(c.leadingBaseCoefficient()));
541            b = e.subtract(c);
542            a = a.multiply(b);
543        }
544        //System.out.println("a = " + a);
545
546        RealRootsAbstract<BigRational> rr = new RealRootsSturm<BigRational>();
547
548        eps = eps.multiply(new BigRational(10000));
549        //System.out.println("eps = " + eps);
550        BigDecimal eps1 = new BigDecimal(eps);
551        BigDecimal eps2 = eps1.multiply(new BigDecimal("100"));
552        //System.out.println("eps1 = " + eps1);
553        //System.out.println("eps2 = " + eps2);
554
555        List<BigDecimal> R = null;
556        R = rr.approximateRoots(a, eps);
557        //System.out.println("R = " + R);
558        assertTrue("#roots = " + N + " ", R.size() == N);
559
560        int i = 0;
561        for (BigDecimal dd : R) {
562            //System.out.println("dd = " + dd);
563            BigDecimal di = Rn.get(i++).toDecimal();
564            //System.out.println("di = " + di);
565            assertTrue("|dd - di| < eps ", dd.subtract(di).abs().compareTo(eps2) <= 0);
566        }
567        boolean t = rr.isApproximateRoot(R, a, eps);
568        assertTrue("some |a(dd)| < eps ", t);
569    }
570
571
572    /**
573     * Test real root isolation and decimal refinement of Wilkinson polynomials,
574     * inverse roots, all roots. p = (x-1)*(x-1/2)*(x-1/3)*...*(x-1/n)
575     */
576    public void testRealRootIsolationDecimalWilkinsonInverseAll() {
577        final int N = 10;
578        d = dfac.getONE();
579        e = dfac.univariate(0);
580
581        List<Interval<BigRational>> Rn = new ArrayList<Interval<BigRational>>(N);
582        a = d;
583        for (int i = 1; i < N; i++) { // use only for i > 0, since reverse
584            c = dfac.fromInteger(i);
585            if (i != 0) {
586                c = d.divide(c);
587            }
588            Rn.add(new Interval<BigRational>(c.leadingBaseCoefficient()));
589            b = e.subtract(c);
590            a = a.multiply(b);
591        }
592        //System.out.println("a = " + a);
593        //System.out.println("Rn = " + Rn);
594        Collections.reverse(Rn);
595        //System.out.println("Rn = " + Rn);
596
597        RealRootsAbstract<BigRational> rr = new RealRootsSturm<BigRational>();
598
599        eps = eps.multiply(new BigRational(1000000));
600        //System.out.println("eps = " + eps);
601        BigDecimal eps1 = new BigDecimal(eps);
602        BigDecimal eps2 = eps1.multiply(new BigDecimal("10"));
603        //System.out.println("eps1 = " + eps1);
604        //System.out.println("eps2 = " + eps2);
605
606        List<BigDecimal> R = null;
607        R = rr.approximateRoots(a, eps);
608        //System.out.println("R = " + R);
609        assertTrue("#roots = " + (N - 1) + " ", R.size() == (N - 1));
610
611        int i = 0;
612        for (BigDecimal dd : R) {
613            //System.out.println("dd = " + dd);
614            BigDecimal di = Rn.get(i++).toDecimal();
615            //System.out.println("di = " + di);
616            assertTrue("|dd - di| < eps ", dd.subtract(di).abs().compareTo(eps2) <= 0);
617        }
618        boolean t = rr.isApproximateRoot(R, a, eps);
619        assertTrue("some |a(dd)| < eps ", t);
620    }
621
622}