docs: def _log_phi(z):
      """Stable computation of the log of the Normal CDF and its derivative."""
      # Adapted from the GPML function `logphi.m`.
      if z * z < 0.0492:
          # First case: z close to zero.
          coef = -z / SQRT2PI
          val = functools.reduce(lambda acc, c: coef * (c + acc), CS, 0)
          res = -2 * val - log(2)
          dres = exp(-(z * z) / 2 - res) / SQRT2PI
      elif z < -11.3137:
          # Second case: z very small.
          num = functools.reduce(
                  lambda acc, r: -z * acc / SQRT2 + r, RS, 0.5641895835477550741)
          den = functools.reduce(lambda acc, q: -z * acc / SQRT2 + q, QS, 1.0)
          res = log(num / (2 * den)) - (z * z) / 2
          dres = abs(den / num) * sqrt(2.0 / pi)
      else:
          res = log(normal_cdf(z))
          dres = exp(-(z * z) / 2 - res) / SQRT2PI
      return res, dres