human
Complete this python program:
def _mpmath_kraft_burrows_nousek(N, B, CL):
assistant
#Here is the Python code you requested.

    """Upper limit on a poisson count rate
    The implementation is based on Kraft, Burrows and Nousek in
    `ApJ 374, 344 (1991) <http://adsabs.harvard.edu/abs/1991ApJ...374..344K>`_.
    The XMM-Newton upper limit server used the same formalism.
    Parameters
    N : int
    Total observed count number
    B : float
    Background count rate (assumed to be known with negligible error
    from a large background area).
    CL : float
    Confidence level (number between 0 and 1)
    Returns
    S : source count limit
    Notes
    Requires the `mpmath <http://mpmath.org/>`_ library.  See
    `~astropy.stats.scipy_poisson_upper_limit` for an implementation
    that is based on scipy and evaluates faster, but runs only to about
    N = 100.
    """
    from mpmath import mpf, factorial, findroot, fsum, power, exp, quad
    N = mpf(N)
    B = mpf(B)
    CL = mpf(CL)
    def eqn8(N, B):
        sumterms = [(power(B, n) / factorial(n)) for n in range((int(N) + 1))]
        return (1.0 / (exp((- B)) * fsum(sumterms)))
    eqn8_res = eqn8(N, B)
    factorial_N = factorial(N)
    def eqn7(S, N, B):
        SpB = (S + B)
        return (eqn8_res * ((exp((- SpB)) * (SpB ** N)) / factorial_N))
    def eqn9_left(S_min, S_max, N, B):
        def eqn7NB(S):
            return eqn7(S, N, B)
        return quad(eqn7NB, [S_min, S_max])
    def find_s_min(S_max, N, B):
        u'\n                                Kraft,    Burrows    and    Nousek    suggest    to    integrate    from    N-B    in    both\n                                directions    at    once,    so    that    S_min    and    S_max    move    similarly    (see\n                                the    article    for    details).    Here,    this    is    implemented    differently:\n                                Treat    S_max    as    the    optimization    parameters    in    func    and    then\n                                calculate    the    matching    s_min    that    has    has    eqn7(S_max)    =\n                                eqn7(S_min)    here.\n                                '
        y_S_max = eqn7(S_max, N, B)
        if (eqn7(0, N, B) >= y_S_max):
            return 0.0
        else:
            def eqn7ysmax(x):
                return (eqn7(x, N, B) - y_S_max)
            return findroot(eqn7ysmax, ((N - B) / 2.0))
    def func(s):
        s_min = find_s_min(s, N, B)
        out = eqn9_left(s_min, s, N, B)
        return (out - CL)
    S_max = findroot(func, (N - B), tol=0.0001)
    S_min = find_s_min(S_max, N, B)
    return (float(S_min), float(S_max))