class SwaptionBermudaFixedIborHullWhiteNumericalIntegrationMethod {

  public MultipleCurrencyAmount presentValue(final SwaptionBermudaFixedIbor swaption, final HullWhiteOneFactorProviderInterface hullWhite) {
    ArgumentChecker.notNull(swaption, "Swaption");
    ArgumentChecker.notNull(hullWhite, "Hull-White provider");
    MulticurveProviderInterface multicurves = hullWhite.getMulticurveProvider();
    HullWhiteOneFactorPiecewiseConstantParameters parameters = hullWhite.getHullWhiteParameters();
    Currency ccy = swaption.getCurrency();
    final int nbExpiry = swaption.getExpiryTime().length;
    Validate.isTrue(nbExpiry > 1, "At least two expiry dates required for this method");

    double tmpdb;

    final double[] theta = new double[nbExpiry + 1]; // Extended expiry time (with 0).
    theta[0] = 0.0;
    System.arraycopy(swaption.getExpiryTime(), 0, theta, 1, nbExpiry);
    final AnnuityPaymentFixed[] cashflow = new AnnuityPaymentFixed[nbExpiry];
    for (int loopexp = 0; loopexp < nbExpiry; loopexp++) {
      cashflow[loopexp] = swaption.getUnderlyingSwap()[loopexp].accept(CFEC, multicurves);
    }
    final int[] n = new int[nbExpiry];
    final double[][][] alpha = new double[nbExpiry][][];
    final double[][][] alpha2 = new double[nbExpiry][][]; // alpha^2

    for (int loopexp = 0; loopexp < nbExpiry; loopexp++) {
      n[loopexp] = cashflow[loopexp].getNumberOfPayments();
      alpha[loopexp] = new double[loopexp + 1][];
      alpha2[loopexp] = new double[loopexp + 1][];
      for (int k = 0; k <= loopexp; k++) {
        alpha[loopexp][k] = new double[n[loopexp]];
        alpha2[loopexp][k] = new double[n[loopexp]];
        for (int l = 0; l < alpha[loopexp][k].length; l++) {
          alpha[loopexp][k][l] = MODEL.alpha(parameters, theta[k], theta[k + 1], theta[k + 1], cashflow[loopexp].getNthPayment(l).getPaymentTime());
          alpha2[loopexp][k][l] = alpha[loopexp][k][l] * alpha[loopexp][k][l];
        }
      }
    }

    final int nbPoint2 = 2 * NB_POINT + 1;
    final int[] startInt = new int[nbExpiry - 1];
    final int[] endInt = new int[nbExpiry - 1];
    for (int i = 1; i < nbExpiry - 1; i++) {
      startInt[i] = 0;
      endInt[i] = nbPoint2 - 1;
    }
    startInt[0] = NB_POINT;
    endInt[0] = NB_POINT;

    final double[][] t = new double[nbExpiry][]; // payment time
    final double[][] dfS = new double[nbExpiry][]; // discount factor
    final double[] beta = new double[nbExpiry];
    final double[][] h = new double[nbExpiry][];
    final double[][] sa2 = new double[nbExpiry][];

    for (int loopexp = 0; loopexp < nbExpiry; loopexp++) {
      beta[loopexp] = MODEL.beta(parameters, theta[loopexp], theta[loopexp + 1]);
      t[loopexp] = new double[n[loopexp]];
      dfS[loopexp] = new double[n[loopexp]];
      h[loopexp] = new double[n[loopexp]];
      sa2[loopexp] = new double[n[loopexp]];
      for (int loopcf = 0; loopcf < n[loopexp]; loopcf++) {
        t[loopexp][loopcf] = cashflow[loopexp].getNthPayment(loopcf).getPaymentTime();
        dfS[loopexp][loopcf] = multicurves.getDiscountFactor(ccy, t[loopexp][loopcf]);
        h[loopexp][loopcf] = (1 - Math.exp(-parameters.getMeanReversion() * t[loopexp][loopcf])) / parameters.getMeanReversion();
        tmpdb = 0.0;
        for (int k = 0; k <= loopexp; k++) {
          tmpdb += alpha2[loopexp][k][loopcf];
        }
        sa2[loopexp][loopcf] = tmpdb;
      }
    }
    final double[] discountedCashFlowN = new double[n[nbExpiry - 1]];
    for (int loopcf = 0; loopcf < n[nbExpiry - 1]; loopcf++) {
      discountedCashFlowN[loopcf] = dfS[nbExpiry - 1][loopcf] * cashflow[nbExpiry - 1].getNthPayment(loopcf).getAmount();
    }
    final double lambda = MODEL.lambda(discountedCashFlowN, sa2[nbExpiry - 1], h[nbExpiry - 1]);
    final double[] betaSort = new double[nbExpiry];
    System.arraycopy(beta, 0, betaSort, 0, nbExpiry);
    Arrays.sort(betaSort);
    final double minbeta = betaSort[0];
    final double maxbeta = betaSort[nbExpiry - 1];

    final double b = Math.min(10 * minbeta, maxbeta);
    final double epsilon = -2.0 / NB_POINT * NORMAL.getInverseCDF(1.0 / (200.0 * NB_POINT)) * b; // <-
    final double[] bX = new double[nbPoint2];
    for (int looppt = 0; looppt < nbPoint2; looppt++) {
      bX[looppt] = -NB_POINT * epsilon + looppt * epsilon;
    }
    final double[] bX2 = new double[4 * NB_POINT + 1];
    for (int looppt = 0; looppt < 4 * NB_POINT + 1; looppt++) {
      bX2[looppt] = -2 * NB_POINT * epsilon + looppt * epsilon;
    }
    final double[] htheta = new double[nbExpiry];
    for (int loopexp = 0; loopexp < nbExpiry; loopexp++) {
      htheta[loopexp] = (1 - Math.exp(-parameters.getMeanReversion() * theta[loopexp + 1])) / parameters.getMeanReversion();
    }

    final double[][] vZ = new double[nbExpiry - 1][nbPoint2];
    for (int i = nbExpiry - 2; i >= 0; i--) {
      for (int looppt = 0; looppt < nbPoint2; looppt++) {
        vZ[i][looppt] = Math.exp(bX[looppt] * htheta[i]);
      }
    }

    final double[][] vW = new double[nbExpiry][]; // Swaption
    final double[][] vT = new double[nbExpiry - 1][]; // Swap

    final double omega = -Math.signum(cashflow[nbExpiry - 1].getNthPayment(0).getAmount());
    final double[] kappaL = new double[nbPoint2];
    for (int looppt = 0; looppt < nbPoint2; looppt++) {
      kappaL[looppt] = (lambda - bX[looppt]) / beta[nbExpiry - 1];
    }

    final double[] sa2N1 = new double[n[nbExpiry - 1]];
    for (int i = 0; i < n[nbExpiry - 1]; i++) {
      tmpdb = 0;
      for (int k = 0; k <= nbExpiry - 2; k++) {
        tmpdb += alpha2[nbExpiry - 1][k][i];
      }
      sa2N1[i] = tmpdb;
    }

    vW[nbExpiry - 1] = new double[4 * NB_POINT + 1];
    for (int j = 0; j < n[nbExpiry - 1]; j++) {
      for (int looppt = 0; looppt < nbPoint2; looppt++) {
        vW[nbExpiry - 1][NB_POINT + looppt] += discountedCashFlowN[j] * Math.exp(-sa2N1[j] / 2.0 - h[nbExpiry - 1][j] * bX[looppt])
            * NORMAL.getCDF(omega * (kappaL[looppt] + alpha[nbExpiry - 1][nbExpiry - 1][j]));
      }
    }
    for (int looppt = 0; looppt < NB_POINT; looppt++) {
      vW[nbExpiry - 1][looppt] = vW[nbExpiry - 1][NB_POINT];
    }
    for (int looppt = 0; looppt < NB_POINT; looppt++) {
      vW[nbExpiry - 1][3 * NB_POINT + 1 + looppt] = vW[nbExpiry - 1][3 * NB_POINT];
    }

    final double c1sqrt2pi = 1.0 / Math.sqrt(2 * Math.PI);
    final double[][] pvcfT = new double[nbExpiry - 1][];
    double[] vL; // Left side of intersection
    double[] vR; // Right side of intersection
    double[][] labc;
    double[][] rabc;
    final double[][] labcM = new double[3][4 * NB_POINT + 1];
    final double[][] rabcM = new double[3][4 * NB_POINT + 1];

    final double[] dabc = new double[3];
    final int[] indSwap = new int[nbExpiry - 1]; // index of the intersection
    double xroot;
    final double[][] xN = new double[nbExpiry - 1][nbPoint2];
    double ci;
    double coi;
    int is;
    final double[] ncdf0 = new double[nbPoint2];
    final double[] ncdf1 = new double[nbPoint2];
    final double[] ncdf2 = new double[nbPoint2];
    final double[] ncdf0X = new double[nbPoint2 + 1];
    final double[] ncdf1X = new double[nbPoint2 + 1];
    final double[] ncdf2X = new double[nbPoint2 + 1];
    double ncdf0x;
    double ncdf1x;
    double ncdf2x;
    double ncdfinit;

    // Main loop for the different expiry dates (except the last one)
    for (int i = nbExpiry - 2; i >= 0; i--) {
      vW[i] = new double[4 * NB_POINT + 1];
      vT[i] = new double[4 * NB_POINT + 1];
      // T: swap
      pvcfT[i] = new double[n[i]];
      for (int j = 0; j < n[i]; j++) {
        pvcfT[i][j] = cashflow[i].getNthPayment(j).getAmount() * dfS[i][j];
        for (int looppt = 0; looppt < 4 * NB_POINT + 1; looppt++) {
          vT[i][looppt] += pvcfT[i][j] * Math.exp(-sa2[i][j] / 2.0 - h[i][j] * bX2[looppt]);
        }
      }
      // Preparation
      for (int looppt = 0; looppt < nbPoint2; looppt++) {
        xN[i][looppt] = bX[looppt] / beta[i];
      }
      ci = htheta[i] * beta[i];
      coi = Math.exp(ci * ci / 2);

      // Left/Right
      if (omega < 0) {
        vL = vW[i + 1];
        vR = vT[i];
      } else {
        vR = vW[i + 1];
        vL = vT[i];
      }
      indSwap[i] = 0;
      while (vL[indSwap[i] + 1] >= vR[indSwap[i] + 1]) {
        indSwap[i]++;
      }
      // Parabola fit
      labc = parafit(epsilon / beta[i], vL);
      rabc = parafit(epsilon / beta[i], vR);
      for (int k = 0; k < 3; k++) {
        dabc[k] = labc[k][indSwap[i]] - rabc[k][indSwap[i]];
        System.arraycopy(labc[k], 0, labcM[k], 0, indSwap[i] + 1);
        System.arraycopy(labc[k], indSwap[i], labcM[k], indSwap[i] + 1, labc[k].length - indSwap[i]);
        System.arraycopy(rabc[k], 0, rabcM[k], 0, indSwap[i] + 1);
        System.arraycopy(rabc[k], indSwap[i], rabcM[k], indSwap[i] + 1, rabc[k].length - indSwap[i]);
      }

      for (int looppt = 0; looppt < 4 * NB_POINT + 1; looppt++) {
        labcM[1][looppt] = labcM[1][looppt] + labcM[0][looppt] * 2 * ci;
        labcM[2][looppt] = labcM[2][looppt] + labcM[1][looppt] * ci - labcM[0][looppt] * ci * ci;
        rabcM[1][looppt] = rabcM[1][looppt] + rabcM[0][looppt] * 2 * ci;
        rabcM[2][looppt] = rabcM[2][looppt] + rabcM[1][looppt] * ci - rabcM[0][looppt] * ci * ci;
      }
      xroot = (-dabc[1] - Math.sqrt(dabc[1] * dabc[1] - 4 * dabc[0] * dabc[2])) / (2 * dabc[0]);

      ncdfinit = NORMAL.getCDF(xN[i][0]);

      for (int looppt = 0; looppt < nbPoint2; looppt++) {
        ncdf0[looppt] = NORMAL.getCDF(xN[i][looppt] - ci);
        ncdf1[looppt] = -c1sqrt2pi * Math.exp(-(xN[i][looppt] - ci) * (xN[i][looppt] - ci) / 2.0);
        ncdf2[looppt] = ncdf1[looppt] * (xN[i][looppt] - ci) + ncdf0[looppt];
      }

      for (int j = startInt[i]; j <= endInt[i]; j++) {
        is = indSwap[i] - j + 1;
        // % all L
        if (j + 2 * NB_POINT <= indSwap[i]) {
          final double[][] xabc = new double[3][2 * NB_POINT];
          for (int k = 0; k < 3; k++) {
            System.arraycopy(labcM[k], j, xabc[k], 0, 2 * NB_POINT);
          }
          for (int looppt = 0; looppt < 2 * NB_POINT; looppt++) {
            xabc[1][looppt] = xabc[1][looppt] + xabc[0][looppt] * 2 * xN[i][j];
            xabc[2][looppt] = xabc[2][looppt] + xabc[1][looppt] * xN[i][j] - xabc[0][looppt] * xN[i][j] * xN[i][j];
          }
          vW[i][j + NB_POINT] = 0;
          vW[i][j + NB_POINT] = vW[i][j + NB_POINT] + coi * ni2ncdf(ncdf2, ncdf1, ncdf0, xabc);
        } else if (j < indSwap[i]) {
          final double[][] xabc = new double[3][2 * NB_POINT + 1];
          tmpdb = xroot - xN[i][j] - ci;
          ncdf0x = NORMAL.getCDF(tmpdb);
          ncdf1x = -Math.exp(-(tmpdb * tmpdb) / 2) * c1sqrt2pi;
          ncdf2x = ncdf1x * tmpdb + ncdf0x;
          for (int k = 0; k < 3; k++) {
            //            System.arraycopy(rabcM[k], j, xabc[k], 0, 2 * _nbPoint + 1); // Swap
            System.arraycopy(labcM[k], j, xabc[k], 0, 2 * NB_POINT + 1);
            System.arraycopy(rabcM[k], indSwap[i] + 1, xabc[k], indSwap[i] + 1 - j, j + 2 * NB_POINT - indSwap[i]);
          }
          for (int looppt = 0; looppt < 2 * NB_POINT; looppt++) {
            xabc[1][looppt] = xabc[1][looppt] + xabc[0][looppt] * 2 * xN[i][j];
            xabc[2][looppt] = xabc[2][looppt] + xabc[1][looppt] * xN[i][j] - xabc[0][looppt] * xN[i][j] * xN[i][j];
          }
          System.arraycopy(ncdf0, 0, ncdf0X, 0, is);
          ncdf0X[is] = ncdf0x;
          System.arraycopy(ncdf0, is, ncdf0X, is + 1, ncdf0.length - is);
          System.arraycopy(ncdf1, 0, ncdf1X, 0, is);
          ncdf1X[is] = ncdf1x;
          System.arraycopy(ncdf1, is, ncdf1X, is + 1, ncdf1.length - is);
          System.arraycopy(ncdf2, 0, ncdf2X, 0, is);
          ncdf2X[is] = ncdf2x;
          System.arraycopy(ncdf2, is, ncdf2X, is + 1, ncdf2.length - is);
          vW[i][j + NB_POINT] = vL[j] * vZ[i][0] * ncdfinit;
          vW[i][j + NB_POINT] += coi * ni2ncdf(ncdf2X, ncdf1X, ncdf0X, xabc);
          vW[i][j + NB_POINT] += vR[j + 2 * NB_POINT] * vZ[i][vZ[i].length - 1] * ncdfinit;
        } else {
          final double[][] xabc = new double[3][2 * NB_POINT];
          for (int k = 0; k < 3; k++) {
            System.arraycopy(rabcM[k], j + 1, xabc[k], 0, 2 * NB_POINT);
            //            System.arraycopy(labcM[k], j + 1, xabc[k], 0, 2 * _nbPoint); // Swaption
          }
          for (int looppt = 0; looppt < 2 * NB_POINT; looppt++) {
            xabc[1][looppt] = xabc[1][looppt] + xabc[0][looppt] * 2 * xN[i][j];
            xabc[2][looppt] = xabc[2][looppt] + xabc[1][looppt] * xN[i][j] - xabc[0][looppt] * xN[i][j] * xN[i][j];
          }
          vW[i][j + NB_POINT] = vR[j] * vZ[i][0] * ncdfinit;
          vW[i][j + NB_POINT] += coi * ni2ncdf(ncdf2, ncdf1, ncdf0, xabc);
          vW[i][j + NB_POINT] += vR[j + 2 * NB_POINT] * vZ[i][vZ[i].length - 1] * ncdfinit;
        }
      }
      for (int j = 0; j < NB_POINT; j++) { // Flat extrapolation
        vW[i][j] = vW[i][NB_POINT];
        vW[i][3 * NB_POINT + 1 + j] = vW[i][3 * NB_POINT];
      }
    } // End main loop

    return MultipleCurrencyAmount.of(swaption.getUnderlyingSwap()[0].getFixedLeg().getCurrency(), vW[0][2 * NB_POINT] * (swaption.isLong() ? 1.0 : -1.0));
  }

  private  SwaptionBermudaFixedIborHullWhiteNumericalIntegrationMethod();

}