001/*
002 * $Id$
003 */
004
005package edu.jas.commons.math;
006
007
008import java.math.MathContext;
009import java.util.ArrayList;
010import java.util.List;
011
012import org.apache.commons.math3.linear.Array2DRowRealMatrix;
013import org.apache.commons.math3.linear.EigenDecomposition;
014import org.apache.commons.math3.linear.RealMatrix;
015
016import edu.jas.arith.BigDecimal;
017import edu.jas.arith.Rational;
018import edu.jas.poly.Complex;
019import edu.jas.poly.ComplexRing;
020import edu.jas.poly.GenPolynomial;
021import edu.jas.ps.UnivPowerSeries;
022import edu.jas.ps.UnivPowerSeriesRing;
023import edu.jas.structure.RingElem;
024
025
026/**
027 * Algorithms related to polynomial roots. Conversion to commons-math classes
028 * and delegation to commons-math algorithms.
029 * @param <C> ring element type
030 * @author Heinz Kredel
031 */
032
033public class Roots<C extends RingElem<C> & Rational> {
034
035
036    /**
037     * Complex root approximation using companion matrix.
038     * @param a squarefree univariate polynomial
039     * @return list of approximations of complex roots of the polynomial
040     */
041    public List<Complex<BigDecimal>> complexRoots(GenPolynomial<C> a) {
042        List<Complex<BigDecimal>> r = new ArrayList<Complex<BigDecimal>>();
043        ComplexRing<BigDecimal> cr = new ComplexRing<BigDecimal>(new BigDecimal(0.0, MathContext.DECIMAL64));
044        Complex<BigDecimal> cc = new Complex<BigDecimal>(cr);
045
046        if (a == null || a.isZERO()) {
047            r.add(cc);
048            return r;
049        }
050        if (a.isConstant()) {
051            return r;
052        }
053        a = a.monic();
054        UnivPowerSeriesRing<C> pr = new UnivPowerSeriesRing<C>(a.ring);
055        UnivPowerSeries<C> ps = pr.fromPolynomial(a);
056
057        // Construct the companion matrix
058        int N = (int) a.degree();
059        RealMatrix A = new Array2DRowRealMatrix(N, N);
060        for (int i = 0; i < N; i++) {
061            A.setEntry(i, N - 1, -ps.coefficient(i).getRational().doubleValue());
062        }
063        for (int i = 1; i < N; i++) {
064            A.setEntry(i, i - 1, 1.0);
065        }
066        //System.out.println("A = " + A);
067
068        // compute eigenvalues
069        EigenDecomposition ed = new EigenDecomposition(A);
070        double[] realValues = ed.getRealEigenvalues();
071        double[] imagValues = ed.getImagEigenvalues();
072
073        //RealMatrix V = ed.getV();
074        //System.out.println("V = " + V);
075        //RealMatrix D = ed.getD();
076        //System.out.println("D = " + D); 
077
078        // construct root list
079        for (int i = 0; i < N; i++) {
080            cc = new Complex<BigDecimal>(cr, new BigDecimal(realValues[i], MathContext.DECIMAL64),
081                            new BigDecimal(imagValues[i], MathContext.DECIMAL64));
082            //System.out.println("cc = " + cc + ", re = " + realValues[i] + ", im = " + imagValues[i]);
083            r.add(cc);
084        }
085        return r;
086    }
087
088}