Makepoly.java

/*
 * $Id: Makepoly.java,v 1.29 2008/07/16 15:40:00 koga Exp $
 * 
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.matrix;

import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.ComplexPolynomialMatrix;
import org.mklab.nfc.matrix.ComplexRationalPolynomialMatrix;
import org.mklab.nfc.matrix.DoubleComplexMatrix;
import org.mklab.nfc.matrix.DoubleComplexPolynomialMatrix;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.matrix.DoublePolynomialMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.matrix.RealPolynomialMatrix;
import org.mklab.nfc.matrix.RealRationalPolynomialMatrix;
import org.mklab.nfc.scalar.AnyComplexPolynomial;
import org.mklab.nfc.scalar.AnyRealPolynomial;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.ComplexPolynomial;
import org.mklab.nfc.scalar.ComplexRationalPolynomial;
import org.mklab.nfc.scalar.DoubleComplexPolynomial;
import org.mklab.nfc.scalar.DoubleNumber;
import org.mklab.nfc.scalar.DoublePolynomial;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.nfc.scalar.RealPolynomial;
import org.mklab.nfc.scalar.RealRationalPolynomial;


/**
 * 特性多項式を求めるクラスです。
 * 
 * <p>Characteristic polynomial
 * 
 * @author koga
 * @version $Revision: 1.29 $
 * 
 */
public class Makepoly {

  /**
   * もし<code>A</code>が正方行列なら、<code>det(sI - A)</code>、すなわち <code>A</code>の特性多項式を返します。
   * 
   * <p> もし<code>A</code>がベクトルなら、 <code>A</code>の成分を根とする多項式を返します。
   * 
   * @param A 多項式の根(正方行列)
   * @return 特性多項式 (characteristic polynomial)
   */
  public static DoublePolynomial makepoly(DoubleMatrix A) {
    final int rowSize = A.getRowSize();
    final int columnSize = A.getColumnSize();

    if (rowSize == 0 && columnSize == 0) {
      return new DoublePolynomial(1);
    }

//    if (rowSize == 1 || columnSize == 1) {
//      return createPolynomialWithRoots(new DoubleComplexMatrix(A));
//    }

    if (rowSize != columnSize) {
      throw new IllegalArgumentException(Messages.getString("Makepoly.0")); //$NON-NLS-1$
    }

    final DoublePolynomial s = createPolynomialVariable();

    final DoubleComplexMatrix roots = A.eigenValue();
    DoublePolynomial polynomial1 = new DoubleComplexPolynomialMatrix(roots.unaryMinus()).addElementWise(new DoubleComplexPolynomial(s)).product().getRealPart(); 

    if (rowSize <= 6) {
      final DoublePolynomial polynomial2 = ((new DoublePolynomialMatrix(A.createUnit()).multiply(s).subtract(new DoublePolynomialMatrix(A)))).determinant();
      final DoubleNumber size2 = polynomial2.evaluate(A).frobNorm();
      final DoubleNumber size1 = polynomial1.evaluate(A).frobNorm();
      if (size2.isLessThan(size1)) {
        return polynomial2;
      }
    }

    return polynomial1;
  }

