HornorVerifier.java

/*
 * Created on 2009/02/06
 * Copyright (C) 2009 Koga Laboratory. All rights reserved.
 *
 */
package org.mklab.cga.polynomial;

import org.mklab.cga.interval.matrix.IntervalNumericalMatrix;
import org.mklab.cga.interval.scalar.IntervalNumericalScalar;
import org.mklab.nfc.matrix.NumericalMatrix;
import org.mklab.nfc.scalar.NumericalScalar;
import org.mklab.nfc.util.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/**
 * 多項式の誤差評価を行うクラスです。
 * 
 * 多項式は、ホーナー法によって解かれます。 つまりこのクラスではホーナー法の誤差評価を行います。
 * 
 * @author yano
 * @version $Revision$, 2009/02/06
 * @param <IS> 区間スカラーの型
 * @param <IM> 区間行列の型
 * @param <S> 成分の型
 * @param <M> 行列の型
 */
public class HornorVerifier<IS extends IntervalNumericalScalar<IS,IM,S,M>, IM extends IntervalNumericalMatrix<IS,IM,S,M>, S extends NumericalScalar<S,M>, M extends NumericalMatrix<S,M>> {
  /** vector of coefficients */
  private M A;
  /** value for evaluation */
  private S v;
  /** interval for evaluation */
  private IS type;
  /** solution */
  private IS solution;
  
  /**
   * Creates {@link HornorVerifier}.
   * @param A 多項式の係数行列
   * @param p 多項式に代入する値
   * @param type 多項式に代入する値の点区間
   */
  public HornorVerifier(final M A, final S p, final IS type) {
    this.A = A;
    this.v = p;
    this.type = type;
  }

  /**
   * 行列<code>A</code>を係数としてもつ多項式に<code>x</code>を代入したときの値の誤差評価を行います。
   */
  public void  solve() {
    int size = this.A.getColumnSize();
    S unit = this.A.getElement(1, 1);
    S[] a = unit.createArray(size);
    for (int i = 0; i < size; i++) {
      a[i] = this.A.getElement(1, i + 1);
    }

    this.solution = solve(a, this.v);

//    if (A.getElement(1).isZero()) {
//      a = unit.createArray(size - 1);
//      for (int i = 0; i < size - 1; i++) {
//        a[i] = A.getElement(1, i + 2);
//      }
//    } else {
//      a = unit.createArray(size);
//      for (int i = 0; i < size; i++) {
//        a[i] = A.getElement(1, i + 1);
//      }
//    }
  }
  
  /**
   * Returns solution.
   * 
   * @return solution
   */
  public IS getSolution() {
    return this.solution;
  }

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

    int size = coefficient.length;
    S unit = coefficient[0];
    S[] s = unit.createArray(size);
    S[] p = unit.createArray(size);
    S[] r = unit.createArray(size);
    S[] q = unit.createArray(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] = x.createZero();
    }

    S t = x.clone();

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

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    S t2 = x.subtract(x.abs().multiply(x.getMachineEpsilon()));

    manager.setRoundMode(RoundMode.ROUND_UP);

    S t1 = x.add(x.abs().multiply(x.getMachineEpsilon()));

    S tc = t1.add(t2).divide(2.0);
    S d = tc.subtract(t2);
    S ta = t1.abs();
    S tn = ta.createUnit();

    r[size - 1] = p[size - 1].subtract(s[size - 1]);

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

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    q[size - 1] = p[size - 1].subtract(s[size - 1]);

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

    S tb = ta.createUnit().subtract(ta);

    manager.setRoundMode(RoundMode.ROUND_UP);

    S g = d.multiply(tn.createUnit().subtract(tn)).divide(tb);
    S a = (tn.createUnit().subtract(tn.multiply(ta))).divide(tb);

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    S gg = g.createUnit().subtract(g);

    manager.setRoundMode(RoundMode.ROUND_UP);

    a = a.divide(gg);

    S pl = a.multiply(q[0]);
    S sl = s[0].add(pl);

    S pu = a.multiply(r[0]);
    S su = s[0].add(pu);

    IS ans = this.type.createInfSup(sl, su);

    manager.setRoundMode(oldRoundMode);
    return ans;
  }
  
  /**
   * Returns A. 
   * @return A
   */
  protected M getA() {
    return this.A;
  }
  
  /**
   * Returns v.
   * 
   * @return v
   */
  protected S getV() {
    return this.v;
  }
  
  /**
   * Returns iv.
   * 
   * @return iv
   */
  protected IS getIV() {
    return this.type;
  }
}