Download code

Jump to: navigation, search

Back to Factorials_with_prime_factorization_(Python)

Download for Windows: single file, zip

Download for UNIX: single file, zip, tar.gz, tar.bz2

prime_factorial.py

 1 # The authors of this work have released all rights to it and placed it
 2 # in the public domain under the Creative Commons CC0 1.0 waiver
 3 # (http://creativecommons.org/publicdomain/zero/1.0/).
 4 # 
 5 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 6 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 7 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 8 # IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 9 # CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10 # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
11 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12 # 
13 # Retrieved from: http://en.literateprograms.org/Factorials_with_prime_factorization_(Python)?oldid=16897
14 
15 def multiplicity(n, p):
16     """Return the power of the prime number p in the
17     factorization of n!"""
18     if p > n: return 0
19     if p > n//2: return 1
20     q, m = n, 0
21     while q >= p:
22         q //= p
23         m += q
24     return m
25 
26 def primes(n):
27     """Generate a list of the prime numbers [2, 3, ... m] where
28     m is the largest prime <= n."""
29     n = n + 1
30     sieve = range(n)
31     sieve[:2] = [0, 0]
32     for i in xrange(2, int(n**0.5)+1):
33         if sieve[i]:
34             for j in xrange(i**2, n, i):
35                 sieve[j] = 0
36     # Filter out the composites, which have been replaced by 0's
37     return [p for p in sieve if p]
38 def powproduct(ns):
39     """Compute the explicit value of a factored integer
40     given as a list of (base, exponent) pairs."""
41     if not ns:
42         return 1
43     units = 1
44     multi = []
45     for base, exp in ns:
46         if exp == 0:
47             continue
48         elif exp == 1:
49             units *= base
50         else:
51             if exp % 2:
52                 units *= base
53             multi.append((base, exp//2))
54     return units * powproduct(multi)**2
55 def factorial(n):
56     return powproduct((p, multiplicity(n, p)) for p in primes(n))


hijacker
hijacker
hijacker
hijacker