1 from sympy
.core
import *
2 from sympy
.core
.basic
import S
4 from sympy
.utilities
.memoization
import recurrence_memo
6 # Lanczos approximation for low-precision numerical factorial
7 _lanczos_coef
= [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
8 771.32342877765313, -176.61502916214059, 12.507343278686905,
9 -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]
12 from cmath
import pi
, sin
, log
, exp
14 return pi
*z
/ (sin(pi
*z
) * _lanczos(-z
))
18 x
+= _lanczos_coef
[i
]/(z
+i
)
19 logw
= 0.91893853320467267+(z
+0.5)*log(z
+7.5)+log(x
)-z
-7.5
22 class _Factorial(SingleValuedFunction
):
24 Factorials and multiple factorials
28 factorial(x) gives the factorial of x, defined as
29 x! = 1*2*3*...*x if x is a positive integer. If x is not a
30 positive integer, the value of x! is defined in terms of the
33 factorial(x, m) returns the m-th order multifactorial of x;
34 i.e., 1 * ... * (x-2*m) * (x-m) * x. In particular,
35 factorial(x, 2) returns the double factorial of x,
37 * x!! = 2*4*...*(x-2)*x if x is a positive even integer
38 * x!! = 1*3*...*(x-2)*x if x is a positive odd integer
40 The argument m must be a positive integer.
48 >>> factorial(Rational(5,2))
63 # generators for multifactorials
67 def _eval_apply(cls
, x
, m
=1):
70 if m
== 1 and x
.is_integer
:
71 # handle factorial poles, also for symbols with assumptions
74 if isinstance(x
, Integer
):
75 return Integer(cls
._fac
1(int(x
)))
77 # half-integer case of the ordinary factorial
78 if m
== 1 and isinstance(x
, Rational
) and x
.q
== 2:
81 return (-1)**(-n
+1) * pi
* x
/ factorial(-x
)
82 return Basic
.sqrt(pi
) * Rational(1, 2**n
) * factorial(2*n
-1, 2)
84 # multifactorials are only defined for integers
85 if (not isinstance(m
, Integer
) and m
> 0) or not \
86 isinstance(x
, Integer
):
95 # e.g., for double factorials, the start value is 1 for odd x
101 # we can define the m-multifactorials for negative numbers
102 # to satisfy f(x,m) = f(x+m,m) / (x+m)
106 return factorial(x
+m
, m
) / (x
+m
)
109 if not (m
, start_value
) in cls
._generators
:
110 @recurrence_memo([start_value
])
112 return (k
*m
+ start_value
) * prev
[-1]
113 cls
._generators
[m
, start_value
] = _f
115 return Integer(cls
._generators
[m
, start_value
]((x
-start_value
) // m
))
117 # This should give a series expansion around x = oo. Needs fixing
118 # def series(self, x, n):
119 # return sqrt(2*pi*x) * x**x * exp(-x) * (1 + O(1/x))
121 # Derivatives are given in terms of polygamma functions
122 # (XXX: only for order m = 1)
123 def fdiff(self
, argindex
=1):
125 from sympy
.functions
import gamma
, polygamma
126 return gamma(self
[0]+1)*polygamma(0,self
[0]+1)
128 raise ArgumentIndexError(self
, argindex
)
131 def _eval_apply_evalf(cls
, x
, m
=1):
132 """Return a low-precision numerical approximation."""
134 a
, b
= x
.as_real_imag()
135 y
= _lanczos(complex(a
, b
))
139 Real(y
.real
) + I
*Real(y
.imag
)
141 factorial
= _Factorial
143 class UnevaluatedFactorial(_Factorial
):
146 def _eval_apply(cls
, x
, m
=1):
149 unfac
= UnevaluatedFactorial
153 # factorial_simplify helpers; could use refactoring
155 def _isfactorial(expr
):
156 #return isinstance(expr, Apply) and isinstance(expr[0], Factorial)
157 return isinstance(expr
, _Factorial
)
159 def _collect_factors(expr
):
160 assert isinstance(expr
, Mul
)
165 if isinstance(x
, Mul
):
166 n
, d
, o
= _collect_factors(x
)
170 elif isinstance(x
, Pow
):
172 if _isfactorial(base
) and \
173 isinstance(exp
, Rational
) and exp
.is_integer
:
175 for i
in xrange(exp
.p
): numer_args
.append(base
[0])
177 for i
in xrange(-exp
.p
): denom_args
.append(base
[0])
180 elif _isfactorial(x
):
181 numer_args
.append(x
[0])
184 return numer_args
, denom_args
, other
187 def _simplify_quotient(na
, da
, other
):
190 for i
, y
in enumerate(na
):
191 for j
, x
in enumerate(da
):
192 #delta = simplify(y - x)
194 if isinstance(delta
, Rational
) and delta
.is_integer
:
195 candidates
.append((delta
, i
, j
))
197 # There may be multiple candidates. Choose the quotient pair
198 # that minimizes the work.
199 candidates
.sort(key
=lambda x
: abs(x
[0]))
200 delta
, i
, j
= candidates
[0]
205 for k
in xrange(1, int(delta
)+1):
208 for k
in xrange(1, -int(delta
)+1):
215 # handle x!*(x+1) and x!/x
216 def _simplify_recurrence(facs
, other
, reciprocal
=False):
217 # this should probably be rewritten more elegantly
221 while j
< len(other
):
226 if othr
- fact
== 1: facs
[i
] += 1; del other
[j
]; j
-= 1
227 elif -othr
- fact
== 1: facs
[i
] += 1; del other
[j
]; other
.append(-1); j
-= 1
228 elif 1/othr
- fact
== 0: facs
[i
] -= 1; del other
[j
]; j
-= 1
229 elif -1/othr
- fact
== 0: facs
[i
] -= 1; del other
[j
]; other
.append(-1); j
-= 1
233 def factorial_simplify(expr
):
235 This function takes an expression containing factorials
236 and removes as many of them as possible by combining
237 products and quotients of factorials into polynomials
238 or other simpler expressions.
240 TODO: handle reflection formula, duplication formula
245 if isinstance(expr
, Add
):
246 return Add(*(factorial_simplify(x
) for x
in expr
))
248 if isinstance(expr
, Pow
):
249 return Pow(factorial_simplify(expr
[0]), expr
[1])
251 if isinstance(expr
, Mul
):
252 na
, da
, other
= _collect_factors(expr
)
253 _simplify_quotient(na
, da
, other
)
254 _simplify_recurrence(na
, other
)
255 _simplify_recurrence(da
, other
, reciprocal
=True)
258 for n
in na
: result
*= factorial(n
)
259 for d
in da
: result
/= factorial(d
)
260 for o
in other
: result
*= o
263 expr
= expr
.subs(unfac
, factorial
)
267 class Rising_factorial(SingleValuedFunction
):
271 Calculate the rising factorial (x)^(n) = x(x+1)...(x+n-1)
272 as a quotient of factorials
276 >>> rising_factorial(3, 2)
283 def _eval_apply(cls
, x
, n
):
284 return factorial_simplify(unfac(x
+n
-1) / unfac(x
-1))
287 rising_factorial
= Rising_factorial
290 class Falling_factorial(SingleValuedFunction
):
294 Calculate the falling factorial (x)_(n) = x(x-1)...(x-n+1)
295 as a quotient of factorials
299 >>> falling_factorial(5, 3)
306 def _eval_apply(cls
, x
, n
):
307 return factorial_simplify(unfac(x
) / unfac(x
-n
))
310 falling_factorial
= Falling_factorial
313 class Binomial2(SingleValuedFunction
):
317 Calculate the binomial coefficient C(n,k) = n!/(k!(n-k)!)
321 When n and k are positive integers, the result is always
324 If k is a positive integer, the result is a polynomial in n
325 that is evaluated explicitly.
331 >>> # Building Pascal's triangle
332 >>> [binomial2(0,k) for k in range(1)]
334 >>> [binomial2(1,k) for k in range(2)]
336 >>> [binomial2(2,k) for k in range(3)]
338 >>> [binomial2(3,k) for k in range(4)]
340 >>> # n can be arbitrary if k is a positive integer
341 >>> binomial2(Rational(5,4), 3)
345 (1/6)*x*(1 - x)*(2 - x)
351 def _eval_apply(cls
, n
, k
):
353 # TODO: move these two cases to factorial_simplify as well
354 if n
== 0 and k
!= 0:
355 return Basic
.sin(pi
*k
)/(pi
*k
)
357 return factorial_simplify(unfac(n
) / unfac(k
) / unfac(n
-k
))
359 binomial2
= Binomial2