PrimeSwingFactorial

1    import bisect 
2     
3    def Primes(n) : 
4        primes = [2, 3] 
5        lim, tog = n // 3, False 
6        composite = [False for i in range(lim)] 
7     
8        d1 = 8; d2 = 8; p1 = 3; p2 = 7; s = 7; s2 = 3; m = -1 
9     
10       while s < lim :             # --  scan the sieve 
11           m += 1                  # --  if a prime is found 
12           if not composite[m] :   # --  cancel its multiples 
13               inc = p1 + p2 
14               for k in range(s,      lim, inc) : composite[k] = True 
15               for k in range(s + s2, lim, inc) : composite[k] = True 
16    
17               tog = not tog 
18               if tog: s += d2; d1 += 16; p1 += 2; p2 += 2; s2 = p2 
19               else:   s += d1; d2 +=  8; p1 += 2; p2 += 6; s2 = p1 
20    
21       k, p, tog = 0, 5, False 
22       while p <= n : 
23           if not composite[k] : primes.append(p) 
24           k += 1; 
25           tog = not tog 
26           p += 2 if tog else 4 
27    
28       return primes 
29    
30   def isqrt(x): 
31       ''' 
32       Writing your own square root function
33       ''' 
34       if x < 0: raise ValueError('square root not defined for negative numbers') 
35       n = int(x) 
36       if n == 0: return 0 
37       a, b = divmod(n.bit_length(), 2) 
38       x = 2**(a + b) 
39       while True: 
40           y = (x + n // x) // 2 
41           if y >= x: return x 
42           x = y 
43            
44   def product(s, n, m): 
45       if n > m: return 1 
46       if n == m: return s[n] 
47       k = (n + m) // 2 
48       return product(s, n, k) * product(s, k + 1, m) 
49    
50   def factorialPS(n): 
51    
52       small_swing = [1,1,1,3,3,15,5,35,35,315,63,693,231,3003,429,6435,6435, 
53           109395,12155,230945,46189,969969,88179,2028117,676039,16900975, 
54           1300075,35102025,5014575,145422675,9694845,300540195,300540195] 
55    
56       def swing(m, primes): 
57           if m < 33: return small_swing[m] 
58            
59           s = bisect.bisect_left(primes, 1 + isqrt(m)) 
60           d = bisect.bisect_left(primes, 1 + m // 3) 
61           e = bisect.bisect_left(primes, 1 + m // 2) 
62           g = bisect.bisect_left(primes, 1 + m) 
63            
64           factors = primes[e:g] 
65           factors += filter(lambda x: (m // x) & 1 == 1, primes[s:d]) 
66           for prime in primes[1:s]:   
67               p, q = 1, m 
68               while True: 
69                   q //= prime 
70                   if q == 0: break 
71                   if q & 1 == 1: 
72                       p *= prime 
73               if p > 1: factors.append(p) 
74    
75           return product(factors, 0, len(factors) - 1) 
76    
77       def odd_factorial(n, primes): 
78           if n < 2: return 1 
79           return (odd_factorial(n // 2, primes)**2) * swing(n, primes) 
80    
81       def eval(n): 
82           if n < 0: 
83               raise ValueError('factorial not defined for negative numbers') 
84    
85           if n == 0: return 1 
86           if n < 20: return product(range(2, n + 1), 0, n-2) 
87    
88           N, bits = n, n 
89           while N != 0: 
90               bits -= N & 1 
91               N >>= 1 
92    
93           primes = Primes(n) 
94           return odd_factorial(n, primes) * 2**bits 
95    
96       return eval(n) 
97    
98   factorialPS(100) 

---

Back to the Homepage of Factorial Algorithms.