IntervalUtil.java

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

import org.mklab.cga.interval.matrix.IntervalComplexNumericalMatrix;
import org.mklab.cga.interval.matrix.IntervalMatrix;
import org.mklab.cga.interval.matrix.IntervalNumericalMatrix;
import org.mklab.cga.interval.matrix.IntervalRealNumericalMatrix;
import org.mklab.cga.interval.scalar.DoubleComplexIntervalNumber;
import org.mklab.cga.interval.scalar.DoubleIntervalNumber;
import org.mklab.cga.interval.scalar.IntervalComplexNumericalScalar;
import org.mklab.cga.interval.scalar.IntervalNumericalScalar;
import org.mklab.cga.interval.scalar.IntervalRealNumericalScalar;
import org.mklab.cga.interval.scalar.IntervalScalar;
import org.mklab.nfc.matrix.BooleanMatrix;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.IntMatrix;
import org.mklab.nfc.matrix.Matrix;
import org.mklab.nfc.matrix.NormType;
import org.mklab.nfc.matrix.NumericalMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.DoubleComplexNumber;
import org.mklab.nfc.scalar.NumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.nfc.scalar.Scalar;
import org.mklab.nfc.util.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/**
 * 区間に関するユーティリティークラスです。
 * 
 * @author hiroki
 */
public final class IntervalUtil {

  /**
   * コンストラクタ。
   */
  private IntervalUtil() {
    //
  }

  /**
   * 中心を求めます。
   * 
   * @param <IS> 区間スカラーの型
   * @param <IM> 区間行列の型
   * @param <S> 成分の型
   * @param <M> 行列の型
   * @param a 区間行列
   * @return 中心
   */
  public static <IS extends IntervalScalar<IS, IM, S, M>, IM extends IntervalMatrix<IS, IM, S, M>, S extends Scalar<S, M>, M extends Matrix<S, M>> M mid(final IM a) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    manager.setRoundMode(RoundMode.ROUND_UP);
    M mid = a.getSupremum().add(a.getInfimum().subtract(a.getSupremum()));

    manager.setRoundMode(oldRoundMode);
    return mid;
  }