  /**
   * もし<code>A</code>が正方行列なら、<code>det(sI - A)</code>、すなわち <code>A</code>の特性多項式を返します。
   * 
   * <p> もし<code>A</code>がベクトルなら、 <code>A</code>の成分を根とする多項式を返します。
   * 
   * @param <RPS> 実多項式の型
   * @param <RPM> 実多項式行列の型
   * @param <CPS> 複素多項式の型
   * @param <CPM> 複素多項式行列の型
   * @param <RRS> 実有理多項式の型
   * @param <RRM> 実有理多項式行列の型
   * @param <CRS> 複素有理多項式の型
   * @param <CRM> 複素有理多項式行列の型
   * @param <RES> 実係数スカラーの型
   * @param <REM> 実係数行列の型
   * @param <CES> 複素係数スカラーの型
   * @param <CEM> 複素係数行列の型
   * @param A 多項式の根(正方行列)
   * @param rs real symbolic  
   * @return 特性多項式 (characteristic polynomial)
   */
  public static <RPS extends RealPolynomial<RPS,RPM,CPS,CPM,RRS,RRM,CRS,CRM,RES,REM,CES,CEM>, RPM extends RealPolynomialMatrix<RPS,RPM,CPS,CPM,RRS,RRM,CRS,CRM,RES,REM,CES,CEM>, CPS extends ComplexPolynomial<RPS,RPM,CPS,CPM,RRS,RRM,CRS,CRM,RES,REM,CES,CEM>, CPM extends ComplexPolynomialMatrix<RPS,RPM,CPS,CPM,RRS,RRM,CRS,CRM,RES,REM,CES,CEM>, RRS extends RealRationalPolynomial<RPS, RPM, CPS, CPM, RRS, RRM, CRS, CRM, RES, REM, CES, CEM>, RRM extends RealRationalPolynomialMatrix<RPS, RPM, CPS, CPM, RRS, RRM, CRS, CRM, RES, REM, CES, CEM>, CRS extends ComplexRationalPolynomial<RPS, RPM, CPS, CPM, RRS, RRM, CRS, CRM, RES, REM, CES, CEM>, CRM extends ComplexRationalPolynomialMatrix<RPS, RPM, CPS, CPM, RRS, RRM, CRS, CRM, RES, REM, CES, CEM>, RES extends RealNumericalScalar<RES,REM,CES,CEM>, REM extends RealNumericalMatrix<RES,REM,CES,CEM>, CES extends ComplexNumericalScalar<RES,REM,CES,CEM>, CEM extends ComplexNumericalMatrix<RES,REM,CES,CEM>> RPS makepoly(REM A, RPS rs) {
    final int rowSize = A.getRowSize();
    final int columnSize = A.getColumnSize();

    if (rowSize == 0 && columnSize == 0) {
      return rs.create(1);
    }

//    if (rowSize == 1 || columnSize == 1) {
//      return createPolynomialWithRoots(new DoubleComplexMatrix(A));
//    }

    if (rowSize != columnSize) {
      throw new IllegalArgumentException(Messages.getString("Makepoly.0")); //$NON-NLS-1$
    }

    final RPS s =  rs.create("s"); //$NON-NLS-1$

    final CEM roots = A.eigenValue();
    RPS polynomial1 =  rs.toComplex().createGrid(roots.unaryMinus()).addElementWise(s.toComplex()).product().getRealPart(); 

    if (rowSize <= 6) {
      final RPS polynomial2 = ((rs.createGrid(A.createUnit()).multiply(s).subtract(rs.createGrid(A)))).determinant();
      final RES size2 = polynomial2.evaluate(A).frobNorm();
      final RES size1 = polynomial1.evaluate(A).frobNorm();
      if (size2.isLessThan(size1)) {
        return polynomial2;
      }
    }

    return polynomial1;
  }

  
  /**
   * 多項式の変数(1次の係数のみ1)を返します。
   * 
   * @return 多項式の変数(1次の係数のみ1)
   */
  private static DoublePolynomial createPolynomialVariable() {
    return new DoublePolynomial("s"); //$NON-NLS-1$
  }

  /**
   * 与えられた根をもつ多項式を返します。
   * 
   * @param <RES> 実係数スカラーの型
   * @param <REM> 実係数行列の型
   * @param <CES> 複素係数スカラーの型
   * @param <CEM> 複素係数行列の型
   * @param matrix 根
   * @return 与えられた根をもつ多項式
   */
  public static <RES extends RealNumericalScalar<RES,REM,CES,CEM>, REM extends RealNumericalMatrix<RES,REM,CES,CEM>, CES extends ComplexNumericalScalar<RES,REM,CES,CEM>, CEM extends ComplexNumericalMatrix<RES,REM,CES,CEM>> AnyRealPolynomial<RES,REM,CES,CEM> makepoly(REM matrix) {
    final int rowSize = matrix.getRowSize();
    final int columnSize = matrix.getColumnSize();

    if (rowSize == 0 || columnSize == 0) {
      throw new IllegalArgumentException(Messages.getString("Makepoly.1")); //$NON-NLS-1$
    }
    
    CEM roots = matrix.eigenValue();

    return makepoly(roots);
  }

  /**
   * 与えられた根をもつ多項式を返します。
   * 
   * @param <RES> 実係数スカラーの型
   * @param <REM> 実係数行列の型
   * @param <CES> 複素係数スカラーの型
   * @param <CEM> 複素係数行列の型
   * @param roots 根
   * @return 与えられた根をもつ多項式
   */
  public static <RES extends RealNumericalScalar<RES,REM,CES,CEM>, REM extends RealNumericalMatrix<RES,REM,CES,CEM>, CES extends ComplexNumericalScalar<RES,REM,CES,CEM>, CEM extends ComplexNumericalMatrix<RES,REM,CES,CEM>> AnyRealPolynomial<RES,REM,CES,CEM> makepoly(CEM roots) {
    final int rowSize = roots.getRowSize();
    final int columnSize = roots.getColumnSize();

    if (rowSize == 0 || columnSize == 0) {
      throw new IllegalArgumentException(Messages.getString("Makepoly.1")); //$NON-NLS-1$
    }

    if (rowSize != 1 && columnSize != 1) {
      throw new IllegalArgumentException(Messages.getString("Makepoly.1")); //$NON-NLS-1$
    }
    
    final CEM columnRoots = org.mklab.nfc.matrix.util.Makecolv.makecolv(roots);
    
    final AnyRealPolynomial<RES,REM,CES,CEM> s =  new AnyComplexPolynomial<>(roots).getRealPart().create("s"); //$NON-NLS-1$
    final AnyRealPolynomial<RES,REM,CES,CEM> polynomial =  s.toComplex().createGrid(columnRoots.unaryMinus()).addElementWise(s.toComplex()).product().getRealPart();

    return polynomial;
  }

