LogVerifier.java

/*
 * Created on 2004/12/20
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 *
 */
package org.mklab.cga.interval.function;

import org.mklab.cga.interval.matrix.DoubleIntervalMatrix;
import org.mklab.cga.interval.matrix.IntervalNumericalMatrix;
import org.mklab.cga.interval.scalar.DoubleIntervalNumber;
import org.mklab.cga.interval.scalar.IntervalNumericalScalar;
import org.mklab.cga.polynomial.HornorVerifier;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.matrix.NumericalMatrix;
import org.mklab.nfc.scalar.DoubleNumber;
import org.mklab.nfc.scalar.NumericalScalar;
import org.mklab.nfc.util.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/**
 * <code>log(t)</code> を精度保証付きで計算するプログラムです。 未完成。
 * 
 * @author hiroki
 * @version $Revision: 1.10 $.2004/12/20
 * @param <IS> 区間スカラーの型
 * @param <IM> 区間行列の型
 * @param <S> 成分の型
 * @param <M> 行列の型
 */
public class LogVerifier<IS extends IntervalNumericalScalar<IS,IM,S,M>, IM extends IntervalNumericalMatrix<IS,IM,S,M>, S extends NumericalScalar<S,M>, M extends NumericalMatrix<S,M>> {

  /** */
  private int n;

  /**
   * コンストラクタ
   * 
   * @param n 対象となる値
   */
  public LogVerifier(int n) {
    if (n % 2 == 0) {
      System.err.println("奇数を設定してください"); //$NON-NLS-1$
      throw new RuntimeException("nには奇数値を設定してください。"); //$NON-NLS-1$
    } else if (n == 3 || n == 5 || n == 7) {//3, 5, 7に関しては打切り誤差も評価
      this.n = n + 1;
    } else {//打切り誤差を評価していない場合
      System.err.println("打切り誤差は無視されています。"); //$NON-NLS-1$
    }
  }

  /**
   * 対数を精度保証付きで求めます。
   * 
   * @param t 対数を求めたい値
   * @return 対数
   */
  public DoubleIntervalNumber solve(double t) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    //log(t)=log(s)+m*log(2) 1 <= s <=2
    int m = 0;
    double s = t;

    //sは[1, 2]の数 t = 2^m*sとあらわす
    while (s >= 2.0) {
      m++;
      s /= 2.0;
    }

    //ckのリストの生成
    double[] ck = new double[64];

    for (int i = 0; i < 64; i++) {
      ck[i] = 1 + i / 64.0;
    }

    //|s-ck| < 1/128となるようなckを探す
    int ck_index = 0;
    for (int i = 0; i < 64; i++) {
      if (Math.abs(s - ck[i]) <= 1 / 128.0) {
        ck_index = i;
        break;
      }
    }

    double r = 2 * (s - ck[ck_index]) / (s + ck[ck_index]);

    DoubleNumber[] pp = new DoubleNumber[this.n + 1];
    //テイラー展開の係数を作成
    for (int i = 0; i < this.n; i++) {
      if (i % 2 == 0) {
        pp[i] = new DoubleNumber(0);
      } else {
        pp[i] = new DoubleNumber(1.0 / (Math.pow(-2, (this.n - i) - 1) * (this.n - i)));
      }
    }

    //log(ck)の計算
    pp[this.n] = new DoubleNumber(Math.log(ck[ck_index]));
    //ホーナー法

    //精度保証1
    manager.setRoundMode(RoundMode.ROUND_DOWN);
    double l = Math.log(2) * m;

    manager.setRoundMode(RoundMode.ROUND_UP);
    double u = Math.log(2) * m;

    DoubleMatrix pm = new DoubleMatrix(pp);
    HornorVerifier<DoubleIntervalNumber, DoubleIntervalMatrix,DoubleNumber,DoubleMatrix> verifer = new HornorVerifier<>(pm, new DoubleNumber(r), new DoubleIntervalNumber(r).add(new DoubleIntervalNumber(l, u)));
    verifer.solve();
    DoubleIntervalNumber log = verifer.getSolution();

    if (this.n == 4) {
      manager.setRoundMode(RoundMode.ROUND_DOWN);
      double inf = log.getInfimum().doubleValue() - (3.7e-13);

      manager.setRoundMode(RoundMode.ROUND_UP);
      double sup = log.getInfimum().doubleValue() + (3.7e-13);

      manager.setRoundMode(oldRoundMode);
      return new DoubleIntervalNumber(inf, sup);
    }
    if (this.n == 6) {
      manager.setRoundMode(RoundMode.ROUND_DOWN);
      double inf = log.getInfimum().doubleValue() - (4.0e-18);

      manager.setRoundMode(RoundMode.ROUND_UP);
      double sup = log.getInfimum().doubleValue() + (4.0e-18);

      manager.setRoundMode(oldRoundMode);
      return new DoubleIntervalNumber(inf, sup);
    }
    if (this.n == 8) {
      manager.setRoundMode(RoundMode.ROUND_DOWN);
      double inf = log.getInfimum().doubleValue() - (4.8e-23);

      manager.setRoundMode(RoundMode.ROUND_UP);
      double sup = log.getInfimum().doubleValue() + (4.8e-23);

      manager.setRoundMode(oldRoundMode);
      return new DoubleIntervalNumber(inf, sup);
    }

