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.