DKA.java

/*
 * Created on 2004/03/11
 *
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.cga.polynomial;

import org.mklab.nfc.matrix.DoubleComplexMatrix;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.scalar.DoubleComplexNumber;
import org.mklab.nfc.scalar.DoubleNumber;
import org.mklab.nfc.util.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/**
 * 多項式の根の近似解を求めるクラス。 ここではDurand-Kerner-Aberth法を用いて計算している。
 * 
 * @author Hiroki
 */
public class DKA {

  /** 設定した行列 */
  private DoubleMatrix matrix;
  /** 設定した行列の要素数 */
  private int length;
  /** 多項式の次数 */
  private int order;
  /** 複素数-i */
  private DoubleComplexNumber complex_i = new DoubleComplexNumber(0, 1);
  /** 設定した行列の要素 */
  private DoubleComplexNumber[] matrix_elements;
  /** 代数方程式の解の要素 */
  private DoubleComplexNumber[] solution_elements;

  /**
   * 多項式の係数を最も大きい次数から順にMatrix型で設定します。
   * 
   * @param m 任意の行列。
   */
  public DKA(final DoubleMatrix m) {
    this.matrix = lotate(m);
  }

  /**
   * Durand-Kerner-Aberthの方法を用いて代数方程式の 近似解を求めます。
   * 
   * @param errorSpan 許容誤差。
   * @return 近似解。
   */
  public DoubleComplexMatrix dka(final double errorSpan) {
    double es = errorSpan;
    this.length = this.matrix.count();
    this.matrix_elements = new DoubleComplexNumber[this.length];
    this.solution_elements = new DoubleComplexNumber[this.length];
    this.order = this.length - 1;

    DoubleComplexNumber[] dx = new DoubleComplexNumber[this.order];
    DoubleComplexNumber[] newx = new DoubleComplexNumber[this.order];

    System.out.println(Messages.getString("DKA.0") + this.order); //$NON-NLS-1$

    for (int i = 0; i < this.length; i++) {
      this.matrix_elements[i] = new DoubleComplexNumber((this.matrix).getDoubleElement(1, this.length - i), 0);
      System.out.println(Messages.getString("DKA.1") + this.matrix_elements[i]); //$NON-NLS-1$
    }

    // Aberthの初期値を配置する円の決定
    DoubleComplexNumber g = this.matrix_elements[1].unaryMinus().divide(this.order);
    System.out.println(Messages.getString("DKA.2") + g); //$NON-NLS-1$

    double max = 0, r0, r;
    for (int i = 1; i <= this.order; i++) {
      if (this.matrix_elements[i].abs().getRealPart().doubleValue() > max) {
        max = this.matrix_elements[i].abs().getRealPart().doubleValue();
      }
    }

    System.out.println("max|a_i|: " + max); //$NON-NLS-1$

    r0 = g.abs().getRealPart().doubleValue() + 1 + max;

    r = r0;
    if (r0 > r) {
      r = r0;
    }

    // Aberthの初期値
    for (int i = 0; i < this.order; i++) {
      double theta = 2 * i * Math.PI / this.order + Math.PI / (2 * this.order);
      this.solution_elements[i] = exp(this.complex_i.multiply(theta)).multiply(r0).add(g); // ここまで3/11
    }

    // 反復(Durand-Kerner法)
    // zk(v+1) = zk - P(zk(v) / Π(zk(v) - zj(v)) (1 <= k <= n) (v = 0, 1, 2,
    // ....)
    for (int k = 0; k < 1000; k++) {
      for (int i = 0; i < this.order; i++) {
        DoubleComplexNumber polynomial = polynomial(this.order, this.matrix_elements, this.solution_elements[i]);
        DoubleComplexNumber bunbo = bunbo(this.order, this.solution_elements, i);
        dx[i] = polynomial.divide(bunbo);
        newx[i] = this.solution_elements[i].subtract(dx[i]);
      }
      // 更新
      for (int i = 0; i < this.order; i++) {
        this.solution_elements[i] = newx[i];
      }

      // 変化量を計算する
      double err = 0.0;
      for (int i = 0; i < this.order; i++) {
        err += dx[i].abs().getRealPart().doubleValue();
      }

      // 変化量が小さくなれば収束したと判断して反復を終了します。
      if (err < es) {
        break;
      }
    }

    DoubleComplexMatrix kai = new DoubleComplexMatrix(this.order, 1);

    for (int i = 0; i < this.order; i++) {
      kai.setElement(i + 1, 1, this.solution_elements[i]);
    }
    return kai;
  }

//  /**
//   * コンパニオン行列の固有値を求めることで 多項式の根の近似解を求めます。 固有値の精度保証を用いて根の計算の精度保証をします。 この方法はたぶんよくない。
//   * 
//   * @return 近似解と誤差。
//   */
//  public DoubleComplexIntervalMatrix getRootErrorWithEigenError() {
//    DoubleMatrix m = CompanionMatrix.create(this.matrix.getSubMatrix(1, 1, 1, this.matrix.getColumnSize() - 1));
//    RealEigenVerifier<DoubleIntervalNumber,DoubleIntervalMatrix,DoubleComplexIntervalNumber,DoubleComplexIntervalMatrix,DoubleNumber,DoubleMatrix,DoubleComplexNumber,DoubleComplexMatrix> verifier = new RealRumpMethod<>();
//    verifier.solve(new DoubleIntervalMatrix(m));
//    verifier.getEigenValue();
//  }

//  public Matrix[] getRootErrorWithEigenSmith() { 
//    Matrix m = matrix.getSubMatrix(1, 1, 1, matrix.getColSize() - 1).getCompanionMatrix();
//    m.print("コンパニオン行列"); 
//    Eigen e = new Eigen(m); 
//  }

