def multifactorial(self, int k):
r"""
Computes the k-th factorial `n!^{(k)}` of self. For k=1
this is the standard factorial, and for k greater than one it is
the product of every k-th terms down from self to k. The recursive
definition is used to extend this function to the negative
integers.
EXAMPLES::
sage: 5.multifactorial(1)
120
sage: 5.multifactorial(2)
15
sage: 23.multifactorial(2)
316234143225
sage: prod([1..23, step=2])
316234143225
sage: (-29).multifactorial(7)
1/2640
"""
if k <= 0:
raise ValueError, "multifactorial only defined for positive
values of k"
if not mpz_fits_sint_p(self.value):
raise ValueError, "multifactorial not implemented for n >=
2^32.\nThis is probably OK, since the answer would have billions of digits."
cdef int n = mpz_get_si(self.value)
# base case
if 0 < n < k:
return one
# easy to calculate
elif n % k == 0:
factorial = Integer(n/k).factorial()
if k == 2:
return factorial << (n/k)
else:
return factorial * Integer(k)**(n/k)
# negative base case
elif -k < n < 0:
return one / (self+k)
# reflection case
elif n < -k:
if (n/k) % 2:
sign = -one
else:
sign = one
return sign / Integer(-k-n).multifactorial(k)
# compute the actual product, optimizing the number of large
# multiplications
cdef int i,j
# we need (at most) log_2(#factors) concurrent sub-products
cdef int prod_count = <int>ceil_c(log_c(n/k+1)/log_c(2))
cdef mpz_t* sub_prods = <mpz_t*>check_allocarray(prod_count,
sizeof(mpz_t))
for i from 0 <= i < prod_count:
mpz_init(sub_prods[i])
sig_on()
cdef residue = n % k
cdef int tip = 0
for i from 1 <= i <= n//k:
mpz_set_ui(sub_prods[tip], k*i + residue)
# for the i-th terms we use the bits of i to calculate how many
# times we need to multiply "up" the stack of sub-products
for j from 0 <= j < 32:
if i & (1 << j):
break
tip -= 1
mpz_mul(sub_prods[tip], sub_prods[tip], sub_prods[tip+1])
tip += 1
cdef int last = tip-1
for tip from last > tip >= 0:
mpz_mul(sub_prods[tip], sub_prods[tip], sub_prods[tip+1])
sig_off()
cdef Integer z = PY_NEW(Integer)
mpz_swap(z.value, sub_prods[0])
for i from 0 <= i < prod_count:
mpz_clear(sub_prods[i])
sage_free(sub_prods)
return z
On Wednesday, October 28, 2015 at 2:41:59 PM UTC+5:30, prateek sharma wrote:
>
> Hi,
> I am looking for multifactorial function in the source code but unable
> to.Help me out...
>
--
You received this message because you are subscribed to the Google Groups
"sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sage-devel+unsubscr...@googlegroups.com.
To post to this group, send email to sage-devel@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-devel.
For more options, visit https://groups.google.com/d/optout.
