DoubleHornorVerifier.java

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

import org.mklab.cga.interval.scalar.DoubleIntervalNumber;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.util.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/**
 * 倍精度で多項式の誤差評価を行うクラスです。
 * 
 * 多項式は、ホーナー法によって解かれます。 つまりこのクラスではホーナー法の誤差評価を行います。
 * 
 * @author hiroki
 * @version $Revision: 1.14 $.2004/11/25
 */
public class DoubleHornorVerifier {

  /**
   * コンストラクタ
   * 
   */
  public DoubleHornorVerifier() {
    //
  }

  /**
   * 行列<code>A</code>を係数としてもつ多項式に<code>x</code>を代入したときの値の誤差評価を行います。
   * 
   * 
   * @param A 多項式の係数行列
   * @param x 多項式に代入する値
   * @return 評価結果
   */
  public DoubleIntervalNumber solve(final DoubleMatrix A, final double x) {
    int size = A.getColumnSize();
    double[] a;

    if (A.getDoubleElement(1) == 0) {
      a = new double[size - 1];
      for (int i = 0; i < size - 1; i++) {
        a[i] = A.getDoubleElement(1, i + 2);
      }
    } else {
      a = new double[size];
      for (int i = 0; i < size; i++) {
        a[i] = A.getDoubleElement(1, i + 1);
      }
    }
    return solve(a, x);
  }

  /**
   * 配列<code>coefficient</code>を係数としてもつ多項式に<code>x</code>を代入したときの値の誤差評価を行います。
   * 
   * @param coefficient 係数の配列
   * @param x 多項式に代入する値
   * @return 評価結果
   */
  public DoubleIntervalNumber solve(final double[] coefficient, final double x) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    int size = coefficient.length;
    double[] s = new double[size];
    double[] p = new double[size];
    double[] r = new double[size];
    double[] q = new double[size];

    // N = x-1.0/10;

    for (int n = 0; n < size; n++) {
      p[n] = coefficient[(size - 1) - n];
    }

    for (int i = 0; i < size; i++) {
      s[i] = 0;
    }

    double t = x;

    s[size - 1] = p[size - 1];
    for (int j = 0; j < size - 1; j++) {
      int i = (size - 2) - j;
      s[i] = s[i + 1] * t + p[i];
    }

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    double t2 = x - Double.MIN_VALUE * Math.abs(x);

    manager.setRoundMode(RoundMode.ROUND_UP);

    double t1 = x + Double.MIN_VALUE * Math.abs(x);

    double tc = (t1 + t2) / 2.0;
    double d = tc - t2;
    double ta = Math.abs(t1);
    double tn = 1.0;

    r[size - 1] = p[size - 1] - s[size - 1];

    for (int j = 0; j < size - 1; j++) {
      int i = (size - 2) - j;
      r[i] = s[i + 1] * tc + p[i] - s[i] + d * Math.abs(s[i + 1]);
      tn *= ta;
    }

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    q[size - 1] = p[size - 1] - s[size - 1];

    for (int j = 0; j < size - 1; j++) {
      int i = (size - 2) - j;
      q[i] = s[i + 1] * tc + p[i] - s[i] + (-d) * Math.abs(s[i + 1]);
    }

    double tb = 1.0 - ta;

    manager.setRoundMode(RoundMode.ROUND_UP);

    double g = d * (1.0 - tn) / tb;
    double a = (1.0 - tn * ta) / tb;

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    double gg = 1.0 - g;

    manager.setRoundMode(RoundMode.ROUND_UP);

    a /= gg;

    double pl = a * q[0];
    double sl = s[0] + pl;

    double pu = a * r[0];
    double su = s[0] + pu;

    DoubleIntervalNumber ans = new DoubleIntervalNumber(sl, su);

    manager.setRoundMode(oldRoundMode);
    return ans;
  }
}