Python scipy.special.factorial() Examples

def test_odd_cat_state(self, tol):
        """test correct even cat state returned"""
        a = 0.212
        cutoff = 10
        p = 1

        state = utils.cat_state(a, p, fock_dim=cutoff)

        n = np.arange(cutoff)
        expected = np.exp(-0.5 * np.abs(a) ** 2) * a ** n / np.sqrt(fac(n)) - np.exp(
            -0.5 * np.abs(-a) ** 2
        ) * (-a) ** n / np.sqrt(fac(n))
        expected /= np.linalg.norm(expected)

        assert np.allclose(state, expected, atol=tol, rtol=0)


# ===================================================================================
# Random matrix tests
# =================================================================================== 

def _alphaK_T_xs(self, temperature, volume, params): # eq. 3.21 of de Koker thesis
        f = self._finite_strain(temperature, volume, params)
        theta = self._theta(temperature, volume, params)
        sum_factors = 0.
        for i in range(len(params['a'])):
            ifact=factorial(i, exact=False)
            if i > 0:
                for j in range(len(params['a'][0])):
                    if j > 0:
                        jfact=factorial(j, exact=False)
                        sum_factors += float(i)*float(j)*params['a'][i][j] \
                            * np.power(f, float(i-1)) * np.power(theta, float(j-1)) \
                            / ifact / jfact
                            
        return -self._dfdV(temperature, volume, params) \
            * self._dthetadT(temperature, volume, params) \
            * sum_factors 

def _K_T_xs(self, temperature, volume, params): # K_T_xs, eq. 3.20 of de Koker thesis
        f = self._finite_strain(temperature, volume, params)
        theta = self._theta(temperature, volume, params)
        K_ToverV=0.
        for i in range(len(params['a'])):
            ifact=factorial(i, exact=False)
            for j in range(len(params['a'][0])):
                if i > 0:
                    jfact=factorial(j, exact=False)
                    prefactor = float(i) * params['a'][i][j] \
                        * np.power(theta, float(j)) / ifact / jfact 
                    K_ToverV += prefactor*self._d2fdV2(temperature, volume, params) \
                        * np.power(f, float(i-1))
                if i > 1:
                    dfdV = self._dfdV(temperature, volume, params)
                    K_ToverV += prefactor * dfdV * dfdV \
                        * float(i-1) * np.power(f, float(i-2))
        return volume*K_ToverV 

def F_q2d(n, m):
    if n == 0:
        num = m ** 2 * factorial2(2 * m - 3)
        den = 2 ** (m + 1) * factorial(m - 1)
        return num / den
    elif n > 0 and m == 1:
        t1num = 4 * (n - 1) ** 2 * n ** 2 + 1
        t1den = 8 * (2 * n - 1) ** 2
        term1 = t1num / t1den
        term2 = 11 / 32 * kronecker(n, 1)
        return term1 + term2
    else:
        Chi = m + n - 2
        nt1 = 2 * n * Chi * (3 - 5 * m + 4 * n * Chi)
        nt2 = m ** 2 * (3 - m + 4 * n * Chi)
        num = nt1 + nt2

        dt1 = (m + 2 * n - 3) * (m + 2 * n - 2)
        dt2 = (m + 2 * n - 1) * (2 * n - 1)
        den = dt1 * dt2

        term1 = num / den
        return term1 * gamma(n, m) 

def _get_M(self, delta_t):
        n = len(self.a)
        A = np.zeros(n)
        for i, ai in enumerate(self.a):
            Ae = [self.a[i] * (-1)**j * factorial(i) / factorial(j) / factorial(
                i-j) / ((delta_t)**i) for j in range(i + 1)]  # Elementary A to assemblate in A
            for j, aej in enumerate(Ae):
                A[j] += aej

        n = len(self.b)
        B = np.zeros(n)
        for i, ai in enumerate(self.b):
            Be = [self.b[i] * (-1)**j * factorial(i) / factorial(j) / factorial(
                i-j) / ((delta_t)**i) for j in range(i + 1)]  # Elementary B to assemblate in B
            for j, bej in enumerate(Be):
                B[j] += bej

        Mo = [-x/B[0] for x in B[1:][::-1]]
        Mi = [x/B[0] for x in A[::-1]]
        return (Mi, Mo) 