  /**
   * Durand-Kerner-Aberthの方法を用いて多項式の 近似解を求め、Smithの定理によって近似解の誤差 を計算します。
   * 
   * @return 近似解と誤差。
   */
  public DoubleComplexMatrix[] getRootErrorWithSmith() {
    return smith();
  }

  /**
   * Smithの定理を用いて代数方程式の近似解と誤差を返します。
   * 
   * @return 近似解とその誤差。
   */
  public DoubleComplexMatrix[] smith() {
    DoubleComplexMatrix solution = dka(1e-15);
    DoubleMatrix error = getError(this.order, this.solution_elements, this.matrix_elements);

    DoubleComplexMatrix[] answer = {solution, error.toComplex()};
    return answer;
  }

  /**
   * 複素数の指数関数を求めます。
   * 
   * @param c 任意の複素数。
   * @return cの指数関数。
   */
  private DoubleComplexNumber exp(final DoubleComplexNumber c) {
    DoubleComplexNumber ans = (this.complex_i.multiply(Math.sin(c.getImaginaryPart().doubleValue())).add(Math.cos(c.getImaginaryPart().doubleValue()))).multiply(Math
        .exp(c.getRealPart().doubleValue()));
    return ans;
  }

  /**
   * Horner法を用いて多項式の値の計算を行います。
   * 
   * @param polynomialOrder 次数。
   * @param a 設定した行列の要素列。
   * @param x 設定した多項式の解。
   * @return 計算結果。
   */
  private DoubleComplexNumber polynomial(final int polynomialOrder, final DoubleComplexNumber[] a, final DoubleComplexNumber x) {
    DoubleComplexNumber w = a[0];

    for (int i = 1; i <= polynomialOrder; i++) {
      w = (w.multiply(x).add(a[i]));
    }
    return w;
  }

  /**
   * zk(v+1) = zk - P(zk(v) / Π(zk(v) - zj(v)) (1 &lt;= k &lt;= n) (v = 0, 1, 2, ....) のΠ(zk(v) - zj(v))の計算を行います。
   * 
   * @param polynomialOrder 次数。
   * @param x 設定した多項式の解。
   * @param i 反復数。
   * @return 計算結果。
   */
  private DoubleComplexNumber bunbo(final int polynomialOrder, final DoubleComplexNumber[] x, final int i) {
    DoubleComplexNumber w = new DoubleComplexNumber(1, 0);
    for (int j = 0; j < polynomialOrder; j++) {
      if (j != i) {
        w = w.multiply(x[i].subtract(x[j]));
      }
    }
    return w;
  }

