Function bivariate_norm_cdf 

pub fn bivariate_norm_cdf(a: f64, b: f64, rho: f64) -> f64

Cumulative bivariate normal distribution P(X <= a, Y <= b) for standard normals with correlation rho.

Genz (2004) Gauss-Legendre quadrature on the arcsine form for |rho| <= 0.925 (accurate to ~1e-14); an adaptive-free Simpson integration of phi(x) N((b - rho x)/sqrt(1-rho^2)) for the highly correlated tail. Used by the Bjerksund-Stensland (2002) two-boundary American approximation.