Source code

001/*
002 * $Id$
003 */
004
005package edu.jas.root;
006
007
008import java.util.ArrayList;
009import java.util.List;
010
011import edu.jas.arith.BigDecimal;
012import edu.jas.arith.BigInteger;
013import edu.jas.arith.BigRational;
014import edu.jas.poly.GenPolynomial;
015import edu.jas.poly.GenPolynomialRing;
016import edu.jas.poly.PolyUtil;
017import edu.jas.structure.NotInvertibleException;
018import edu.jas.structure.Power;
019
020import junit.framework.Test;
021import junit.framework.TestCase;
022import junit.framework.TestSuite;
023
024
025/**
026 * RealAlgebraicNumber Test using JUnit.
027 * @author Heinz Kredel
028 */
029
030public class RealAlgebraicTest 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>RealAlgebraicTest</CODE> object.
043     * @param name String.
044     */
045    public RealAlgebraicTest(String name) {
046        super(name);
047    }
048
049
050    /**
051     * suite.
052     */
053    public static Test suite() {
054        TestSuite suite = new TestSuite(RealAlgebraicTest.class);
055        return suite;
056    }
057
058
059    RealAlgebraicRing<BigRational> fac;
060
061
062    GenPolynomialRing<BigRational> mfac;
063
064
065    RealAlgebraicNumber<BigRational> a, b, c, d, e;
066
067
068    RealAlgebraicNumber<BigRational> alpha;
069
070
071    int rl = 1;
072
073
074    int kl = 10;
075
076
077    int ll = 10;
078
079
080    int el = ll;
081
082
083    float q = 0.5f;
084
085
086    @Override
087    protected void setUp() {
088        a = b = c = d = e = null;
089        BigRational l = new BigRational(1);
090        BigRational r = new BigRational(2);
091        Interval<BigRational> positive = new Interval<BigRational>(l, r);
092        String[] vars = new String[] { "alpha" };
093        mfac = new GenPolynomialRing<BigRational>(new BigRational(1), rl, vars);
094        GenPolynomial<BigRational> mo = mfac.univariate(0, 2);
095        mo = mo.subtract(mfac.fromInteger(2)); // alpha^2 -2 
096        fac = new RealAlgebraicRing<BigRational>(mo, positive);
097        alpha = fac.getGenerator();
098    }
099
100
101    @Override
102    protected void tearDown() {
103        a = b = c = d = e = null;
104        fac = null;
105        alpha = null;
106    }
107
108
109    /**
110     * Test constructor and toString.
111     * 
112     */
113    public void testConstruction() {
114        c = fac.getONE();
115        //System.out.println("c = " + c);
116        //System.out.println("c.getVal() = " + c.getVal());
117        assertTrue("length( c ) = 1", c.number.getVal().length() == 1);
118        assertTrue("isZERO( c )", !c.isZERO());
119        assertTrue("isONE( c )", c.isONE());
120
121        d = fac.getZERO();
122        //System.out.println("d = " + d);
123        //System.out.println("d.getVal() = " + d.getVal());
124        assertTrue("length( d ) = 0", d.number.getVal().length() == 0);
125        assertTrue("isZERO( d )", d.isZERO());
126        assertTrue("isONE( d )", !d.isONE());
127    }
128
129
130    /**
131     * Test random polynomial.
132     */
133    public void testRandom() {
134        for (int i = 0; i < 7; i++) {
135            a = fac.random(el);
136            //System.out.println("a = " + a);
137            if (a.isZERO() || a.isONE()) {
138                continue;
139            }
140            // fac.random(kl*(i+1), ll+2*i, el+i, q );
141            assertTrue("length( a" + i + " ) <> 0", a.number.getVal().length() >= 0);
142            assertTrue(" not isZERO( a" + i + " )", !a.isZERO());
143            assertTrue(" not isONE( a" + i + " )", !a.isONE());
144        }
145    }
146
147
148    /**
149     * Test addition.
150     */
151    public void testAddition() {
152        a = fac.random(ll);
153        b = fac.random(ll);
154
155        c = a.sum(b);
156        d = c.subtract(b);
157        assertEquals("a+b-b = a", a, d);
158
159        c = a.sum(b);
160        d = b.sum(a);
161        assertEquals("a+b = b+a", c, d);
162
163        c = fac.random(ll);
164        d = c.sum(a.sum(b));
165        e = c.sum(a).sum(b);
166        assertEquals("c+(a+b) = (c+a)+b", d, e);
167
168
169        c = a.sum(fac.getZERO());
170        d = a.subtract(fac.getZERO());
171        assertEquals("a+0 = a-0", c, d);
172
173        c = fac.getZERO().sum(a);
174        d = fac.getZERO().subtract(a.negate());
175        assertEquals("0+a = 0+(-a)", c, d);
176    }
177
178
179    /**
180     * Test multiplication.
181     */
182    public void testMultiplication() {
183        a = fac.random(ll);
184        assertTrue("not isZERO( a )", !a.isZERO());
185
186        b = fac.random(ll);
187        assertTrue("not isZERO( b )", !b.isZERO());
188
189        c = b.multiply(a);
190        d = a.multiply(b);
191        assertTrue("not isZERO( c )", !c.isZERO());
192        assertTrue("not isZERO( d )", !d.isZERO());
193
194        //System.out.println("a = " + a);
195        //System.out.println("b = " + b);
196        e = d.subtract(c);
197        assertTrue("isZERO( a*b-b*a ) " + e, e.isZERO());
198
199        assertTrue("a*b = b*a", c.equals(d));
200        assertEquals("a*b = b*a", c, d);
201
202        c = fac.random(ll);
203        //System.out.println("c = " + c);
204        d = a.multiply(b.multiply(c));
205        e = (a.multiply(b)).multiply(c);
206
207        //System.out.println("d = " + d);
208        //System.out.println("e = " + e);
209
210        //System.out.println("d-e = " + d.subtract(c) );
211
212        assertEquals("a(bc) = (ab)c", d, e);
213        assertTrue("a(bc) = (ab)c", d.equals(e));
214
215        c = a.multiply(fac.getONE());
216        d = fac.getONE().multiply(a);
217        assertEquals("a*1 = 1*a", c, d);
218
219        c = a.inverse();
220        d = c.multiply(a);
221        //System.out.println("a = " + a);
222        //System.out.println("c = " + c);
223        //System.out.println("d = " + d);
224        assertEquals("a*1/a = 1", fac.getONE(), d);
225
226        try {
227            a = fac.getZERO().inverse();
228            fail("0 invertible");
229        } catch (NotInvertibleException expected) {
230            // pass
231        }
232    }
233
234
235    /**
236     * Test distributive law.
237     */
238    public void testDistributive() {
239        a = fac.random(ll);
240        b = fac.random(ll);
241        c = fac.random(ll);
242
243        d = a.multiply(b.sum(c));
244        e = a.multiply(b).sum(a.multiply(c));
245
246        assertEquals("a(b+c) = ab+ac", d, e);
247    }
248
249
250    /**
251     * Test sign of real algebraic numbers.
252     */
253    public void testSignum() {
254        a = fac.random(ll);
255        b = fac.random(ll);
256        c = fac.random(ll);
257
258        int sa = a.signum();
259        int sb = b.signum();
260        int sc = c.signum();
261
262        d = a.multiply(b);
263        e = a.multiply(c);
264
265        int sd = d.signum();
266        int se = e.signum();
267
268        assertEquals("sign(a*b) = sign(a)*sign(b) ", sa * sb, sd);
269        assertEquals("sign(a*c) = sign(a)*sign(c) ", sa * sc, se);
270
271        b = a.negate();
272        sb = b.signum();
273        assertEquals("sign(-a) = -sign(a) ", -sa, sb);
274    }
275
276
277    /**
278     * Test compareTo of real algebraic numbers.
279     */
280    public void testCompare() {
281        a = fac.random(ll).abs();
282        b = a.sum(fac.getONE());
283        c = b.sum(fac.getONE());
284
285        int ab = a.compareTo(b);
286        int bc = b.compareTo(c);
287        int ac = a.compareTo(c);
288
289        assertTrue("a < a+1 ", ab < 0);
290        assertTrue("a+1 < a+2 ", bc < 0);
291        assertTrue("a < a+2 ", ac < 0);
292
293        a = a.negate();
294        b = a.sum(fac.getONE());
295        c = b.sum(fac.getONE());
296
297        ab = a.compareTo(b);
298        bc = b.compareTo(c);
299        ac = a.compareTo(c);
300
301        assertTrue("a < a+1 ", ab < 0);
302        assertTrue("a+1 < a+2 ", bc < 0);
303        assertTrue("a < a+2 ", ac < 0);
304    }
305
306
307    /**
308     * Test arithmetic of magnitude of real algebraic numbers.
309     */
310    public void testMagnitude() {
311        a = fac.random(ll);
312        b = fac.random(ll);
313        c = fac.random(ll);
314
315        //BigDecimal ad = new BigDecimal(a.magnitude());
316        //BigDecimal bd = new BigDecimal(b.magnitude());
317        //BigDecimal cd = new BigDecimal(c.magnitude());
318
319        d = a.multiply(b);
320        e = a.sum(b);
321
322        BigDecimal dd = new BigDecimal(d.magnitude());
323        BigDecimal ed = new BigDecimal(e.magnitude());
324
325        BigDecimal dd1 = new BigDecimal(a.magnitude().multiply(b.magnitude()));
326        BigDecimal ed1 = new BigDecimal(a.magnitude().sum(b.magnitude()));
327
328        //System.out.println("ad  = " + ad);
329        //System.out.println("bd  = " + bd);
330        //System.out.println("cd  = " + cd);
331        //System.out.println("dd  = " + dd);
332        //System.out.println("dd1 = " + dd1);
333        //System.out.println("ed  = " + ed);
334        //System.out.println("ed1 = " + ed1);
335
336        //BigRational eps = Power.positivePower(new BigRational(1L,10L),BigDecimal.DEFAULT_PRECISION);
337        BigRational eps = Power.positivePower(new BigRational(1L, 10L), 8);
338        BigDecimal epsd = new BigDecimal(eps);
339
340        BigDecimal rel = dd.abs().sum(dd1.abs());
341        BigDecimal err = dd.subtract(dd1).abs().divide(rel);
342        //System.out.println("rel = " + rel);
343        //System.out.println("|dd-dd1|/rel = " + err + ", eps = " + epsd);
344        assertTrue("mag(a*b) = mag(a)*mag(b): " + dd + ", " + dd1, err.compareTo(epsd) <= 0);
345        assertTrue("mag(a+b) = mag(a)+mag(b): " + ed + ", " + ed1,
346                        ed.subtract(ed1).abs().compareTo(epsd) <= 0);
347
348
349        d = a.divide(b);
350        e = a.subtract(b);
351
352        dd = new BigDecimal(d.magnitude());
353        ed = new BigDecimal(e.magnitude());
354
355        dd1 = new BigDecimal(a.magnitude().divide(b.magnitude()));
356        ed1 = new BigDecimal(a.magnitude().subtract(b.magnitude()));
357
358        //System.out.println("dd  = " + dd);
359        //System.out.println("dd1 = " + dd1);
360        //System.out.println("ed  = " + ed);
361        //System.out.println("ed1 = " + ed1);
362
363        assertTrue("mag(a/b) = mag(a)/mag(b)", dd.subtract(dd1).abs().compareTo(epsd) <= 0);
364        assertTrue("mag(a-b) = mag(a)-mag(b)", ed.subtract(ed1).abs().compareTo(epsd) <= 0);
365    }
366
367
368    /**
369     * Test real root isolation.
370     */
371    public void testRealRootIsolation() {
372        //System.out.println();
373        GenPolynomialRing<RealAlgebraicNumber<BigRational>> dfac;
374        dfac = new GenPolynomialRing<RealAlgebraicNumber<BigRational>>(fac, 1);
375
376        GenPolynomial<RealAlgebraicNumber<BigRational>> ar;
377        ar = dfac.random(3, 5, 7, q);
378        //System.out.println("ar = " + ar);
379        if (ar.degree() % 2 == 0) { // ensure existence of a root
380            ar = ar.multiply(dfac.univariate(0));
381        }
382        //System.out.println("ar = " + ar);
383
384        RealRoots<RealAlgebraicNumber<BigRational>> rrr = new RealRootsSturm<RealAlgebraicNumber<BigRational>>();
385        List<Interval<RealAlgebraicNumber<BigRational>>> R = rrr.realRoots(ar);
386        //System.out.println("R = " + R);
387        assertTrue("#roots >= 0 ", R.size() > 0);
388
389        BigRational eps = Power.positivePower(new BigRational(1L, 10L), BigDecimal.DEFAULT_PRECISION);
390        BigDecimal epsd = new BigDecimal(Power.positivePower(new BigRational(1L, 10L), 14 - ar.degree())); // less
391        //System.out.println("eps = " + epsd);
392
393        R = rrr.refineIntervals(R, ar, eps);
394        //System.out.println("R = " + R);
395        for (Interval<RealAlgebraicNumber<BigRational>> v : R) {
396            RealAlgebraicNumber<BigRational> m = v.middle();
397            //System.out.println("m = " + m);
398            RealAlgebraicNumber<BigRational> n;
399            n = PolyUtil.<RealAlgebraicNumber<BigRational>> evaluateMain(fac, ar, m);
400            //System.out.println("n = " + n);
401            BigRational nr = n.magnitude();
402            //System.out.println("nr = " + nr);
403            BigDecimal nd = new BigDecimal(nr);
404            //System.out.println("nd = " + nd);
405            assertTrue("|nd| < eps: " + nd + " " + epsd, nd.abs().compareTo(epsd) <= 0);
406            //no: assertTrue("n == 0: " + n, n.isZERO());
407        }
408    }
409
410
411    /**
412     * Test real root isolation for Wilkinson like polynomials.
413     * product_{i=1..n}( x - i * alpha )
414     */
415    public void testRealRootIsolationWilkinson() {
416        //System.out.println();
417        GenPolynomialRing<RealAlgebraicNumber<BigRational>> dfac;
418        dfac = new GenPolynomialRing<RealAlgebraicNumber<BigRational>>(fac, 1);
419
420        GenPolynomial<RealAlgebraicNumber<BigRational>> ar, br, cr, dr, er;
421
422        RealRoots<RealAlgebraicNumber<BigRational>> rrr = new RealRootsSturm<RealAlgebraicNumber<BigRational>>();
423
424        final int N = 3;
425        dr = dfac.getONE();
426        er = dfac.univariate(0);
427
428        List<Interval<RealAlgebraicNumber<BigRational>>> Rn;
429        Rn = new ArrayList<Interval<RealAlgebraicNumber<BigRational>>>(N);
430        ar = dr;
431        for (int i = 0; i < N; i++) {
432            cr = dfac.fromInteger(i).multiply(alpha); // i * alpha
433            Rn.add(new Interval<RealAlgebraicNumber<BigRational>>(cr.leadingBaseCoefficient()));
434            br = er.subtract(cr); // ( x - i * alpha )
435            ar = ar.multiply(br);
436        }
437        //System.out.println("ar = " + ar);
438
439        List<Interval<RealAlgebraicNumber<BigRational>>> R = rrr.realRoots(ar);
440        //System.out.println("R = " + R);
441
442        assertTrue("#roots = " + N + " ", R.size() == N);
443
444        BigRational eps = Power.positivePower(new BigRational(1L, 10L), BigDecimal.DEFAULT_PRECISION);
445        eps = eps.multiply(new BigRational("1/10"));
446        //System.out.println("eps = " + eps);
447
448        R = rrr.refineIntervals(R, ar, eps);
449        //System.out.println("R = " + R);
450        int i = 0;
451        for (Interval<RealAlgebraicNumber<BigRational>> v : R) {
452            BigDecimal dd = v.toDecimal();//.sum(eps1);
453            BigDecimal di = Rn.get(i++).toDecimal();
454            //System.out.println("v  = " + dd);
455            //System.out.println("vi = " + di);
456            assertTrue("|dd - di| < eps ", dd.compareTo(di) == 0);
457        }
458    }
459
460
461    /**
462     * Test continued fraction.
463     */
464    public void testContinuedFraction() {
465        final int M = 100;
466        RealAlgebraicNumber<BigRational> x = fac.random(15);
467        //x = fac.parse("5/12");
468        //x = fac.parse("-1");
469        //x = fac.getGenerator(); // alpha = sqrt(2)
470        //System.out.println("x = " + x + ", mag(x) = " + x.decimalMagnitude());
471
472        List<BigInteger> cf = RealArithUtil.continuedFraction(x, M);
473        //System.out.println("cf(" + x + ") = " + cf);
474
475        BigRational a = RealArithUtil.continuedFractionApprox(cf);
476        RealAlgebraicNumber<BigRational> ax = fac.fromRational(a);
477        //System.out.println("ax = " + ax + ", mag(ax) = " + ax.decimalMagnitude());
478
479        BigRational eps = (new BigRational("1/10")).power(M / 2);
480        RealAlgebraicNumber<BigRational> dx = x.subtract(ax).abs();
481        // check ax > 1:
482        dx = dx.divide(ax.abs());
483        //System.out.println("dx = " + dx + ", mag(dx) = " + dx.decimalMagnitude());
484        BigRational dr = dx.magnitude();
485        assertTrue("|a - approx(cf(a))| = " + dr, dr.compareTo(eps) < 1);
486    }
487
488}