mpmath icon indicating copy to clipboard operation
mpmath copied to clipboard

Bell incomplete polynomials

Open lrnv opened this issue 5 years ago • 0 comments

Hi,

The library (which is perfect btw) provides the bell completes polynomials.

Do you think you could also add incompletes ones ? The following code is an exemple of possible implementation, based on https://github.com/antiphon/BellB/blob/master/R/bellB.R and wikipedia :

import numpy as np
import mpmath as mp

def __rf(m, l, x):
    if l < 1:
        return x[m - 1]
    if (l == m - 1):
        return mp.binomial(m, l) * x[m - l - 1] * __rf(l, l - 1, x)
    else:
        rez = mp.mpf(0)
        for i in range(l, m):
            rez = rez + mp.binomial(m, i) * x[m - i - 1] * __rf(i, l - 1, x)
        return rez

def __bell_incomplete(n, k, x):
    assert (len(x) == n - k + 1)
    if (n == 0):
        return [mp.mpf(1) for i in np.arange(k)]
    K = k - 1
    v = __rf(n, k - 1, x)
    return v / mp.factorial(k)

def __bell_complete(n, x):
    if (n == 0):
        return 1
    else:
        if (n == 1):
            return __bell_incomplete(n, 1, x[:n])
        else:
            rez = 0
            for k in np.arange(n) + 1:
                rez = rez + __bell_incomplete(n, k, x[:(n - k + 1)])

This code is fairly not perfect, but does the job. A good test will be to check that this function returns the same as mp.bell for complete bell polynomials.

lrnv avatar Oct 15 '19 14:10 lrnv