  /**
   * 与えられた根をもつ多項式を返します。
   * 
   * @param <RPS> 実多項式の型
   * @param <RPM> 実多項式行列の型
   * @param <CPS> 複素多項式の型
   * @param <CPM> 複素多項式行列の型
   * @param <RRS> 実有理多項式の型
   * @param <RRM> 実有理多項式行列の型
   * @param <CRS> 複素有理多項式の型
   * @param <CRM> 複素有理多項式行列の型
   * @param <RES> 実係数スカラーの型
   * @param <REM> 実係数行列の型
   * @param <CES> 複素係数スカラーの型
   * @param <CEM> 複素係数行列の型
   * @param roots 根
   * @param rs real polynomial 
   * @return 与えられた根をもつ多項式
   */
  public static <RPS extends RealPolynomial<RPS,RPM,CPS,CPM,RRS,RRM,CRS,CRM,RES,REM,CES,CEM>, RPM extends RealPolynomialMatrix<RPS,RPM,CPS,CPM,RRS,RRM,CRS,CRM,RES,REM,CES,CEM>, CPS extends ComplexPolynomial<RPS,RPM,CPS,CPM,RRS,RRM,CRS,CRM,RES,REM,CES,CEM>, CPM extends ComplexPolynomialMatrix<RPS,RPM,CPS,CPM,RRS,RRM,CRS,CRM,RES,REM,CES,CEM>, RRS extends RealRationalPolynomial<RPS, RPM, CPS, CPM, RRS, RRM, CRS, CRM, RES, REM, CES, CEM>, RRM extends RealRationalPolynomialMatrix<RPS, RPM, CPS, CPM, RRS, RRM, CRS, CRM, RES, REM, CES, CEM>, CRS extends ComplexRationalPolynomial<RPS, RPM, CPS, CPM, RRS, RRM, CRS, CRM, RES, REM, CES, CEM>, CRM extends ComplexRationalPolynomialMatrix<RPS, RPM, CPS, CPM, RRS, RRM, CRS, CRM, RES, REM, CES, CEM>, RES extends RealNumericalScalar<RES,REM,CES,CEM>, REM extends RealNumericalMatrix<RES,REM,CES,CEM>, CES extends ComplexNumericalScalar<RES,REM,CES,CEM>, CEM extends ComplexNumericalMatrix<RES,REM,CES,CEM>> RPS  makepoly(CEM roots, RPS rs) {
    final int rowSize = roots.getRowSize();
    final int columnSize = roots.getColumnSize();

    if (rowSize == 0 || columnSize == 0) {
      throw new IllegalArgumentException(Messages.getString("Makepoly.1")); //$NON-NLS-1$
    }

    if (rowSize != 1 && columnSize != 1) {
      throw new IllegalArgumentException(Messages.getString("Makepoly.1")); //$NON-NLS-1$
    }
    
    final CEM columnRoots = org.mklab.nfc.matrix.util.Makecolv.makecolv(roots);
    
    final RPS s =  rs.create("s"); //$NON-NLS-1$
    final RPS polynomial =  s.toComplex().createGrid(columnRoots.unaryMinus()).addElementWise(s.toComplex()).product().getRealPart();

    return polynomial;
  }

  /**
   * 与えられた根をもつ多項式を返します。
   * 
   * @param roots 根
   * @return 与えられた根をもつ多項式
   */
  public static DoublePolynomial makepoly(DoubleComplexMatrix roots) {
    final int rowSize = roots.getRowSize();
    final int columnSize = roots.getColumnSize();

    if (rowSize == 0 || columnSize == 0) {
      return new DoublePolynomial(1);
    }

    if (rowSize != 1 && columnSize != 1) {
      throw new IllegalArgumentException(Messages.getString("Makepoly.1")); //$NON-NLS-1$
    }

    final DoubleComplexMatrix columnRoots = Makecolv.makecolv(roots);

    final DoublePolynomial s = createPolynomialVariable();
    final DoublePolynomial polynomial =  new DoubleComplexPolynomialMatrix(columnRoots.unaryMinus()).addElementWise(new DoubleComplexPolynomial(s)).product().getRealPart();

    return polynomial;
  }

}