def _falling_factorial(x, n):
    r"""
    Return the factorial of `x` to the `n` falling.

    This is defined as:

    .. math::   x^\underline n = (x)_n = x (x-1) \cdots (x-n+1)

    This can more efficiently calculate ratios of factorials, since:

    n!/m! == falling_factorial(n, n-m)

    where n >= m

    skipping the factors that cancel out

    the usual factorial n! == ff(n, n)
    """
    val = 1
    for k in range(x - n + 1, x + 1):
        val *= k
    return val 

def G_q2d(n, m):
    if n == 0:
        num = factorial2(2 * m - 1)
        den = 2 ** (m + 1) * factorial(m - 1)
        return num / den
    elif n > 0 and m == 1:
        t1num = (2 * n ** 2 - 1) * (n ** 2 - 1)
        t1den = 8 * (4 * n ** 2 - 1)
        term1 = -t1num / t1den
        term2 = 1 / 24 * kronecker(n, 1)
        return term1 + term2  # this is minus in the paper
    else:
        # nt1 = numerator term 1, d = denominator...
        nt1 = 2 * n * (m + n - 1) - m
        nt2 = (n + 1) * (2 * m + 2 * n - 1)
        num = nt1 * nt2
        dt1 = (m + 2 * n - 2) * (m + 2 * n - 1)
        dt2 = (m + 2 * n) * (2 * n + 1)
        den = dt1 * dt2

        term1 = num / den  # there is a leading negative in the paper
        return term1 * gamma(n, m) 

def __init__(self, xi, yi, axis=0):
        _Interpolator1DWithDerivatives.__init__(self, xi, yi, axis)

        self.xi = np.asarray(xi)
        self.yi = self._reshape_yi(yi)
        self.n, self.r = self.yi.shape

        c = np.zeros((self.n+1, self.r), dtype=self.dtype)
        c[0] = self.yi[0]
        Vk = np.zeros((self.n, self.r), dtype=self.dtype)
        for k in xrange(1,self.n):
            s = 0
            while s <= k and xi[k-s] == xi[k]:
                s += 1
            s -= 1
            Vk[0] = self.yi[k]/float(factorial(s))
            for i in xrange(k-s):
                if xi[i] == xi[k]:
                    raise ValueError("Elements if `xi` can't be equal.")
                if s == 0:
                    Vk[i+1] = (c[i]-Vk[i])/(xi[i]-xi[k])
                else:
                    Vk[i+1] = (Vk[i+1]-Vk[i])/(xi[i]-xi[k])
            c[k] = Vk[k-s]
        self.c = c 

def imp_lsc(self, gamma, sigma, w, dz):
        """
        gamma - energy
        sigma - transverse RMS size of the beam
        w - omega = 2*pi*f
        """
        eps = 1e-16
        ass = 40.0

        alpha = w * sigma / (gamma * speed_of_light)
        alpha2 = alpha * alpha

        inda = np.where(alpha2 > ass)[0]
        ind = np.where((alpha2 <= ass) & (alpha2 >= eps))[0]

        T = np.zeros(w.shape)
        T[ind] = np.exp(alpha2[ind]) *exp1(alpha2[ind])

        x= alpha2[inda]
        k = 0
        for i in range(10):
            k += (-1) ** i * factorial(i) / (x ** (i + 1))
        T[inda] = k
        Z = 1j * Z0 / (4 * pi * speed_of_light*gamma**2) * w * T * dz
        return Z # --> Omm/m 

def torontonian_analytical(l, nbar):
    r"""Return the value of the Torontonian of the O matrices generated by gen_omats

    Args:
        l (int): number of modes
        nbar (float): mean photon number of the first mode (the only one not prepared in vacuum)

    Returns:
        float: Value of the torontonian of gen_omats(l,nbar)
    """
    if np.allclose(l, nbar, atol=1e-14, rtol=0.0):
        return 1.0
    beta = -(nbar / (l * (1 + nbar)))
    pref = factorial(l) / beta
    p1 = pref * l / poch(1 / beta, l + 2)
    p2 = pref * beta / poch(2 + 1 / beta, l)
    return (p1 + p2) * (-1) ** l 

def _falling_factorial(x, n):
    r"""
    Return the factorial of `x` to the `n` falling.

    This is defined as:

    .. math::   x^\underline n = (x)_n = x (x-1) \cdots (x-n+1)

    This can more efficiently calculate ratios of factorials, since:

    n!/m! == falling_factorial(n, n-m)

    where n >= m

    skipping the factors that cancel out

    the usual factorial n! == ff(n, n)
    """
    val = 1
    for k in range(x - n + 1, x + 1):
        val *= k
    return val 