//  /**
//   * 区間aが区間b内に含まれるかどうかを判定します。
//   * 
//   * @param <IS> 区間スカラーの型
//   * @param <IM> 区間行列の型
//   * @param <S> 成分の型
//   * @param <M> 行列の型
//   * @param a 区間。
//   * @param b 区間。
//   * @return 含まれていればtrue, そうでなければfalse。
//   */
//  public static <IS extends IntervalNumericalScalar<IS, IM, S, M>, IM extends IntervalNumericalMatrix<IS, IM, S, M>, S extends NumericalScalar<S, M>, M extends NumericalMatrix<S, M>> BooleanMatrix in0(
//      final IM a, final IM b) {
//    RoundModeManager manager = RoundModeManager.getManager();
//    RoundMode oldRoundMode = manager.getRoundMode();
//
//    BooleanMatrix ans = new BooleanMatrix(a.getRowSize(), a.getColumnSize());
//    for (int i = 0; i < a.getRowSize(); i++) {
//      for (int j = 0; j < a.getColumnSize(); j++) {
//        IS bb = b.getElement(i + 1, j + 1);
//        IS aa = a.getElement(i + 1, j + 1);
//        boolean res = (bb.getInfimum().isLessThan(aa.getInfimum())) & (aa.getSupremum().isLessThan(bb.getSupremum()));
//        ans.setElement(i + 1, j + 1, res);
//      }
//    }
//
//    manager.setRoundMode(oldRoundMode);
//    return ans;
//  }

  /**
   * 区間aが区間b内に含まれるかどうかを判定します。
   * 
   * @param <RIS> 実区間スカラーの型
   * @param <RIM> 実区間行列の型
   * @param <CIS> 複素区間スカラーの型
   * @param <CIM> 複素区間行列の型
   * @param <RS> 実スカラーの型
   * @param <RM> 実行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @param a 区間。
   * @param b 区間。
   * @return 含まれていればtrue, そうでなければfalse。
   */
  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> boolean in0(
      final RIS a, final RIS b) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    boolean res = (b.getInfimum().isLessThan(a.getInfimum())) & (a.getSupremum().isLessThan(b.getSupremum()));
    manager.setRoundMode(oldRoundMode);
    return res;
  }

  /**
   * 区間aが区間b内に含まれるかどうかを判定します。
   * 
   * @param <RIS> 実区間スカラーの型
   * @param <RIM> 実区間行列の型
   * @param <CIS> 複素区間スカラーの型
   * @param <CIM> 複素区間行列の型
   * @param <RS> 実スカラーの型
   * @param <RM> 実行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @param a 区間。
   * @param b 区間。
   * @return 含まれていればtrue, そうでなければfalse。
   */
  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> boolean in0(
      final CIS a, final CIS b) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    manager.setRoundMode(RoundMode.ROUND_DOWN);
    RS r1 = (b.getMiddle().getRealPart()).subtract((a.getMiddle().getRealPart()));
    RS l1 = (b.getMiddle().getImaginaryPart()).subtract((a.getMiddle().getImaginaryPart()));

    manager.setRoundMode(RoundMode.ROUND_UP);
    RS r2 = (b.getMiddle().getRealPart()).subtract((a.getMiddle().getRealPart()));
    RS l2 = (b.getMiddle().getImaginaryPart()).subtract((a.getMiddle().getImaginaryPart()));

    RS r = r1.abs().max(r2.abs());
    RS l = l1.abs().max(l2.abs());
    CS c = r1.toComplex().create(r, l);
    RS d = c.abs().getRealPart();

    boolean res;
    if (a.getRadius().isZero()) {
      res = b.getRadius().isGreaterThan(d);
    } else {
      res = b.getRadius().isGreaterThan((d.add(a.getRadius())));
    }

    manager.setRoundMode(oldRoundMode);
    return res;
  }

  /**
   * 点aが区間b内に含まれるかどうかを判定します。
   * 
   * @param a 点
   * @param b 区間
   * @return 含まれていればtrue, そうでなければfalse
   */
  public static boolean in0(final double a, final DoubleIntervalNumber b) {
    DoubleIntervalNumber interval_a = new DoubleIntervalNumber(a);
    return in0(interval_a, b);
  }

  /**
   * 点aが区間b内に含まれるかどうかを判定します。
   * 
   * @param a 点
   * @param b 区間
   * @return 含まれていればtrue, そうでなければfalse
   */
  public static boolean in0(final DoubleComplexNumber a, final DoubleComplexIntervalNumber b) {
    DoubleComplexIntervalNumber interval_a = new DoubleComplexIntervalNumber(a);
    return in0(interval_a, b);
  }

  /**
   * 区間行列Aが区間行列B内に含まれるかどうかを判定します。
   * 
   * @param A 区間
   * @param B 区間
   * @param <RIS> 実区間スカラーの型
   * @param <RIM> 実区間行列の型
   * @param <CIS> 複素区間スカラーの型
   * @param <CIM> 複素区間行列の型
   * @param <RS> 実スカラーの型
   * @param <RM> 実行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @return 区間に含まれる要素はtrue, 含まれない要素はfalse
   */
  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BooleanMatrix in0(
      final RIM A, final RIM B) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    BooleanMatrix res = B.getInfimum().compareElementWise(".<", A.getInfimum()).andElementWise(A.getSupremum().compareElementWise(".<", B.getSupremum())); //$NON-NLS-1$ //$NON-NLS-2$
    manager.setRoundMode(oldRoundMode);
    return res;
  }

  /**
   * 区間行列Aが区間行列B内に含まれるかどうかを判定します。
   * 
   * @param A 区間
   * @param B 区間
   * @param <RIS> 実区間スカラーの型
   * @param <RIM> 実区間行列の型
   * @param <CIS> 複素区間スカラーの型
   * @param <CIM> 複素区間行列の型
   * @param <RS> 実スカラーの型
   * @param <RM> 実行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @return 区間に含まれる要素はtrue, 含まれない要素はfalse
   */
  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BooleanMatrix in0(
      final CIM A, final CIM B) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    manager.setRoundMode(RoundMode.ROUND_DOWN);
    RM R1 = B.getMiddle().getRealPart().subtract(A.getMiddle().getRealPart());
    RM L1 = B.getMiddle().getImaginaryPart().subtract(A.getMiddle().getImaginaryPart());

    manager.setRoundMode(RoundMode.ROUND_UP);
    RM R2 = B.getMiddle().getRealPart().subtract(A.getMiddle().getRealPart());
    RM L2 = B.getMiddle().getImaginaryPart().subtract(A.getMiddle().getImaginaryPart());

    RM R = (R1).absElementWise().maxElementWise((R2).absElementWise());
    RM L = (L1).absElementWise().maxElementWise((L2).absElementWise());
    CM C = R.createComplex(R, L);
    RM d = (C).absElementWise().getRealPart();

    BooleanMatrix res;
    if (A.getRadius().isZero()) {
      res = d.compareElementWise(".<", B.getRadius()); //$NON-NLS-1$
    } else {
      res = d.add(A.getRadius()).compareElementWise(".<", B.getRadius()); //$NON-NLS-1$
    }

    manager.setRoundMode(oldRoundMode);
    return res;
  }

  /**
   * @param A A
   * @param B B
   * @param <RIS> 実区間スカラーの型
   * @param <RIM> 実区間行列の型
   * @param <CIS> 複素区間スカラーの型
   * @param <CIM> 複素区間行列の型
   * @param <RS> 実スカラーの型
   * @param <RM> 実行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @return boolean
   */
  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BooleanMatrix in0forRump(
      final RIM A, final RIM B) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    BooleanMatrix res = B.getInfimum().compareElementWise(".<=", A.getInfimum()).andElementWise(A.getSupremum().compareElementWise(".<=", B.getSupremum())); //$NON-NLS-1$ //$NON-NLS-2$
    manager.setRoundMode(oldRoundMode);
    return res;
  }

