Python scipy.special.factorial2() Examples

def error_for_ke_cutoff(cell, ke_cutoff):
    kmax = np.sqrt(ke_cutoff*2)
    errmax = 0
    for i in range(cell.nbas):
        l = cell.bas_angular(i)
        es = cell.bas_exp(i)
        cs = abs(cell.bas_ctr_coeff(i)).max(axis=1)
        fac = (256*np.pi**4*cs**4 * factorial2(l*4+3)
               / factorial2(l*2+1)**2)
        efac = np.exp(-ke_cutoff/(2*es))
        err1 = .5*fac/(4*es)**(2*l+1) * kmax**(4*l+3) * efac
        errmax = max(errmax, err1.max())
        if np.any(ke_cutoff < 5*es):
            err2 = (1.41*efac+2.51)*fac/(4*es)**(2*l+2) * kmax**(4*l+5)
            errmax = max(errmax, err2[ke_cutoff<5*es].max())
        if np.any(ke_cutoff < es):
            err2 = (1.41*efac+2.51)*fac/2**(2*l+2) * np.sqrt(2*es)
            errmax = max(errmax, err2[ke_cutoff<es].max())
    return errmax 

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 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 _estimate_ke_cutoff(alpha, l, c, precision=INTEGRAL_PRECISION, weight=1.):
    '''Energy cutoff estimation based on cubic lattice'''
    # This function estimates the energy cutoff for (ii|ii) type of electron
    # repulsion integrals. The energy cutoff for nuclear attraction is larger
    # than the energy cutoff for ERIs.  The estimated error is roughly
    #     error ~ 64 pi^3 c^2 /((2l+1)!!(4a)^l) (2Ecut)^{l+.5} e^{-Ecut/4a}
    # log_k0 = 3 + np.log(alpha) / 2
    # l2fac2 = factorial2(l*2+1)
    # log_rest = np.log(precision*l2fac2*(4*alpha)**l / (16*np.pi**2*c**2))
    # Enuc_cut = 4*alpha * (log_k0*(2*l+1) - log_rest)
    # Enuc_cut[Enuc_cut <= 0] = .5
    # log_k0 = .5 * np.log(Ecut*2)
    # Enuc_cut = 4*alpha * (log_k0*(2*l+1) - log_rest)
    # Enuc_cut[Enuc_cut <= 0] = .5
    #
    # However, nuclear attraction can be evaluated with the trick of Ewald
    # summation which largely reduces the requirements to the energy cutoff.
    # In practice, the cutoff estimation for ERIs as below should be enough.
    log_k0 = 3 + np.log(alpha) / 2
    l2fac2 = factorial2(l*2+1)
    log_rest = np.log(precision*l2fac2**2*(4*alpha)**(l*2+1) / (128*np.pi**4*c**4))
    Ecut = 2*alpha * (log_k0*(4*l+3) - log_rest)
    Ecut[Ecut <= 0] = .5
    log_k0 = .5 * np.log(Ecut*2)
    Ecut = 2*alpha * (log_k0*(4*l+3) - log_rest)
    Ecut[Ecut <= 0] = .5
    return Ecut 

def test_factorial2(self):
        assert_array_almost_equal([105., 384., 945.],
                                  special.factorial2([7, 8, 9], exact=False))
        assert_equal(special.factorial2(7, exact=True), 105) 

def low_rank_hafnian(G):
    r"""Returns the hafnian of the low rank matrix :math:`\bm{A} = \bm{G} \bm{G}^T` where :math:`\bm{G}` is rectangular of size
    :math:`n \times r`  with :math:`r \leq n`.

    Note that the rank of :math:`\bm{A}` is precisely :math:`r`.

    The hafnian is calculated using the algorithm described in Appendix C of
    *A faster hafnian formula for complex matrices and its benchmarking on a supercomputer*,
    :cite:`bjorklund2018faster`.

    Args:
        G (array): factorization of the low rank matrix A = G @ G.T.

    Returns:
        (complex): hafnian of A.
    """
    n, r = G.shape
    if n % 2 != 0:
        return 0
    if r == 1:
        return factorial2(n - 1) * np.prod(G)
    poly = 1
    x = symbols("x0:" + str(r))
    for k in range(n):
        term = 0
        for j in range(r):
            term += G[k, j] * x[j]
        poly = expand(poly * term)

    comb = partitions(r, n // 2)
    haf_val = 0.0
    for c in comb:
        monomial = 1
        facts = 1
        for i, pi in enumerate(c):
            monomial *= x[i] ** (2 * pi)
            facts = facts * factorial2(2 * pi - 1)
        haf_val += complex(poly.coeff(monomial) * facts)
    return haf_val 

def test_rank_one(n):
    """ Test the hafnian of rank one matrices so that it is within
    10% of the exact value """
    x = np.random.rand(n)
    A = np.outer(x, x)
    exact = factorial2(n - 1) * np.prod(x)
    approx = haf_real(A, approx=True, nsamples=10000)
    assert np.allclose(approx, exact, rtol=2e-1, atol=0) 

def test_pmp(self, n):
        r"""Checks that the number of elements in pmp(n) is precisely the (n-1)!! for even n"""
        length = len(list(pmp(tuple(range(n)))))
        assert np.allclose(length, factorial2(n - 1)) 

def test_hafnian(self, n):
        r"""Checks that the hafnian of the all ones matrix of size n is (n-1)!!"""
        M = np.ones([n, n])
        if n % 2 == 0:
            assert np.allclose(factorial2(n - 1), hafnian(M))
        else:
            assert np.allclose(0, hafnian(M)) 

def _mom(self, k):
        return .5*special.factorial2(k-1)*(1+(-1)**k) 

def normalize(self):
        ''' Routine to normalize the basis functions, in case they
            do not integrate to unity.
        '''
        l,m,n = self.shell
        L = l+m+n
        # self.norm is a list of length equal to number primitives
        # normalize primitives first (PGBFs)
        self.norm = np.sqrt(np.power(2,2*(l+m+n)+1.5)*
                        np.power(self.exps,l+m+n+1.5)/
                        fact2(2*l-1)/fact2(2*m-1)/
                        fact2(2*n-1)/np.power(np.pi,1.5))

        # now normalize the contracted basis functions (CGBFs)
        # Eq. 1.44 of Valeev integral whitepaper
        prefactor = np.power(np.pi,1.5)*\
            fact2(2*l - 1)*fact2(2*m - 1)*fact2(2*n - 1)/np.power(2.0,L)

        N = 0.0
        num_exps = len(self.exps)
        for ia in range(num_exps):
            for ib in range(num_exps):
                N += self.norm[ia]*self.norm[ib]*self.coefs[ia]*self.coefs[ib]/\
                         np.power(self.exps[ia] + self.exps[ib],L+1.5)

        N *= prefactor
        N = np.power(N,-0.5)
        for ia in range(num_exps):
            self.coefs[ia] *= N