def H(n, x):
    """Explicit expression for Hermite polynomials."""
    prefactor = factorial(n)
    terms = tf.reduce_sum(
        tf.stack(
            [
                _numer_safe_power(-1, m)
                / (factorial(m) * factorial(n - 2 * m))
                * _numer_safe_power(2 * x, n - 2 * m)
                for m in range(int(np.floor(n / 2)) + 1)
            ],
            axis=0,
        ),
        axis=0,
    )
    return prefactor * terms 

def __init__(self, xi, yi, axis=0):
        _Interpolator1DWithDerivatives.__init__(self, xi, yi, axis)

        self.xi = np.asarray(xi)
        self.yi = self._reshape_yi(yi)
        self.n, self.r = self.yi.shape

        c = np.zeros((self.n+1, self.r), dtype=self.dtype)
        c[0] = self.yi[0]
        Vk = np.zeros((self.n, self.r), dtype=self.dtype)
        for k in xrange(1,self.n):
            s = 0
            while s <= k and xi[k-s] == xi[k]:
                s += 1
            s -= 1
            Vk[0] = self.yi[k]/float(factorial(s))
            for i in xrange(k-s):
                if xi[i] == xi[k]:
                    raise ValueError("Elements if `xi` can't be equal.")
                if s == 0:
                    Vk[i+1] = (c[i]-Vk[i])/(xi[i]-xi[k])
                else:
                    Vk[i+1] = (Vk[i+1]-Vk[i])/(xi[i]-xi[k])
            c[k] = Vk[k-s]
        self.c = c 

def test_loss_channel_on_coherent_states(
        self, setup_backend, T, mag_alpha, phase_alpha, cutoff, tol
    ):
        """Tests various loss channels on coherent states (result should be coherent state with amplitude weighted by sqrt(T)."""

        rootT_alpha = np.sqrt(T) * mag_alpha * np.exp(1j * phase_alpha)

        backend = setup_backend(1)

        backend.prepare_coherent_state(mag_alpha, phase_alpha, 0)
        backend.loss(T, 0)
        s = backend.state()
        if s.is_pure:
            numer_state = s.ket()
        else:
            numer_state = s.dm()
        n = np.arange(cutoff)
        ref_state = (
            np.exp(-0.5 * np.abs(rootT_alpha) ** 2)
            * rootT_alpha ** n
            / np.sqrt(factorial(n))
        )
        ref_state = np.outer(ref_state, np.conj(ref_state))
        assert np.allclose(numer_state, ref_state, atol=tol, rtol=0.0) 

def test_displaced_squeezed_state_fock(self, r_d, phi_d, r_s, phi_s, hbar, cutoff, tol):
        """test displaced squeezed state returns correct Fock basis state vector"""
        state = utils.displaced_squeezed_state(r_d, phi_d, r_s, phi_s, basis="fock", fock_dim=cutoff, hbar=hbar)
        a = r_d * np.exp(1j * phi_d)

        if r_s == 0:
            pytest.skip("test only non-zero squeezing")

        n = np.arange(cutoff)
        gamma = a * np.cosh(r_s) + np.conj(a) * np.exp(1j * phi_s) * np.sinh(r_s)
        coeff = np.diag(
            (0.5 * np.exp(1j * phi_s) * np.tanh(r_s)) ** (n / 2) / np.sqrt(fac(n) * np.cosh(r_s))
        )

        expected = H(gamma / np.sqrt(np.exp(1j * phi_s) * np.sinh(2 * r_s)), coeff)
        expected *= np.exp(
            -0.5 * np.abs(a) ** 2 - 0.5 * np.conj(a) ** 2 * np.exp(1j * phi_s) * np.tanh(r_s)
        )

        assert np.allclose(state, expected, atol=tol, rtol=0) 

def orbit_cardinality(orbit: list, modes: int) -> int:
    """Gives the number of samples belonging to the input orbit.

    For example, there are three possible samples in the orbit [2, 1, 1] with three modes: [1, 1,
    2], [1, 2, 1], and [2, 1, 1]. With four modes, there are 12 samples in total.

    **Example usage:**

    >>> orbit_cardinality([2, 1, 1], 4)
    12

    Args:
        orbit (list[int]): orbit; we count how many samples are contained in it
        modes (int): number of modes in the samples

    Returns:
        int: number of samples in the orbit
    """
    sample = orbit + [0] * (modes - len(orbit))
    counts = list(Counter(sample).values())
    return int(factorial(modes) / np.prod(factorial(counts))) 

