Is there any good way to calculate cumulative density function of a binary normal distribution

Question

$I(x, y, dx, dy) = \int^{x+dx}_{x}\int^{y+dy}_{y} e^{-\frac{1}{2}(ax^2+by^2+2cxy)}$, how can I estimate it with as little computation as possible?

Besides, I wonder if it is always hard to calculate or estimate the integration of $f(|ax^2+by^2+2cxy|)$ where $f(\cdot)$ is a monotonically decreasing function.