  /**
   * 代数方程式の誤差を計算する部分。
   * 
   * @param polynomialOrder 次数。
   * @param x 設定した多項式の解。
   * @param a 多項式の要素列。
   * @return 誤差。
   */
  private DoubleMatrix getError(final int polynomialOrder, final DoubleComplexNumber[] x, final DoubleComplexNumber[] a) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    DoubleComplexNumber polynomial = null;
    DoubleComplexNumber bunbo = null;
    DoubleComplexMatrix polynomial_up = new DoubleComplexMatrix(polynomialOrder, 1);
    DoubleComplexMatrix polynomial_down = new DoubleComplexMatrix(polynomialOrder, 1);
    DoubleComplexMatrix bunbo_up = new DoubleComplexMatrix(polynomialOrder, 1);
    DoubleComplexMatrix bunbo_down = new DoubleComplexMatrix(polynomialOrder, 1);

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    for (int i = 0; i < polynomialOrder; i++) {
      bunbo = bunbo(polynomialOrder, x, i);
      polynomial = polynomial(polynomialOrder, a, x[i]);
      bunbo_up.setElement(i + 1, 1, bunbo);
      polynomial_up.setElement(i + 1, 1, polynomial);
    }

    manager.setRoundMode(RoundMode.ROUND_UP);

    for (int i = 0; i < polynomialOrder; i++) {
      bunbo = bunbo(polynomialOrder, x, i);
      polynomial = polynomial(polynomialOrder, a, x[i]);
      bunbo_down.setElement(i + 1, 1, bunbo);
      polynomial_down.setElement(i + 1, 1, polynomial);
    }

    DoubleMatrix polynomial_max = (polynomial_up.absElementWise()).maxElementWise(polynomial_down.absElementWise()).getRealPart();
    DoubleMatrix bunbo_min = (bunbo_up.absElementWise()).minElementWise(bunbo_down.absElementWise()).getRealPart();

    DoubleMatrix error = polynomial_max.divideElementWise(bunbo_min).multiply(polynomialOrder);

    manager.setRoundMode(oldRoundMode);
    return error;
  }

  /**
   * 行列の要素を入れ替え、最大次数の係数を1にします。
   * 
   * @param m 任意行列。
   * @return mの要素を入れ替え、最大次数の係数を1にした行列。
   */
  private DoubleMatrix lotate(final DoubleMatrix m) {
    this.matrix = m.flipLeftRight();

    for (int i = this.matrix.getColumnSize(); i > 0; i--) {
      DoubleNumber c = this.matrix.getElement(1, i);
      if (c.isUnit()) {
        break;
      } else if (c.isZero() == false) {
        DoubleNumber cc = c.inverse();
        this.matrix = this.matrix.multiply(cc);
        break;
      }
    }

    return this.matrix;
  }

  /**
   * メインメソッド
   * 
   * @param args コマンドライン引数
   */
  public static void main(String[] args) {
    double time_start = System.currentTimeMillis();

    for (int i = 0; i < 1; i++) {
      DoubleMatrix A = new DoubleMatrix(new double[][] {{1, 3, 2, 6}});
      A.print("A"); //$NON-NLS-1$
      DKA dk = new DKA(A);
      dk.dka(1e-15).print("dka"); //$NON-NLS-1$
      double time1 = System.currentTimeMillis();
      DoubleComplexMatrix[] answer1 = dk.getRootErrorWithSmith();
      double time2 = System.currentTimeMillis();
      answer1[0].print(Messages.getString("DKA.3")); //$NON-NLS-1$
      answer1[1].getRealPart().print(Messages.getString("DKA.4")); //$NON-NLS-1$
      System.out.println(Messages.getString("DKA.5") + (time2 - time1)); //$NON-NLS-1$
    }

    double time_end = System.currentTimeMillis();
    System.out.println(Messages.getString("DKA.6") + (time_end - time_start)); //$NON-NLS-1$
  }
}