def test_nongaussian(self, a, phi, setup_backend, cutoff, tol):
        """Tests that probabilities of particular Fock states |n> are
        correct for a nongaussian state."""
        backend = setup_backend(2)

        alpha = a * np.exp(1j * phi)
        n = np.arange(cutoff)
        ref_state = np.exp(-0.5 * np.abs(alpha) ** 2) * alpha ** n / np.sqrt(fac(n))
        ref_probs = np.abs(ref_state) ** 2

        backend.prepare_coherent_state(np.abs(alpha), np.angle(alpha), 0)
        backend.prepare_fock_state(cutoff // 2, 1)
        state = backend.state()

        for n in range(cutoff):
            prob_n = state.fock_prob([n, cutoff // 2])
            assert np.allclose(prob_n, ref_probs[n], atol=tol, rtol=0) 

def test_pure(self, a, phi, setup_backend, cutoff, batch_size, tol):
        """Tests that the numeric probabilities in the full Fock basis are
        correct for a one-mode pure state."""
        backend = setup_backend(1)

        alpha = a * np.exp(1j * phi)
        n = np.arange(cutoff)
        ref_state = np.exp(-0.5 * np.abs(alpha) ** 2) * alpha ** n / np.sqrt(fac(n))
        ref_probs = np.abs(ref_state) ** 2

        backend.prepare_coherent_state(np.abs(alpha), np.angle(alpha), 0)
        state = backend.state()

        probs = state.all_fock_probs(cutoff=cutoff)
        if isinstance(probs, tf.Tensor):
            probs = probs.numpy()
        probs = probs.flatten()

        if batch_size is not None:
            ref_probs = np.tile(ref_probs, batch_size)

        assert np.allclose(probs, ref_probs, atol=tol, rtol=0) 

def test_reduced_state_fock_probs(self, cutoff, setup_backend, batch_size, tol):
        """Test backend calculates correct fock prob of reduced coherent state"""
        backend = setup_backend(2)

        backend.prepare_coherent_state(np.abs(a), np.angle(a), 0)
        backend.prepare_squeezed_state(r, phi, 1)

        state = backend.state(modes=[0])
        probs = np.array([state.fock_prob([i]) for i in range(cutoff)]).T

        n = np.arange(cutoff)
        ref_state = np.exp(-0.5 * np.abs(a) ** 2) * a ** n / np.sqrt(fac(n))
        ref_probs = np.abs(ref_state) ** 2

        if batch_size is not None:
            ref_probs = np.tile(ref_probs, batch_size)

        assert np.allclose(probs.flatten(), ref_probs.flatten(), atol=tol, rtol=0) 

def test_rdm(self, setup_backend, tol, cutoff, batch_size):
        """Test reduced density matrix of a coherent state is as expected
        This is supported by all backends, since it returns
        the reduced density matrix of a single mode."""
        backend = setup_backend(2)
        backend.prepare_coherent_state(np.abs(a), np.angle(a), 0)
        backend.prepare_coherent_state(0.1, 0, 1)

        state = backend.state()
        rdm = state.reduced_dm(0, cutoff=cutoff)

        n = np.arange(cutoff)
        ket = np.exp(-0.5 * np.abs(a) ** 2) * a ** n / np.sqrt(fac(n))
        rdm_exact = np.outer(ket, ket.conj())

        if batch_size is not None:
            np.tile(rdm_exact, [batch_size, 1])

        assert np.allclose(rdm, rdm_exact, atol=tol, rtol=0) 

def test_coherent_state_fock_elements(self, setup_backend, r, p, cutoff, pure, tol):
        r"""Tests if a range of alpha-displaced states have the correct Fock basis elements:
           |\alpha> = exp(-0.5 |\alpha|^2) \sum_n \alpha^n / \sqrt{n!} |n>
        """

        backend = setup_backend(1)
        backend.displacement(r, p, 0)
        state = backend.state()

        if state.is_pure:
            numer_state = state.ket()
        else:
            numer_state = state.dm()

        n = np.arange(cutoff)
        ref_state = np.exp(-0.5 * r ** 2) * (r*np.exp(1j*p)) ** n / np.sqrt(fac(n))

        if not pure:
            ref_state = np.outer(ref_state, np.conj(ref_state))

        assert np.allclose(numer_state, ref_state, atol=tol, rtol=0) 

def _naive_rsh_ph(l, m, theta, phi):

    if m == 0:
        return np.sqrt((2 * l + 1.) / (4 * np.pi)) * lpmv(m, l, np.cos(theta))
    elif m < 0:
        return np.sqrt(((2 * l + 1.) * factorial(l + m)) /
                       (2 * np.pi * factorial(l - m))) * lpmv(-m, l, np.cos(theta)) * np.sin(-m * phi)
    elif m > 0:
        return np.sqrt(((2 * l + 1.) * factorial(l - m)) /
                       (2 * np.pi * factorial(l + m))) * lpmv(m, l, np.cos(theta)) * np.cos(m * phi) 

def test_displaced_squeezed_with_no_displacement(
        self, setup_backend, r, phi, cutoff, batch_size, pure, tol
    ):
        """Tests if a squeezed coherent state with no displacement is equal to a squeezed state (Eq. (5.5.6) in Loudon)."""
        mag_alpha = 0
        phase_alpha = 0
        backend = setup_backend(1)

        backend.prepare_displaced_squeezed_state(mag_alpha, phase_alpha, r, phi, 0)
        state = backend.state()

        if state.is_pure:
            num_state = state.ket()
        else:
            num_state = state.dm()

        n = np.arange(0, cutoff, 2)
        even_refs = (
            np.sqrt(sech(r))
            * np.sqrt(factorial(n))
            / factorial(n / 2)
            * (-0.5 * np.exp(1j * phi) * np.tanh(r)) ** (n / 2)
        )

        if batch_size is not None:
            if pure:
                even_entries = num_state[:, ::2]
            else:
                even_entries = num_state[:, ::2, ::2]
                even_refs = np.outer(even_refs, np.conj(even_refs))
        else:
            if pure:
                even_entries = num_state[::2]
            else:
                even_entries = num_state[::2, ::2]
                even_refs = np.outer(even_refs, np.conj(even_refs))

        assert np.allclose(even_entries, even_refs, atol=tol, rtol=0) 

def test_ones(self, n):
        """Check all ones matrix has perm(J_n)=n!"""
        A = np.array([[1]])
        p = permanent_repeated(A, [n])
        assert np.allclose(p, fac(n)) 

def _naive_csh_ph(l, m, theta, phi):
    """
    CSH as defined by Pinchon-Hoggan. Same as wikipedia's quantum-normalized SH = naive_Y_quantum()
    """
    if l == 0 and m == 0:
        return 1. / np.sqrt(4 * np.pi)
    else:
        phase = ((1j) ** (m + np.abs(m)))
        normalizer = np.sqrt(((2 * l + 1.) * factorial(l - np.abs(m)))
                             /
                             (4 * np.pi * factorial(l + np.abs(m))))
        P = lpmv(np.abs(m), l, np.cos(theta))
        e = np.exp(1j * m * phi)
        return phase * normalizer * P * e 

def test_block_ones(self, n, dtype, recursive):
        """Check hafnian([[0, I_n], [I_n, 0]])=n!"""
        O = np.zeros([n, n])
        B = np.ones([n, n])
        A = np.vstack([np.hstack([O, B]), np.hstack([B, O])])
        A = dtype(A)
        haf = hafnian(A, recursive=recursive)
        expected = float(fac(n))
        assert np.allclose(haf, expected) 

def pdf(self, grid, dataSegment):
        """
        Probability density function of the Poisson model

        Args:
            grid(list): Parameter grid for discrete rate (lambda) values
            dataSegment(ndarray): Data segment from formatted data (in this case a single number of events)

        Returns:
            ndarray: Discretized Poisson pdf (with same shape as grid)
        """
        return (grid[0] ** dataSegment[0]) * (np.exp(-grid[0])) / (np.math.factorial(dataSegment[0])) 

def test_reference_squeezed_states(
        self, setup_backend, r, theta, batch_size, pure, cutoff, tol
    ):
        """Tests if a range of squeezed vacuum states are equal to the form of Eq. (5.5.6) in Loudon."""
        backend = setup_backend(1)

        backend.prepare_squeezed_state(r, theta, 0)
        s = backend.state()

        if s.is_pure:
            num_state = s.ket()
        else:
            num_state = s.dm()

        n = np.arange(0, cutoff, 2)
        even_refs = (
            np.sqrt(sech(r))
            * np.sqrt(factorial(n))
            / factorial(n / 2)
            * (-0.5 * np.exp(1j * theta) * np.tanh(r)) ** (n / 2)
        )

        if batch_size is not None:
            if pure:
                even_entries = num_state[:, ::2]
            else:
                even_entries = num_state[:, ::2, ::2]
                even_refs = np.outer(even_refs, np.conj(even_refs))
        else:
            if pure:
                even_entries = num_state[::2]
            else:
                even_entries = num_state[::2, ::2]
                even_refs = np.outer(even_refs, np.conj(even_refs))

        assert np.allclose(even_entries, even_refs, atol=tol, rtol=0.0) 

def test_beamsplitter_on_mode_subset(
            self, setup_backend, mag_alpha, t, r_phi, modes, cutoff, tol
    ):
        """Tests applying the beamsplitter on different mode subsets."""

        phase_alpha = np.pi / 5
        alpha = mag_alpha * np.exp(1j * phase_alpha)
        r = np.exp(1j * r_phi) * np.sqrt(1.0 - np.abs(t) ** 2)

        backend = setup_backend(4)

        backend.displacement(mag_alpha, phase_alpha, modes[0])
        backend.beamsplitter(np.arccos(t), r_phi, *modes)
        state = backend.state()

        alpha_outA = t * alpha
        alpha_outB = r * alpha

        n = np.arange(cutoff)
        ref_stateA = (
            np.exp(-0.5 * np.abs(alpha_outA) ** 2)
            * alpha_outA ** n
            / np.sqrt(factorial(n))
        )
        ref_stateB = (
            np.exp(-0.5 * np.abs(alpha_outB) ** 2)
            * alpha_outB ** n
            / np.sqrt(factorial(n))
        )

        numer_state = state.reduced_dm(list(modes))

        ref_state = np.einsum(
            "i,j,k,l->ijkl",
            ref_stateA,
            np.conj(ref_stateA),
            ref_stateB,
            np.conj(ref_stateB),
        )

        assert np.allclose(numer_state, ref_state, atol=tol, rtol=0) 

def wigner_d_naive(l, m, n, beta):
    """
    Numerically naive implementation of the Wigner-d function.
    This is useful for checking the correctness of other implementations.

    :param l: the degree of the Wigner-d function. l >= 0
    :param m: the order of the Wigner-d function. -l <= m <= l
    :param n: the order of the Wigner-d function. -l <= n <= l
    :param beta: the argument. 0 <= beta <= pi
    :return: d^l_mn(beta) in the TODO: what basis? complex, quantum(?), centered, cs(?)
    """
    from scipy.special import eval_jacobi
    try:
        from scipy.misc import factorial
    except:
        from scipy.special import factorial

    from sympy.functions.special.polynomials import jacobi, jacobi_normalized
    from sympy.abc import j, a, b, x
    from sympy import N
    #jfun = jacobi_normalized(j, a, b, x)
    jfun = jacobi(j, a, b, x)
    # eval_jacobi = lambda q, r, p, o: float(jfun.eval(int(q), int(r), int(p), float(o)))
    # eval_jacobi = lambda q, r, p, o: float(N(jfun, int(q), int(r), int(p), float(o)))
    eval_jacobi = lambda q, r, p, o: float(jfun.subs({j:int(q), a:int(r), b:int(p), x:float(o)}))

    mu = np.abs(m - n)
    nu = np.abs(m + n)
    s = l - (mu + nu) / 2
    xi = 1 if n >= m else (-1) ** (n - m)

    # print(s, mu, nu, np.cos(beta), type(s), type(mu), type(nu), type(np.cos(beta)))
    jac = eval_jacobi(s, mu, nu, np.cos(beta))
    z = np.sqrt((factorial(s) * factorial(s + mu + nu)) / (factorial(s + mu) * factorial(s + nu)))

    # print(l, m, n, beta, np.isfinite(mu), np.isfinite(nu), np.isfinite(s), np.isfinite(xi), np.isfinite(jac), np.isfinite(z))
    assert np.isfinite(mu) and np.isfinite(nu) and np.isfinite(s) and np.isfinite(xi) and np.isfinite(jac) and np.isfinite(z)
    assert np.isfinite(xi * z * np.sin(beta / 2) ** mu * np.cos(beta / 2) ** nu * jac)
    return xi * z * np.sin(beta / 2) ** mu * np.cos(beta / 2) ** nu * jac