def __init__(self, n, s):
self.name = 'GrundmannMöller(dim={}, {})'.format(n, s)
d = 2*s + 1
self.degree = d
self.dim = n
exponents = get_all_exponents(n+1, s)
data = [
(
fr((-1)**i * 2**(-2*s) * (d+n-2*i)**d, fact(i) * fact(d+n-i)),
numpy.array([
[fr(2*p + 1, d+n-2*i) for p in part]
for part in exponents[s-i]
])
)
for i in range(s+1)
]
self.bary, self.weights = untangle(data)
self.weights /= sum(self.weights)
self.points = self.bary[:, 1:]
return
python类factorial()的实例源码
def __init__(self, settings):
Trajectory.__init__(self, settings)
# setup symbolic expressions
tau, k = sp.symbols('tau, k')
gamma = self._settings["differential_order"] + 1
alpha = sp.factorial(2 * gamma + 1)
f = sp.binomial(gamma, k) * (-1) ** k * tau ** (gamma + k + 1) / (gamma + k + 1)
phi = alpha / sp.factorial(gamma) ** 2 * sp.summation(f, (k, 0, gamma))
# differentiate phi(tau), index in list corresponds to order
dphi_sym = [phi] # init with phi(tau)
for order in range(self._settings["differential_order"]):
dphi_sym.append(dphi_sym[-1].diff(tau))
# lambdify
self.dphi_num = []
for der in dphi_sym:
self.dphi_num.append(sp.lambdify(tau, der, 'numpy'))
def test_gosper_sum_AeqB_part2():
f2a = n**2*a**n
f2b = (n - r/2)*binomial(r, n)
f2c = factorial(n - 1)**2/(factorial(n - x)*factorial(n + x))
g2a = -a*(a + 1)/(a - 1)**3 + a**(
m + 1)*(a**2*m**2 - 2*a*m**2 + m**2 - 2*a*m + 2*m + a + 1)/(a - 1)**3
g2b = (m - r)*binomial(r, m)/2
ff = factorial(1 - x)*factorial(1 + x)
g2c = 1/ff*(
1 - 1/x**2) + factorial(m)**2/(x**2*factorial(m - x)*factorial(m + x))
g = gosper_sum(f2a, (n, 0, m))
assert g is not None and simplify(g - g2a) == 0
g = gosper_sum(f2b, (n, 0, m))
assert g is not None and simplify(g - g2b) == 0
g = gosper_sum(f2c, (n, 1, m))
assert g is not None and simplify(g - g2c) == 0
def eval(cls, nu, z):
from sympy import (unpolarify, expand_mul, uppergamma, exp, gamma,
factorial)
nu2 = unpolarify(nu)
if nu != nu2:
return expint(nu2, z)
if nu.is_Integer and nu <= 0 or (not nu.is_Integer and (2*nu).is_Integer):
return unpolarify(expand_mul(z**(nu - 1)*uppergamma(1 - nu, z)))
# Extract branching information. This can be deduced from what is
# explained in lowergamma.eval().
z, n = z.extract_branch_factor()
if n == 0:
return
if nu.is_integer:
if (nu > 0) != True:
return
return expint(nu, z) \
- 2*pi*I*n*(-1)**(nu - 1)/factorial(nu - 1)*unpolarify(z)**(nu - 1)
else:
return (exp(2*I*pi*nu*n) - 1)*z**(nu - 1)*gamma(1 - nu) + expint(nu, z)
def _eval_aseries(self, n, args0, x, logx):
point = args0[0]
# Expansion at oo
if point is S.Infinity:
z = self.args[0]
# expansion of S(x) = S1(x*sqrt(pi/2)), see reference[5] page 1-8
p = [(-1)**k * C.factorial(4*k + 1) /
(2**(2*k + 2) * z**(4*k + 3) * 2**(2*k)*C.factorial(2*k))
for k in xrange(0, n)]
q = [1/(2*z)] + [(-1)**k * C.factorial(4*k - 1) /
(2**(2*k + 1) * z**(4*k + 1) * 2**(2*k - 1)*C.factorial(2*k - 1))
for k in xrange(1, n)]
p = [-sqrt(2/pi)*t for t in p] + [C.Order(1/z**n, x)]
q = [-sqrt(2/pi)*t for t in q] + [C.Order(1/z**n, x)]
return S.Half + (C.sin(z**2)*Add(*p) + C.cos(z**2)*Add(*q)).subs(x, sqrt(2/pi)*x)
# All other points are not handled
return super(fresnels, self)._eval_aseries(n, args0, x, logx)
def _eval_aseries(self, n, args0, x, logx):
point = args0[0]
# Expansion at oo
if point is S.Infinity:
z = self.args[0]
l = [ 1/sqrt(S.Pi) * C.factorial(2*k)*(-S(
4))**(-k)/C.factorial(k) * (1/z)**(2*k + 1) for k in xrange(0, n) ]
o = C.Order(1/z**(2*n + 1), x)
# It is very inefficient to first add the order and then do the nseries
return (Add(*l))._eval_nseries(x, n, logx) + o
# Expansion at I*oo
t = point.extract_multiplicatively(S.ImaginaryUnit)
if t is S.Infinity:
z = self.args[0]
# TODO: is the series really correct?
l = [ 1/sqrt(S.Pi) * C.factorial(2*k)*(-S(
4))**(-k)/C.factorial(k) * (1/z)**(2*k + 1) for k in xrange(0, n) ]
o = C.Order(1/z**(2*n + 1), x)
# It is very inefficient to first add the order and then do the nseries
return (Add(*l))._eval_nseries(x, n, logx) + o
# All other points are not handled
return super(_erfs, self)._eval_aseries(n, args0, x, logx)
def test_maxima_functions():
assert parse_maxima('expand( (x+1)^2)') == x**2 + 2*x + 1
assert parse_maxima('factor( x**2 + 2*x + 1)') == (x + 1)**2
assert parse_maxima('2*cos(x)^2 + sin(x)^2') == 2*cos(x)**2 + sin(x)**2
assert parse_maxima('trigexpand(sin(2*x)+cos(2*x))') == \
-1 + 2*cos(x)**2 + 2*cos(x)*sin(x)
assert parse_maxima('solve(x^2-4,x)') == [-2, 2]
assert parse_maxima('limit((1+1/x)^x,x,inf)') == E
assert parse_maxima('limit(sqrt(-x)/x,x,0,minus)') == -oo
assert parse_maxima('diff(x^x, x)') == x**x*(1 + log(x))
assert parse_maxima('sum(k, k, 1, n)', name_dict=dict(
n=Symbol('n', integer=True),
k=Symbol('k', integer=True)
)) == (n**2 + n)/2
assert parse_maxima('product(k, k, 1, n)', name_dict=dict(
n=Symbol('n', integer=True),
k=Symbol('k', integer=True)
)) == factorial(n)
assert parse_maxima('ratsimp((x^2-1)/(x+1))') == x - 1
assert Abs( parse_maxima(
'float(sec(%pi/3) + csc(%pi/3))') - 3.154700538379252) < 10**(-5)
def power_series(z, t, C, spatial_der_order=0):
if not all([isinstance(item, (Number, np.ndarray)) for item in [z, t]]):
raise TypeError
z = np.atleast_1d(z)
t = np.atleast_1d(t)
if not all([len(item.shape) == 1 for item in [z, t]]):
raise ValueError
x = np.nan*np.zeros((len(t), len(z)))
for i in range(len(z)):
sum_x = np.zeros(t.shape[0])
for j in range(len(C)-spatial_der_order):
sum_x += C[j+spatial_der_order][0, :]*z[i]**j/sm.factorial(j)
x[:, i] = sum_x
if any([dim == 1 for dim in x.shape]):
x = x.flatten()
return x
def test_limit_seq_fail():
# improve Summation algorithm or add ad-hoc criteria
e = (harmonic(n)**3 * Sum(1/harmonic(k), (k, 1, n)) /
(n * Sum(harmonic(k)/k, (k, 1, n))))
assert limit_seq(e, n) == 2
# No unique dominant term
e = (Sum(2**k * binomial(2*k, k) / k**2, (k, 1, n)) /
(Sum(2**k/k*2, (k, 1, n)) * Sum(binomial(2*k, k), (k, 1, n))))
assert limit_seq(e, n) == S(3) / 7
# Simplifications of summations needs to be improved.
e = n**3*Sum(2**k/k**2, (k, 1, n))**2 / (2**n * Sum(2**k/k, (k, 1, n)))
assert limit_seq(e, n) == 2
e = (harmonic(n) * Sum(2**k/k, (k, 1, n)) /
(n * Sum(2**k*harmonic(k)/k**2, (k, 1, n))))
assert limit_seq(e, n) == 1
e = (Sum(2**k*factorial(k) / k**2, (k, 1, 2*n)) /
(Sum(4**k/k**2, (k, 1, n)) * Sum(factorial(k), (k, 1, 2*n))))
assert limit_seq(e, n) == S(3) / 16
def eval(cls, nu, z):
from sympy import (unpolarify, expand_mul, uppergamma, exp, gamma,
factorial)
nu2 = unpolarify(nu)
if nu != nu2:
return expint(nu2, z)
if nu.is_Integer and nu <= 0 or (not nu.is_Integer and (2*nu).is_Integer):
return unpolarify(expand_mul(z**(nu - 1)*uppergamma(1 - nu, z)))
# Extract branching information. This can be deduced from what is
# explained in lowergamma.eval().
z, n = z.extract_branch_factor()
if n == 0:
return
if nu.is_integer:
if (nu > 0) != True:
return
return expint(nu, z) \
- 2*pi*I*n*(-1)**(nu - 1)/factorial(nu - 1)*unpolarify(z)**(nu - 1)
else:
return (exp(2*I*pi*nu*n) - 1)*z**(nu - 1)*gamma(1 - nu) + expint(nu, z)
def _eval_aseries(self, n, args0, x, logx):
from sympy import Order
point = args0[0]
# Expansion at oo
if point is S.Infinity:
z = self.args[0]
# expansion of S(x) = S1(x*sqrt(pi/2)), see reference[5] page 1-8
p = [(-1)**k * factorial(4*k + 1) /
(2**(2*k + 2) * z**(4*k + 3) * 2**(2*k)*factorial(2*k))
for k in range(0, n)]
q = [1/(2*z)] + [(-1)**k * factorial(4*k - 1) /
(2**(2*k + 1) * z**(4*k + 1) * 2**(2*k - 1)*factorial(2*k - 1))
for k in range(1, n)]
p = [-sqrt(2/pi)*t for t in p] + [Order(1/z**n, x)]
q = [-sqrt(2/pi)*t for t in q] + [Order(1/z**n, x)]
return S.Half + (sin(z**2)*Add(*p) + cos(z**2)*Add(*q)).subs(x, sqrt(2/pi)*x)
# All other points are not handled
return super(fresnels, self)._eval_aseries(n, args0, x, logx)
def compute_dobrodeev(n, I0, I2, I22, I4, pm_type, i, j, k):
'''Compute some helper quantities used in
L.N. Dobrodeev,
Cubature rules with equal coefficients for integrating functions with
respect to symmetric domains,
USSR Computational Mathematics and Mathematical Physics,
Volume 18, Issue 4, 1978, Pages 27-34,
<https://doi.org/10.1016/0041-5553(78)90064-2>.
'''
t = 1 if pm_type == 'I' else -1
L = binomial(n, i) * 2**i
M = fact(n) // (fact(j) * fact(k) * fact(n-j-k)) * 2**(j+k)
N = L + M
F = I22/I0 - I2**2/I0**2 + (I4/I0 - I22/I0) / n
R = -(j+k-i) / i * I2**2/I0**2 + (j+k-1)/n * I4/I0 - (n-1)/n * I22/I0
H = 1/i * (
(j+k-i) * I2**2/I0**2 + (j+k)/n * ((i-1) * I4/I0 - (n-1)*I22/I0)
)
Q = L/M*R + H - t * 2*I2/I0 * (j+k-i)/i * sqrt(L/M*F)
G = 1/N
a = sqrt(n/i * (I2/I0 + t * sqrt(M/L*F)))
b = sqrt(n/(j+k) * (I2/I0 - t * sqrt(L/M*F) + t * sqrt(k/j*Q)))
c = sqrt(n/(j+k) * (I2/I0 - t * sqrt(L/M*F) - t * sqrt(j/k*Q)))
return G, a, b, c
def _generate_i(n, i):
L = fact(n) // fact(i) // fact(n-i) * 2**i
G = fr(1, L)
a = sqrt(fr(3, n+2))
return G, a
def _generate_jk(n, pm_type, j, k):
M = fact(n) // fact(j) // fact(k) // fact(n-j-k) * 2**(j+k)
G = fr(1, M)
t = 1 if pm_type == 'I' else -1
b = sqrt(fr(1, j+k) * (1 + t * (k/j * sqrt(3*(j+k)/(n+2) - 1))))
c = sqrt(fr(1, j+k) * (1 - t * (j/k * sqrt(3*(j+k)/(n+2) - 1))))
return G, b, c
# pylint: disable=too-many-arguments, too-many-locals
def pdf(self, k):
return self.lamda**k / factorial(k) * exp(-self.lamda)
def test_expressions_failing():
assert residue(1/(x**4 + 1), x, exp(I*pi/4)) == -(S(1)/4 + I/4)/sqrt(2)
n = Symbol('n', integer=True, positive=True)
assert residue(exp(z)/(z - pi*I/4*a)**n, z, I*pi*a) == \
exp(I*pi*a/4)/factorial(n - 1)
assert residue(1/(x**2 + a**2)**2, x, a*I) == -I/4/a**3
def test_factorial():
from sympy import factorial, E
f = factorial(x)
assert limit(f, x, oo) == oo
assert limit(x/f, x, oo) == 0
# see Stirling's approximation:
# http://en.wikipedia.org/wiki/Stirling's_approximation
assert limit(f/(sqrt(2*pi*x)*(x/E)**x), x, oo) == 1
assert limit(f, x, -oo) == factorial(-oo)
assert limit(f, x, x**2) == factorial(x**2)
assert limit(f, x, -x**2) == factorial(-x**2)
def test_gosper_sum():
assert gosper_sum(1, (k, 0, n)) == 1 + n
assert gosper_sum(k, (k, 0, n)) == n*(1 + n)/2
assert gosper_sum(k**2, (k, 0, n)) == n*(1 + n)*(1 + 2*n)/6
assert gosper_sum(k**3, (k, 0, n)) == n**2*(1 + n)**2/4
assert gosper_sum(2**k, (k, 0, n)) == 2*2**n - 1
assert gosper_sum(factorial(k), (k, 0, n)) is None
assert gosper_sum(binomial(n, k), (k, 0, n)) is None
assert gosper_sum(factorial(k)/k**2, (k, 0, n)) is None
assert gosper_sum((k - 3)*factorial(k), (k, 0, n)) is None
assert gosper_sum(k*factorial(k), k) == factorial(k)
assert gosper_sum(
k*factorial(k), (k, 0, n)) == n*factorial(n) + factorial(n) - 1
assert gosper_sum((-1)**k*binomial(n, k), (k, 0, n)) == 0
assert gosper_sum((
-1)**k*binomial(n, k), (k, 0, m)) == -(-1)**m*(m - n)*binomial(n, m)/n
assert gosper_sum((4*k + 1)*factorial(k)/factorial(2*k + 1), (k, 0, n)) == \
(2*factorial(2*n + 1) - factorial(n))/factorial(2*n + 1)
# issue 6033:
assert gosper_sum(
n*(n + a + b)*a**n*b**n/(factorial(n + a)*factorial(n + b)), \
(n, 0, m)) == -a*b*(exp(m*log(a))*exp(m*log(b))*factorial(a)* \
factorial(b) - factorial(a + m)*factorial(b + m))/(factorial(a)* \
factorial(b)*factorial(a + m)*factorial(b + m))
def test_gosper_sum_AeqB_part1():
f1a = n**4
f1b = n**3*2**n
f1c = 1/(n**2 + sqrt(5)*n - 1)
f1d = n**4*4**n/binomial(2*n, n)
f1e = factorial(3*n)/(factorial(n)*factorial(n + 1)*factorial(n + 2)*27**n)
f1f = binomial(2*n, n)**2/((n + 1)*4**(2*n))
f1g = (4*n - 1)*binomial(2*n, n)**2/((2*n - 1)**2*4**(2*n))
f1h = n*factorial(n - S(1)/2)**2/factorial(n + 1)**2
g1a = m*(m + 1)*(2*m + 1)*(3*m**2 + 3*m - 1)/30
g1b = 26 + 2**(m + 1)*(m**3 - 3*m**2 + 9*m - 13)
g1c = (m + 1)*(m*(m**2 - 7*m + 3)*sqrt(5) - (
3*m**3 - 7*m**2 + 19*m - 6))/(2*m**3*sqrt(5) + m**4 + 5*m**2 - 1)/6
g1d = -S(2)/231 + 2*4**m*(m + 1)*(63*m**4 + 112*m**3 + 18*m**2 -
22*m + 3)/(693*binomial(2*m, m))
g1e = -S(9)/2 + (81*m**2 + 261*m + 200)*factorial(
3*m + 2)/(40*27**m*factorial(m)*factorial(m + 1)*factorial(m + 2))
g1f = (2*m + 1)**2*binomial(2*m, m)**2/(4**(2*m)*(m + 1))
g1g = -binomial(2*m, m)**2/4**(2*m)
g1h = 4*pi -(2*m + 1)**2*(3*m + 4)*factorial(m - S(1)/2)**2/factorial(m + 1)**2
g = gosper_sum(f1a, (n, 0, m))
assert g is not None and simplify(g - g1a) == 0
g = gosper_sum(f1b, (n, 0, m))
assert g is not None and simplify(g - g1b) == 0
g = gosper_sum(f1c, (n, 0, m))
assert g is not None and simplify(g - g1c) == 0
g = gosper_sum(f1d, (n, 0, m))
assert g is not None and simplify(g - g1d) == 0
g = gosper_sum(f1e, (n, 0, m))
assert g is not None and simplify(g - g1e) == 0
g = gosper_sum(f1f, (n, 0, m))
assert g is not None and simplify(g - g1f) == 0
g = gosper_sum(f1g, (n, 0, m))
assert g is not None and simplify(g - g1g) == 0
g = gosper_sum(f1h, (n, 0, m))
# need to call rewrite(gamma) here because we have terms involving
# factorial(1/2)
assert g is not None and simplify(g - g1h).rewrite(gamma) == 0
def test_gosper_nan():
a = Symbol('a', positive=True)
b = Symbol('b', positive=True)
n = Symbol('n', integer=True)
m = Symbol('m', integer=True)
f2d = n*(n + a + b)*a**n*b**n/(factorial(n + a)*factorial(n + b))
g2d = 1/(factorial(a - 1)*factorial(
b - 1)) - a**(m + 1)*b**(m + 1)/(factorial(a + m)*factorial(b + m))
g = gosper_sum(f2d, (n, 0, m))
assert simplify(g - g2d) == 0
def test_rsolve_hyper():
assert rsolve_hyper([-1, -1, 1], 0, n) in [
C0*(S.Half - S.Half*sqrt(5))**n + C1*(S.Half + S.Half*sqrt(5))**n,
C1*(S.Half - S.Half*sqrt(5))**n + C0*(S.Half + S.Half*sqrt(5))**n,
]
assert rsolve_hyper([n**2 - 2, -2*n - 1, 1], 0, n) in [
C0*rf(sqrt(2), n) + C1*rf(-sqrt(2), n),
C1*rf(sqrt(2), n) + C0*rf(-sqrt(2), n),
]
assert rsolve_hyper([n**2 - k, -2*n - 1, 1], 0, n) in [
C0*rf(sqrt(k), n) + C1*rf(-sqrt(k), n),
C1*rf(sqrt(k), n) + C0*rf(-sqrt(k), n),
]
assert rsolve_hyper(
[2*n*(n + 1), -n**2 - 3*n + 2, n - 1], 0, n) == C1*factorial(n) + C0*2**n
assert rsolve_hyper(
[n + 2, -(2*n + 3)*(17*n**2 + 51*n + 39), n + 1], 0, n) == None
assert rsolve_hyper([-n - 1, -1, 1], 0, n) == None
assert rsolve_hyper([-1, 1], n, n).expand() == C0 + n**2/2 - n/2
assert rsolve_hyper([-1, 1], 1 + n, n).expand() == C0 + n**2/2 + n/2
assert rsolve_hyper([-1, 1], 3*(n + n**2), n).expand() == C0 + n**3 - n
assert rsolve_hyper([-a, 1],0,n).expand() == C0*a**n
assert rsolve_hyper([-a, 0, 1], 0, n).expand() == (-1)**n*C1*a**(n/2) + C0*a**(n/2)
assert rsolve_hyper([1, 1, 1], 0, n).expand() == \
C0*(-S(1)/2 - sqrt(3)*I/2)**n + C1*(-S(1)/2 + sqrt(3)*I/2)**n
assert rsolve_hyper([1, -2*n/a - 2/a, 1], 0, n) == None
def test_hyper_rewrite_sum():
from sympy import RisingFactorial, factorial, Dummy, Sum
_k = Dummy("k")
assert replace_dummy(hyper((1, 2), (1, 3), x).rewrite(Sum), _k) == \
Sum(x**_k / factorial(_k) * RisingFactorial(2, _k) /
RisingFactorial(3, _k), (_k, 0, oo))
assert hyper((1, 2, 3), (-1, 3), z).rewrite(Sum) == \
hyper((1, 2, 3), (-1, 3), z)
def taylor_term(n, x, *previous_terms):
if n < 0 or n % 2 == 0:
return S.Zero
else:
x = sympify(x)
k = C.floor((n - 1)/S(2))
if len(previous_terms) > 2:
return -previous_terms[-2] * x**2 * (n - 2)/(n*k)
else:
return 2*(-1)**k * x**n/(n*C.factorial(k)*sqrt(S.Pi))
def taylor_term(n, x, *previous_terms):
if n == 0:
return S.One
elif n < 0 or n % 2 == 0:
return S.Zero
else:
x = sympify(x)
k = C.floor((n - 1)/S(2))
if len(previous_terms) > 2:
return -previous_terms[-2] * x**2 * (n - 2)/(n*k)
else:
return -2*(-1)**k * x**n/(n*C.factorial(k)*sqrt(S.Pi))
def taylor_term(n, x, *previous_terms):
if n < 0 or n % 2 == 0:
return S.Zero
else:
x = sympify(x)
k = C.floor((n - 1)/S(2))
if len(previous_terms) > 2:
return previous_terms[-2] * x**2 * (n - 2)/(n*k)
else:
return 2 * x**n/(n*C.factorial(k)*sqrt(S.Pi))
def taylor_term(n, x, *previous_terms):
if n < 0:
return S.Zero
else:
x = sympify(x)
if len(previous_terms) > 1:
p = previous_terms[-1]
return (-pi**2*x**4*(4*n - 1)/(8*n*(2*n + 1)*(4*n + 3))) * p
else:
return x**3 * (-x**4)**n * (S(2)**(-2*n - 1)*pi**(2*n + 1)) / ((4*n + 3)*C.factorial(2*n + 1))
def taylor_term(n, x, *previous_terms):
if n < 0:
return S.Zero
else:
x = sympify(x)
if len(previous_terms) > 1:
p = previous_terms[-1]
return (-pi**2*x**4*(4*n - 3)/(8*n*(2*n - 1)*(4*n + 1))) * p
else:
return x * (-x**4)**n * (S(2)**(-2*n)*pi**(2*n)) / ((4*n + 1)*C.factorial(2*n))
def _eval_aseries(self, n, args0, x, logx):
point = args0[0]
# Expansion at oo
if point is S.Infinity:
z = self.args[0]
# expansion of C(x) = C1(x*sqrt(pi/2)), see reference[5] page 1-8
p = [(-1)**k * C.factorial(4*k + 1) /
(2**(2*k + 2) * z**(4*k + 3) * 2**(2*k)*C.factorial(2*k))
for k in xrange(0, n)]
q = [1/(2*z)] + [(-1)**k * C.factorial(4*k - 1) /
(2**(2*k + 1) * z**(4*k + 1) * 2**(2*k - 1)*C.factorial(2*k - 1))
for k in xrange(1, n)]
p = [-sqrt(2/pi)*t for t in p] + [C.Order(1/z**n, x)]
q = [ sqrt(2/pi)*t for t in q] + [C.Order(1/z**n, x)]
return S.Half + (C.cos(z**2)*Add(*p) + C.sin(z**2)*Add(*q)).subs(x, sqrt(2/pi)*x)
# All other points are not handled
return super(fresnelc, self)._eval_aseries(n, args0, x, logx)
###############################################################################
#################### HELPER FUNCTIONS #########################################
###############################################################################
def eval_levicivita(*args):
"""Evaluate Levi-Civita symbol."""
from sympy import factorial
n = len(args)
return prod(
prod(args[j] - args[i] for j in xrange(i + 1, n))
/ factorial(i) for i in xrange(n))
# converting factorial(i) to int is slightly faster
def eval(cls, arg):
if arg.is_Number:
if arg is S.NaN:
return S.NaN
elif arg is S.Infinity:
return S.Infinity
elif arg.is_Integer:
if arg.is_positive:
return C.factorial(arg - 1)
else:
return S.ComplexInfinity
elif arg.is_Rational:
if arg.q == 2:
n = abs(arg.p) // arg.q
if arg.is_positive:
k, coeff = n, S.One
else:
n = k = n + 1
if n & 1 == 0:
coeff = S.One
else:
coeff = S.NegativeOne
for i in range(3, 2*k, 2):
coeff *= i
if arg.is_positive:
return coeff*sqrt(S.Pi) / 2**n
else:
return 2**n*sqrt(S.Pi) / coeff
