def _impl(n, p):
            """
            Binomial distribution.  Numpy's variant of the BINV algorithm
            is used.
            (Numpy uses BTPE for n*p >= 30, though)
            """
            if n < 0:
                raise ValueError("binomial(): n <= 0")
            if not (0.0 <= p <= 1.0):
                raise ValueError("binomial(): p outside of [0, 1]")
            if p == 0.0:
                return 0
            if p == 1.0:
                return n

            flipped = p > 0.5
            if flipped:
                p = 1.0 - p
            q = 1.0 - p

            niters = 1
            qn = q ** n
            while qn <= 1e-308:
                # Underflow => split into several iterations
                # Note this is much slower than Numpy's BTPE
                niters <<= 2
                n >>= 2
                qn = q ** n
                assert n > 0

            np_prod = n * p
            bound = min(n, np_prod + 10.0 * math.sqrt(np_prod * q + 1))

            total = 0
            while niters > 0:
                X = 0
                U = np.random.random()
                px = qn
                while X <= bound:
                    if U <= px:
                        total += n - X if flipped else X
                        niters -= 1
                        break
                    U -= px
                    X += 1
                    px = ((n - X + 1) * p * px) / (X * q)

            return total

            Binomial distribution.  Numpy's variant of the BINV algorithm
            is used.
            (Numpy uses BTPE for n*p >= 30, though)
            
numba/cpython/randomimpl.py