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}