    DoubleIntervalNumber ans = log;
    manager.setRoundMode(oldRoundMode);
    return ans;
  }

  /**
   * 対数を精度保証付きで求めます。
   * 
   * @param t 対数を求めたい値。
   * @return 対数。
   */
  public IS solve(final IS t) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    S unit = t.getInfimum().createUnit();

    //精度保証2
    int m = 0;

    // Interval i_t = t;
    IS i_s = t;
    IS i_a = t.create(2); // new DoubleIntervalNumber(2.0);//2で割る
    while ((i_s.getSupremum()).isGreaterThan(2)) {
      m++;
      i_s = i_s.divide(i_a);
    }

    IS[] i_ck = t.createArray(64); //new DoubleIntervalNumber[64];
    for (int i = 0; i < 64; i++) {
      i_ck[i] = t.create(i).divide(64).add(1); // new DoubleIntervalNumber(1 + i / 64.0);
    }

    //|s-ck| < 1/128となるようなckを探す
    int ck_index = 0;
    for (int i = 0; i < 64; i++) {
      if ((i_s.getMiddle().subtract(i_ck[i].getMiddle())).abs().isLessThanOrEquals(1 / 128.0)) {
        ck_index = i;
        break;
      }
    }

    IS i_r = t.create(2).multiply((i_s.subtract(i_ck[ck_index])).divide(i_s.add(i_ck[ck_index])));
    IS[] i_p = t.createArray(this.n); // new DoubleIntervalNumber[this.n];
    //テイラー展開係数の作成
    for (int i = 0; i < this.n; i++) {
      if (i % 2 == 0) {
        i_p[i] = t.create(0); //new DoubleIntervalNumber(0);
      } else {
        i_p[i] = t.create(1).divide(t.create(-2).power(i - 1).multiply(i)); // Math.pow(-2, i - 1) * i)
      }
    }

    //i_p[0] = new DoubleIntervalNumber(Math.log(((DoubleNumber)i_ck[ck_index].getInfimum()).doubleValue()), Math.log(((DoubleNumber)i_ck[ck_index].getSupremum()).doubleValue()));
    i_p[0] = t.createInfSup((i_ck[ck_index].getInfimum()).log(), (i_ck[ck_index].getSupremum()).log());

    i_s = horner(i_p, i_r);

    //精度保証1
    manager.setRoundMode(RoundMode.ROUND_DOWN);
    S l = unit.create(2).log().multiply(m); // Math.log(2) * m;

    manager.setRoundMode(RoundMode.ROUND_UP);
    S u = unit.create(2).log().multiply(m); // Math.log(2) * m;

    IS log = i_s.add(t.createInfSup(l, u));

    if (this.n == 4) {
      manager.setRoundMode(RoundMode.ROUND_DOWN);
      S inf = log.getInfimum().subtract(3.7e-13);

      manager.setRoundMode(RoundMode.ROUND_UP);
      S sup = log.getInfimum().add(3.7e-13);

      manager.setRoundMode(oldRoundMode);
      return t.createInfSup(inf, sup);
    }
    if (this.n == 6) {
      manager.setRoundMode(RoundMode.ROUND_DOWN);
      S inf = log.getInfimum().subtract(4.0e-18);

      manager.setRoundMode(RoundMode.ROUND_UP);
      S sup = log.getInfimum().add(4.0e-18);

      manager.setRoundMode(oldRoundMode);
      return t.createInfSup(inf, sup);
    }
    if (this.n == 8) {
      manager.setRoundMode(RoundMode.ROUND_DOWN);
      S inf = log.getInfimum().subtract(4.8e-23);

      manager.setRoundMode(RoundMode.ROUND_UP);
      S sup = log.getInfimum().add(4.8e-23);

      manager.setRoundMode(oldRoundMode);
      return t.createInfSup(inf, sup);
    }

    IS ans = log;
    manager.setRoundMode(oldRoundMode);
    return ans;
  }

  //  /**
  //   * @param p
  //   * @param t
  //   * @return ?
  //   */
  //  private double horner(double[] p, double t) {
  //    double[] s = new double[p.length];
  //    int i;
  //    s[n - 1] = p[n - 1];
  //
  //    for (int j = 0; j < n - 1; j++) {
  //      i = n - 2 - j;
  //      s[i] = s[i + 1] * t + p[i];
  //    }
  //
  //    return s[0];
  //  }

  /**
   * @param p p
   * @param t t
   * @return result
   */
  private IS horner(IS[] p, IS t) {
    IS unit = p[0];
    IS[] s = unit.createArray(p.length); //new DoubleIntervalNumber[p.length];
    int i;
    s[this.n - 1] = p[this.n - 1];

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

    return s[0];
  }

}