class BlackFormulaRepository {

  public static Pair<ValueDerivatives, double[][]> priceAdjoint2(
      double forward,
      double strike,
      double timeToExpiry,
      double lognormalVol,
      boolean isCall) {
    // Forward sweep
    double discountFactor = 1.0;
    double sqrttheta = Math.sqrt(timeToExpiry);
    double omega = isCall ? 1 : -1;
    // Implementation Note: Forward sweep.
    double volPeriod = 0, kappa = 0, d1 = 0, d2 = 0;
    double x = 0;
    double p;
    if (strike < NEAR_ZERO || sqrttheta < NEAR_ZERO) {
      x = omega * (forward - strike);
      p = (x > 0 ? discountFactor * x : 0.0);
      volPeriod = sqrttheta < NEAR_ZERO ? 0 : (lognormalVol * sqrttheta);
    } else {
      volPeriod = lognormalVol * sqrttheta;
      kappa = Math.log(forward / strike) / volPeriod - 0.5 * volPeriod;
      d1 = NORMAL.getCDF(omega * (kappa + volPeriod));
      d2 = NORMAL.getCDF(omega * kappa);
      p = discountFactor * omega * (forward * d1 - strike * d2);
    }
    // Implementation Note: Backward sweep.
    double[][] bsD2 = new double[3][3];
    double pBar = 1.0;
    double density1 = 0.0;
    double d1Bar = 0.0;
    double forwardBar = 0, strikeBar = 0, volPeriodBar = 0, lognormalVolBar = 0, sqrtthetaBar = 0, timeToExpiryBar = 0;
    if (strike < NEAR_ZERO || sqrttheta < NEAR_ZERO) {
      forwardBar = (x > 0 ? discountFactor * omega : 0.0);
      strikeBar = (x > 0 ? -discountFactor * omega : 0.0);
    } else {
      d1Bar = discountFactor * omega * forward * pBar;
      density1 = NORMAL.getPDF(omega * (kappa + volPeriod));
      // Implementation Note: kappa_bar = 0; no need to implement it.
      // Methodology Note: kappa_bar is optimal exercise boundary. The
      // derivative at the optimal point is 0.
      forwardBar = discountFactor * omega * d1 * pBar;
      strikeBar = -discountFactor * omega * d2 * pBar;
      volPeriodBar = density1 * omega * d1Bar;
      lognormalVolBar = sqrttheta * volPeriodBar;
      sqrtthetaBar = lognormalVol * volPeriodBar;
      timeToExpiryBar = 0.5 / sqrttheta * sqrtthetaBar;
    }
    DoubleArray bsD = DoubleArray.of(forwardBar, strikeBar, timeToExpiryBar, lognormalVolBar);
    if (strike < NEAR_ZERO || sqrttheta < NEAR_ZERO) {
      return Pair.of(ValueDerivatives.of(p, bsD), bsD2);
    }
    // Backward sweep: second derivative
    double d2Bar = -discountFactor * omega * strike;
    double density2 = NORMAL.getPDF(omega * kappa);
    double d1Kappa = omega * density1;
    double d1KappaKappa = -(kappa + volPeriod) * d1Kappa;
    double d2Kappa = omega * density2;
    double d2KappaKappa = -kappa * d2Kappa;
    double kappaKappaBar2 = d1KappaKappa * d1Bar + d2KappaKappa * d2Bar;
    double kappaV = -Math.log(forward / strike) / (volPeriod * volPeriod) - 0.5;
    double kappaVV = 2 * Math.log(forward / strike) / (volPeriod * volPeriod * volPeriod);
    double d1TotVV = density1 * omega * (-(kappa + volPeriod) * (kappaV + 1) * (kappaV + 1) + kappaVV);
    double d2TotVV = d2KappaKappa * kappaV * kappaV + d2Kappa * kappaVV;
    double vVbar2 = d1Bar * d1TotVV + d2Bar * d2TotVV;
    double volVolBar2 = vVbar2 * timeToExpiry;
    double kappaStrikeBar2 = -discountFactor * omega * d2Kappa;
    double kappaStrike = -1.0 / (strike * volPeriod);
    double strikeStrikeBar2 = (kappaKappaBar2 * kappaStrike + 2 * kappaStrikeBar2) * kappaStrike;
    double kappaStrikeV = 1.0 / strike / (volPeriod * volPeriod);
    double d1VK = -omega * (kappa + volPeriod) * density1 * (kappaV + 1) * kappaStrike + omega * density1 * kappaStrikeV;
    double d2V = d2Kappa * kappaV;
    double d2VK = -omega * kappa * density2 * kappaV * kappaStrike + omega * density2 * kappaStrikeV;
    double strikeD2Bar2 = -discountFactor * omega;
    double strikeVolblackBar2 = strikeD2Bar2 * d2V + d1Bar * d1VK + d2Bar * d2VK;
    double strikeVolBar2 = strikeVolblackBar2 * sqrttheta;
    double kappaForward = 1.0 / (forward * volPeriod);
    double forwardForwardBar2 = discountFactor * omega * d1Kappa * kappaForward;
    double strikeForwardBar2 = discountFactor * omega * d1Kappa * kappaStrike;
    double volForwardBar2 = discountFactor * omega * d1Kappa * (kappaV + 1) * sqrttheta;
    bsD2[0][0] = forwardForwardBar2;
    bsD2[0][2] = volForwardBar2;
    bsD2[2][0] = volForwardBar2;
    bsD2[0][1] = strikeForwardBar2;
    bsD2[1][0] = strikeForwardBar2;
    bsD2[2][2] = volVolBar2;
    bsD2[1][2] = strikeVolBar2;
    bsD2[2][1] = strikeVolBar2;
    bsD2[1][1] = strikeStrikeBar2;
    return Pair.of(ValueDerivatives.of(p, bsD), bsD2);
  }

  private  BlackFormulaRepository();

}