векторизовать неполную гамма-функцию с отрицательными аргументами

Я использую неполную гамма-функцию mpmaths с отрицательным z и сложными границами интегрирования a, b. См. документацию: документы

Это невозможно сделать с неполной гамма-функцией scipy. Я сравнил отмеченный mpmath с Mathematica и получил те же результаты. К сожалению, функция mpmaths может оценивать только одно значение за раз, а это значит, что мне приходится перебирать массив, содержащий значения z, a, b.

Я считаю, что векторизация сложна, потому что для отрицательного z нужно рекурсивно оценивать интеграл по частям (см. этот пост и соответствующую вики: неполная гамма-функция в python?)

Кто-нибудь знает модуль, который может справиться с этим? Или если такого нет, если это вообще возможно?

Итак по сути:

from mpmath import gammainc
values = np.random.normal(0, 100, 1000) + 1j*np.random.normal(0, 100, 1000)
res_mp = np.array([gammainc(cc, 4, 10) for cc in values])

На самом деле я реализовал то, что, по моему мнению, работает в некоторой степени:

import numpy as np

def new_quad_routine(func, c, a, b,x_list, w_list):
    c_1 = (b-a)/2.0
    c_2 = (b+a)/2.0

    eval_points = c_1*x_list+c_2             
    func_evals = func(eval_points, c)                              
    return c_1 * np.sum(func_evals * w_list[:,np.newaxis], axis=0)

def new_quad_gauss_7(func, c, a, b):   
    """
    integrates a complex function with bounds, a, b with an array of arguments c
    call for instance with the gamma function:
    new_quad_gauss_7(gamma_integrator, 4, 10)

    """
    x_gauss = np.array([-0.949107912342759, -0.741531185599394, -0.405845151377397, 0, 0.405845151377397, 0.741531185599394, 0.949107912342759])            
    w_gauss = np.array([0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, 0.381830050505119, 0.279705391489277,0.129484966168870])
    return new_quad_routine(func, c, a, b, x_gauss, w_gauss)

def gamma_integrator(t, c):
    return t[:, np.newaxis]**c*np.exp(-t[:,np.newaxis])

def gammainc_function(c, a, b):
    if type(c) is not np.ndarray:
        raise ValueError("Need a numpy array for function argument!")
    return new_quad_gauss_7(gamma_integrator, c-1, a, b)

Это основано на численном интеграторе в этом посте.

Это довольно быстро:

values=np.repeat(-2+3j-1, 100)
In [143]: %timeit new_quad_gauss_7(gamma_integrator, 4, 10, values)
107 µs ± 574 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [144]: %timeit [gammainc(cc, 4, 10) for cc in values]
224 ms ± 2.86 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Хотя не знаю, насколько это хорошо. Итак, вопрос немного меняется на то, как я могу сравнить и/или улучшить это? Кроме того, это, кажется, работает для фиксированных границ. Однако я не знаю, как реализовать интеграл с бесконечными границами на обоих концах.