//  /**
//   * @param A A
//   * @param B B
//   * @param <IS> 実区間スカラーの型
//   * @param <IM> 実区間行列の型
//   *  @param <S> 実スカラーの型
//   * @param <M> 実行列の型
//   *  @return boolean
//   */
//  public static <IS extends IntervalNumericalScalar<IS, IM, S, M>, IM extends IntervalNumericalMatrix<IS, IM, S, M>,  S extends NumericalScalar<S, M>, M extends NumericalMatrix<S, M>> BooleanMatrix in0forRump(
//      final IM A, final IM B) {
//    RoundModeManager manager = RoundModeManager.getManager();
//    RoundMode oldRoundMode = manager.getRoundMode();
//
//    BooleanMatrix res = B.getInfimum().compareElementWise(".<=", A.getInfimum()).andElementWise(A.getSupremum().compareElementWise(".<=", B.getSupremum())); //$NON-NLS-1$ //$NON-NLS-2$
//    manager.setRoundMode(oldRoundMode);
//    return res;
//  }
  
  /**
   * @param A A
   * @param B B
   * @param <RIS> 実区間スカラーの型
   * @param <RIM> 実区間行列の型
   * @param <CIS> 複素区間スカラーの型
   * @param <CIM> 複素区間行列の型
   * @param <RS> 実スカラーの型
   * @param <RM> 実行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @return boolean
   */
  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BooleanMatrix in0forRump(
      final CIM A, final CIM B) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    BooleanMatrix res;

    manager.setRoundMode(RoundMode.ROUND_DOWN);
    RM R1 = B.getMiddle().getRealPart().subtract(A.getMiddle().getRealPart());
    RM L1 = B.getMiddle().getImaginaryPart().subtract(A.getMiddle().getImaginaryPart());

    manager.setRoundMode(RoundMode.ROUND_UP);
    RM R2 = B.getMiddle().getRealPart().subtract(A.getMiddle().getRealPart());
    RM L2 = B.getMiddle().getImaginaryPart().subtract(A.getMiddle().getImaginaryPart());

    RM R = R1.absElementWise().maxElementWise(R2.absElementWise());
    RM L = L1.absElementWise().maxElementWise(L2.absElementWise());
    CM C = R.createComplex(R, L);
    RM d = C.absElementWise().getRealPart();

    if (A.getRadius().isZero()) {
      res = d.compareElementWise(".<=", B.getRadius()); //$NON-NLS-1$
    } else {
      res = d.add(A.getRadius()).compareElementWise(".<=", B.getRadius()); //$NON-NLS-1$
    }

    manager.setRoundMode(oldRoundMode);
    return res;
  }

  /**
   * 行列Aが区間行列B内に含まれるかどうかを判定します。
   * 
   * @param A 行列
   * @param B 区間行列
   * @param <RIS> 実区間スカラーの型
   * @param <RIM> 実区間行列の型
   * @param <CIS> 複素区間スカラーの型
   * @param <CIM> 複素区間行列の型
   * @param <RS> 実スカラーの型
   * @param <RM> 実行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @return 区間にに含まれる要素はtrue, 含まれない要素はfalse
   */
  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BooleanMatrix in0(
      final RM A, final RIM B) {
    RIM IA = B.create(A);
    return in0forRump(IA, B);
  }

  /**
   * 行列Aが区間行列B内に含まれるかどうかを判定します。
   * 
   * @param A 行列
   * @param B 区間行列
   * @param <RIS> 実区間スカラーの型
   * @param <RIM> 実区間行列の型
   * @param <CIS> 複素区間スカラーの型
   * @param <CIM> 複素区間行列の型
   * @param <RS> 実スカラーの型
   * @param <RM> 実行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @return 区間にに含まれる要素はtrue, 含まれない要素はfalse
   */
  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> BooleanMatrix in0(
      final CM A, final CIM B) {
    CIM IA = B.create(A);
    return in0forRump(IA, B);
  }

  /**
   * <code>BooleanMatrix</code> の要素のtrueを1,falseを0とした行列を返します。
   * 
   * @param Bool 二値行列
   * @return trueを1にfalseを0に変換した行列
   */
  public static IntMatrix toZOMatrix(final BooleanMatrix Bool) {
    IntMatrix zomatrix = new IntMatrix(Bool.getRowSize(), Bool.getColumnSize());
    for (int i = 0; i < Bool.getRowSize(); i++) {
      for (int j = 0; j < Bool.getColumnSize(); j++) {
        if (Bool.getElement(i + 1, j + 1)) {
          zomatrix.setElement(i + 1, j + 1, 1);
        }
      }
    }

    return zomatrix;
  }

  /**
   * Comparison Matrixを求めます。
   * 
   * @param <IS> 区間スカラーの型
   * @param <IM> 区間行列の型
   * @param <S> 成分の型
   * @param <M> 行列の型
   * @param im 区間行列
   * @return ComparisonMatrix
   */
  public static <IS extends IntervalNumericalScalar<IS, IM, S, M>, IM extends IntervalNumericalMatrix<IS, IM, S, M>, S extends NumericalScalar<S, M>, M extends NumericalMatrix<S, M>> M comparisonMatrix(
      final IM im) {
    M Ac = im.abssElementWise().unaryMinus();
    M D = Ac.diagonalToVector();
    M M = (im.diagonalToVector()).migElementWise();

    S unit = Ac.getElement(1, 1);
    S[] d1 = unit.createArray(D.length());
    S[] d2 = unit.createArray(M.length());

    for (int i = 0; i < D.length(); i++) {
      d1[i] = D.getElement(i + 1, 1);
      d2[i] = M.getElement(i + 1, 1);
    }

    Ac = Ac.subtract(d1[0].createGrid(d1).vectorToDiagonal()).add(d2[0].createGrid(d2).vectorToDiagonal());
    return Ac;
  }

  /**
   * Comparison Matrixを求めます。
   * 
   * @param <RIS> 実区間スカラーの型
   * @param <RIM> 実区間行列の型
   * @param <CIS> 複素区間スカラーの型
   * @param <CIM> 複素区間行列の型
   * @param <RS> 実スカラーの型
   * @param <RM> 実行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @param im 区間行列
   * @return ComparisonMatrix
   */
  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RM comparisonMatrix(
      final RIM im) {
    RM Ac = im.abssElementWise().unaryMinus();
    RM D = Ac.diagonalToVector();
    RM M = (im.diagonalToVector()).migElementWise();

    RS unit = Ac.getElement(1, 1);
    RS[] d1 = unit.createArray(D.length());
    RS[] d2 = unit.createArray(M.length());

    for (int i = 0; i < D.length(); i++) {
      d1[i] = D.getElement(i + 1, 1);
      d2[i] = M.getElement(i + 1, 1);
    }

    Ac = Ac.subtract(d1[0].createGrid(d1).vectorToDiagonal()).add(d2[0].createGrid(d2).vectorToDiagonal());
    return Ac;
  }

  /**
   * Comparison Matrixを求めます。
   * 
   * @param <RIS> 実区間スカラーの型
   * @param <RIM> 実区間行列の型
   * @param <CIS> 複素区間スカラーの型
   * @param <CIM> 複素区間行列の型
   * @param <RS> 実スカラーの型
   * @param <RM> 実行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @param im 区間行列
   * @return ComparisonMatrix
   */
  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> CM comparisonMatrix(
      final CIM im) {
    CM Ac = im.abssElementWise().unaryMinus();
    CM D = Ac.diagonalToVector();
    CM M = (im.diagonalToVector()).migElementWise();

    CS unit = Ac.getElement(1, 1);
    CS[] d1 = unit.createArray(D.length());
    CS[] d2 = unit.createArray(M.length());

    for (int i = 0; i < D.length(); i++) {
      d1[i] = D.getElement(i + 1, 1);
      d2[i] = M.getElement(i + 1, 1);
    }

    Ac = Ac.subtract(d1[0].createGrid(d1).vectorToDiagonal()).add(d2[0].createGrid(d2).vectorToDiagonal());
    return Ac;
  }

  //  /**
  //   * 中心と半径による区間行列を作成します。
  //   * 
  //   * @param <RIS> 実区間スカラーの型
  //   * @param <RIM> 実区間行列の型
  //   * @param <CIS> 複素区間スカラーの型
  //   * @param <CIM> 複素区間行列の型
  //   * @param <RS> 実スカラーの型
  //   * @param <RM> 実行列の型
  //   * @param <CS> 複素スカラーの型
  //   * @param <CM> 複素行列の型
  //   * @param mid 中心
  //   * @param rad 半径
  //   * @return 区間行列
  //   */
  //  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RIM midrad(final RM mid, final RM rad) {
  //    RoundModeManager manager = RoundModeManager.getManager();
  //    RoundMode oldRoundMode = manager.getRoundMode();
  //
  //    manager.setRoundMode(RoundMode.ROUND_DOWN);
  //    RM inf = mid.subtract(rad);
  //
  //    manager.setRoundMode(RoundMode.ROUND_UP);
  //    RM sup = mid.add(rad);
  //
  //    RIM midrad = mid.createI new AbstractIntervalRealNumericalMatrix(inf, sup);
  //
  //    manager.setRoundMode(oldRoundMode);
  //    return midrad;
  //
  //    manager.setRoundMode(oldRoundMode);
  //    throw new RuntimeException("Not implemented."); //$NON-NLS-1$
  //  }

  //  /**
  //   * 中心と半径による区間行列を作成します。
  //   * 
  //   * @param mid 中心
  //   * @param rad 半径
  //   * @return 区間行列
  //   */
  //  public static IntervalMatrix<?> midrad(final NumericalMatrix<?> mid, final NumericalMatrix<?> rad) {
  //    RoundModeManager manager = RoundModeManager.getManager();
  //    RoundMode oldRoundMode = manager.getRoundMode();
  //
  //    if (mid.isReal() && rad.isReal()) {
  //      manager.setRoundMode(RoundMode.ROUND_DOWN);
  //      NumericalMatrix<?> inf = mid.subtract(rad);
  //
  //      manager.setRoundMode(RoundMode.ROUND_UP);
  //      NumericalMatrix<?> sup = mid.add(rad);
  //
  //      IntervalMatrix<?> midrad = new AbstractIntervalRealNumericalMatrix(inf, sup);
  //
  //      manager.setRoundMode(oldRoundMode);
  //      return midrad;
  //    } else if (mid.isComplex() && rad.isReal()) {
  //      IntervalMatrix<?> midrad = new AbstractIntervalComplexNumericalMatrix(mid, rad);
  //      manager.setRoundMode(oldRoundMode);
  //      return midrad;
  //    }
  //
  //    manager.setRoundMode(oldRoundMode);
  //    throw new RuntimeException("Not implemented."); //$NON-NLS-1$
  //  }

  //  /**
  //   * 中心と半径により区間を作成します。
  //   * 
  //   * @param mid 中心。
  //   * @param rad 半径。
  //   * @return 区間。
  //   */
  //  public static IntervalScalar<?> midrad(final double mid, final double rad) {
  //    RoundModeManager manager = RoundModeManager.getManager();
  //    RoundMode oldRoundMode = manager.getRoundMode();
  //
  //    manager.setRoundMode(RoundMode.ROUND_DOWN);
  //    double inf = mid - rad;
  //
  //    manager.setRoundMode(RoundMode.ROUND_UP);
  //    double sup = mid + rad;
  //
  //    DoubleIntervalNumber midrad = new DoubleIntervalNumber(inf, sup);
  //
  //    manager.setRoundMode(oldRoundMode);
  //    return midrad;
  //  }

  /**
   * 中心と半径により区間を作成します。
   * 
   * @param mid 中心
   * @param rad 半径
   * @return 区間
   */
  public static DoubleComplexIntervalNumber midrad(final DoubleComplexNumber mid, final double rad) {
    return new DoubleComplexIntervalNumber(mid, rad);
  }

  /**
   * 自身のノルムを求めます。
   * 
   * @param <IS> 区間スカラーの型
   * @param <IM> 区間行列の型
   * @param <S> スカラーの型
   * @param <M> 行列の型
   * 
   * @param A 区間行列
   * @param type ノルムの種類
   * @return norm(NormType.FROBENIUS)で自身のフロベニウスノルムを、norm(NormType.INFINITY)で無限大ノルム 、 norm(NormType.ONE)で1-ノルムを、norm(NormType.TWO)で最大特異値。
   */
  public static <IS extends IntervalNumericalScalar<IS, IM, S, M>, IM extends IntervalNumericalMatrix<IS, IM, S, M>, S extends NumericalScalar<S, M>, M extends NumericalMatrix<S, M>> S norm(
      final IM A, final NormType type) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    M x = A.abssElementWise();

    manager.setRoundMode(RoundMode.ROUND_UP);
    S res = x.norm(type);

    manager.setRoundMode(oldRoundMode);
    return res;
  }

  //  /**
  //   * AとBの共通の要素を返します。共通でない要素のところにはNaNを入れて返します。
  //   * 
  //   * @param A 行列
  //   * @param B 行列
  //   * @return 共通要素を持った区間行列
  //   */
  //  public static IntervalMatrix<?> intersect(final NumericalMatrix<?> A, final NumericalMatrix<?> B) {
  //    //    if (A instanceof IntervalMatrix && B instanceof IntervalMatrix) {
  //    //      return intersect((IntervalMatrix<?>)A, (IntervalMatrix<?>)B);
  //    //    }
  //
  //    IntervalMatrix<?> intervalA = IntervalMatrixFactory.toInterval(A);
  //    IntervalMatrix<?> intervalB = IntervalMatrixFactory.toInterval(B);
  //    return intersect(intervalA, intervalB);
  //  }

  /**
   * AとBの共通の要素を返します。共通でない要素のところにはNaNを入れて返します。
   * 
   * @param <RIS> 実区間スカラーの型
   * @param <RIM> 実区間行列の型
   * @param <CIS> 複素区間スカラーの型
   * @param <CIM> 複素区間行列の型
   * @param <RS> 実スカラーの型
   * @param <RM> 実行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @param A 区間行列
   * @param B 区間行列
   * @return 共通要素を持った区間行列
   */
  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> RIM intersect(
      final RIM A, final RIM B) {

    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    RM inf = (A.getInfimum()).maxElementWise(B.getInfimum());
    RM sup = (A.getSupremum()).minElementWise(B.getSupremum());

    BooleanMatrix index = inf.compareElementWise(".>", sup).orElementWise(A.getInfimum().isNanElementWise()).orElementWise(A.getSupremum().isNanElementWise()).orElementWise( //$NON-NLS-1$
        B.getInfimum().isNanElementWise()).orElementWise(B.getSupremum().isNanElementWise());
    if (index.anyTrue()) {
      RS unit = A.getMiddle().getElement(1);

      for (int i = 0; i < index.getRowSize(); i++) {
        for (int j = 0; j < index.getColumnSize(); j++) {
          if (index.getElement(i + 1, j + 1)) {
            (inf).setElement(i + 1, j + 1, unit.getNaN());
            (sup).setElement(i + 1, j + 1, unit.getNaN());
          }
        }
      }
    }

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

//  /**
//   * AとBの共通の要素を返します。共通でない要素のところにはNaNを入れて返します。
//   * 
//   * @param <IS> 実区間スカラーの型
//   * @param <IM> 実区間行列の型
//   * @param <S> 実スカラーの型
//   * @param <M> 実行列の型
//   * @param A 区間行列
//   * @param B 区間行列
//   * @return 共通要素を持った区間行列
//   */
//  public static <IS extends IntervalNumericalScalar<IS, IM, S, M>, IM extends IntervalNumericalMatrix<IS, IM, S, M>, S extends NumericalScalar<S, M>, M extends NumericalMatrix<S, M>> IM intersect(
//      final IM A, final IM B) {
//    RoundModeManager manager = RoundModeManager.getManager();
//    RoundMode oldRoundMode = manager.getRoundMode();
//
//    M inf = (A.getInfimum()).maxElementWise(B.getInfimum());
//    M sup = (A.getSupremum()).minElementWise(B.getSupremum());
//
//    BooleanMatrix index = inf.compareElementWise(".>", sup).orElementWise(A.getInfimum().isNanElementWise()).orElementWise(A.getSupremum().isNanElementWise()).orElementWise( //$NON-NLS-1$
//        B.getInfimum().isNanElementWise()).orElementWise(B.getSupremum().isNanElementWise());
//    if (index.anyTrue()) {
//      S unit = A.getMiddle().getElement(1);
//
//      for (int i = 0; i < index.getRowSize(); i++) {
//        for (int j = 0; j < index.getColumnSize(); j++) {
//          if (index.getElement(i + 1, j + 1)) {
//            (inf).setElement(i + 1, j + 1, unit.getNaN());
//            (sup).setElement(i + 1, j + 1, unit.getNaN());
//          }
//        }
//      }
//    }
//
//    manager.setRoundMode(oldRoundMode);
//    return A.createInfSup(inf, sup);
//  }
  
  
  /**
   * AとBの共通の要素を返します。共通でない要素のところにはNaNを入れて返します。
   * 
   * @param <RIS> 実区間スカラーの型
   * @param <RIM> 実区間行列の型
   * @param <CIS> 複素区間スカラーの型
   * @param <CIM> 複素区間行列の型
   * @param <RS> 実スカラーの型
   * @param <RM> 実行列の型
   * @param <CS> 複素スカラーの型
   * @param <CM> 複素行列の型
   * @param A 区間行列
   * @param B 区間行列
   * @return 共通要素を持った区間行列
   */
  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> CIM intersect(
      final CIM A, final CIM B) {

    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    // bounds for distance of midpoints
    CM D = A.getMiddle().subtract(B.getMiddle());
    CM d = D.create(D.getRealPart().absElementWise(), D.getImaginaryPart().absElementWise());

    CS eps = A.getMiddle().getElement(1).getMachineEpsilon();

    manager.setRoundMode(RoundMode.ROUND_DOWN);
    RM distinf = (d.multiply(eps.unaryMinus().add(1))).absElementWise().getRealPart();

    manager.setRoundMode(RoundMode.ROUND_UP);
    RM distsup = (d.multiply(eps.add(1))).absElementWise().getRealPart();

    RIM dist = A.getRealPart().createInfSup(distinf, distsup);
    RIM arad = A.getRealPart().create(A.getRadius());
    RIM brad = B.getRealPart().create(B.getRadius());

    RIM x_d = ((arad.add(brad)).multiplyElementWise((arad.subtract(brad))).divideElementWise((dist.multiplyElementWise(dist))).addElementWise(1)).divide(2);// TODO
    // ほんとは/2をやりたい

    CIM c_mid = A.create(A.getMiddle()).add(x_d.toComplex().multiplyElementWise((B.create(B.getMiddle().subtract(A.getMiddle())))));

    CM cmid = c_mid.getMiddle();
    RM crad = c_mid.getRadius();

    manager.setRoundMode(RoundMode.ROUND_UP);
    RIM x = dist.multiplyElementWise(x_d);
    RIM ccrad = arad.multiplyElementWise(arad).add((x.unaryMinus()).multiplyElementWise(x));

    manager.setRoundMode(RoundMode.ROUND_UP);
    RM cccrad = (ccrad.getSupremum()).sqrtElementWise();
    crad = crad.add(cccrad);

    // entire interval A in B for some indices
    manager.setRoundMode(RoundMode.ROUND_UP);
    BooleanMatrix index = distsup.add(A.getRadius()).compareElementWise(".<=", B.getRadius()); //$NON-NLS-1$

    if (index.anyTrue()) {
      for (int i = 0; i < index.getRowSize(); i++) {
        for (int j = 0; j < index.getColumnSize(); j++) {
          if (index.getElement(i + 1, j + 1)) {
            // TODO 具体的な行列を使わない
            cmid.setElement(i + 1, j + 1, (A.getMiddle()).getElement(i, j));
            crad.setElement(i + 1, j + 1, (A.getRadius()).getElement(i, j));
          }
        }
      }
    }

    // entire interval B in A for some indices
    index = distsup.add(B.getRadius()).compareElementWise(".<=", A.getRadius()); //$NON-NLS-1$

    if (index.anyTrue()) {
      for (int i = 0; i < index.getRowSize(); i++) {
        for (int j = 0; j < index.getColumnSize(); j++) {
          if (index.getElement(i + 1, j + 1)) {
            // TODO 具体的な行列を使わない
            cmid.setElement(i + 1, j + 1, (B.getMiddle()).getElement(i + 1, j + 1));
            crad.setElement(i + 1, j + 1, (B.getRadius()).getElement(i + 1, j + 1));
          }
        }
      }
    }

    index = crad.compareElementWise("!=", crad.createZero()).orElementWise(A.getMiddle().isNanElementWise()) //$NON-NLS-1$
        .orElementWise(A.getRadius().isNanElementWise()).orElementWise(B.getMiddle().isNanElementWise()).orElementWise(B.getRadius().isNanElementWise());

    if (index.anyTrue()) {
      RS runit = A.getMiddle().getElement(1).getRealPart();
      CS cunit = A.getMiddle().getElement(1);

      for (int i = 0; i < index.getRowSize(); i++) {
        for (int j = 0; j < index.getColumnSize(); j++) {
          if (index.getElement(i + 1, j + 1)) {
            cmid.setElement(i + 1, j + 1, cunit.getNaN());
            crad.setElement(i + 1, j + 1, runit.getNaN());
          }
        }
      }
    }

    manager.setRoundMode(oldRoundMode);
    return A.createMidRad(cmid, crad);
  }

  //  /**
  //   * 区間行列Aの負を返します。
  //   * 
  //   * @param <E> 成分の型
  //   * 
  //   * @param A 区間行列
  //   * @return 負の値
  //   */
  //  public static <E extends NumericalScalar<E>> IntervalMatrix<E> unaryMinus(final IntervalMatrix<E> A) {
  //    return A.multiply(-1);
  //  }

  //  /**
  //   * <code>ComplexIntervalMatrix</code> に変換します。
  //   * 
  //   * @param <RIS> 実区間スカラーの型
  //   * @param <RIM> 実区間行列の型
  //   * @param <CIS> 複素区間スカラーの型
  //   * @param <CIM> 複素区間行列の型
  //   * @param <RS> 実スカラーの型
  //   * @param <RM> 実行列の型
  //   * @param <CS> 複素スカラーの型
  //   * @param <CM> 複素行列の型
  //   * @param im 区間行列
  //   * @return 区間行列
  //   */
  //  public static <RIS extends IntervalRealNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RIM extends IntervalRealNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIS extends IntervalComplexNumericalScalar<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, CIM extends IntervalComplexNumericalMatrix<RIS, RIM, CIS, CIM, RS, RM, CS, CM>, RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> CIM toComplex(
  //      final RIM im) {
  //
  //    return im.toComplex();
  //}

  //  /**
  //   * 実数区間を複素数区間に変換します。
  //   * 
  //   * @param i 実数区間
  //   * @return 複素数区間
  //   */
  //  public static IntervalScalar<?> toComplex(final IntervalScalar<?> i) {
  //    RoundModeManager manager = RoundModeManager.getManager();
  //    RoundMode oldRoundMode = manager.getRoundMode();
  //
  //    DoubleComplexNumber mid;
  //    double rad;
  //
  //    if (i instanceof DoubleIntervalNumber) {
  //      if (i.getInfimum().equals(i.getSupremum())) {
  //        mid = new DoubleComplexNumber((i.getInfimum()).doubleValue(), 0.0);
  //        rad = 0;
  //      } else {
  //        manager.setRoundMode(RoundMode.ROUND_UP);
  //        
  //        mid = (DoubleComplexNumber)i.getInfimum().add(i.getSupremum().subtract(i.getInfimum()).multiply(0.5));
  //        rad = (i.getMiddle().subtract(i.getInfimum())).doubleValue();
  //      }
  //    } else {
  //      manager.setRoundMode(oldRoundMode);
  //      return i;
  //    }
  //
  //    manager.setRoundMode(oldRoundMode);
  //    return new DoubleComplexIntervalNumber(mid, rad);
  //  }
}