code: def _cbnd(dh, dk, rho):
      """Computes values for the cumulative standard bivariate normal distribution.
  
      More specifically, compultes `P(x > dh, y > dk)` where `x` and `y` are
      standard normal variables with correlation `rho`.
  
      Args:
        dh: A real `Tensor` representing lower integration limits for `x`.
        dk: A `Tensor` of the same dtype as `dh` and of compatible shape
          representing lower integration limits `y`.
        rho: A `Tensor` of the same dtype as `dh` and of compatible shape
          representing correlation coefficients.
  
      Returns:
        A `Tensor` of cumulative distribution function values.
  
      #### References:
      [1] Genz, A., Numerical Computation of Rectangular Bivariate and Trivariate
            Normal and t Probabilities, 2004
        http://www.math.wsu.edu/faculty/genz/papers/bvnt.pdf
      """
  
      # Python code adapted from Fortran code and naming convention will align
      # variables with the variables named in BVND function.
      # http://www.math.wsu.edu/faculty/genz/software/fort77/tvpack.f
  
      dtype = rho.dtype
  
      h = tf.cast(-dh, dtype=dtype)
      k = tf.cast(-dk, dtype=dtype)
      hk = h * k
  
      # In order to avoid `shape_invariants` argument in the while_loop, set
      # `bvn` to be a Tensor of zeros of the appropriate shape. The shape
      # corresponds to the number of options being priced.
      bvn = tf.zeros_like(hk)
  
      hs = (h * h + k * k) / 2
      asr = tf.math.asin(rho)
  
      # Uses Gauss-Legendre quadrature to estimate the integral.
      # Unlike in the reference paper, always uses 20 points regardless of value
      # of rho, for efficiency reasons on the GPU.
      def transformed_bvn(hk, hs, asr):
          def transformed_bvn_distribution(x):
              # Shape bvn.shape + [1]
              hk_exp = tf.expand_dims(hk, axis=-1)
              hs_exp = tf.expand_dims(hs, axis=-1)
              asr_exp = tf.expand_dims(asr, axis=-1)
              # Shape bvn.shape + x.shape
              sn1 = tf.math.sin(asr_exp * (1 * x + 1) / 2)
              return tf.math.exp(
                  (hk_exp * sn1 - hs_exp) / (1 - sn1 * sn1)
              )
  
          # bvn.shape
          ones = tf.ones_like(hk)
          res = integration.gauss_legendre(
              func=transformed_bvn_distribution,
              lower=-ones,
              upper=ones,
              num_points=20,
              dtype=dtype,
          )
          # bvn.shape
          return res
  
      bvn = bvn + transformed_bvn(hk, hs, asr)
      bvn = bvn * asr / (4 * np.pi)
      bvn = bvn + _ncdf(-h) * _ncdf(-k)
  
      return bvn
  
imports: from typing import Tuple
  import numpy
  import tensorflow.compat.v2
  from tf_quant_finance import types
  from tf_quant_finance.black_scholes import vanilla_prices
  from tf_quant_finance.math import gradient
  from tf_quant_finance.math import integration
  from tf_quant_finance.math.root_search import newton
label_desc: Use image node_id 16 for calling a global function with example usage: _cbnd(dh, dk, rho) and returns: bvn, res, tf