def mc_integrate(lower_lim, upper_lim, nsamps):
            """
            approximate the definite integral of `func` from
            `lower_lim` to `upper_lim`
            """
            out = cuda.to_device(np.zeros(nsamps, dtype="float32"))
            rng_states = create_xoroshiro128p_states(nsamps, seed=42)

            # jit the function for use in CUDA kernels

            mc_integrator_kernel.forall(nsamps)(
                out, rng_states, lower_lim, upper_lim
            )
            # normalization factor to convert
            # to the average: (b - a)/(N - 1)
            factor = (upper_lim - lower_lim) / (nsamps - 1)

            return sum_reduce(out) * factor

            approximate the definite integral of `func` from
            `lower_lim` to `upper_lim`
            
numba/cuda/tests/doc_examples/test_montecarlo.py