I often use prime numbers and divisors in python, but I wrote it because there was no library that handles prime numbers and divisors.
There are "trial division" and "Eratosthenes sieve" as algorithms related to primality test and prime number enumeration that have been around for a long time.
Implementation considered from the definition of prime numbers. Returns True if it is not divisible by all integers greater than or equal to 2 and less than or equal to $ \ sqrt {num} $.
import math
def is_prime_1(num):
if num < 2:
return False
else:
num_sqrt = math.floor(math.sqrt(num))
for prime in range(2, num_sqrt + 1):
if num % prime == 0:
return False
break
return True
Implementation using trial split.
def is_prime_2(num):
if num < 2:
return False
if num == 2 or num == 3 or num == 5:
return True
if num % 2 == 0 or num % 3 == 0 or num % 5 == 0:
return False
#Pseudo-prime number(A number that is not divisible by 2, 3, or 5)Divide one after another
prime = 7
step = 4
num_sqrt = math.sqrt(num)
while prime <= num_sqrt:
if num % prime == 0:
return False
prime += step
step = 6 - step
return True
Comparing the execution speed, the trial split is certainly faster.
%timeit is_prime_1(random.randrange(10 ** 5))
%timeit is_prime_2(random.randrange(10 ** 5))
100000 loops, best of 3: 3.67 µs per loop
100000 loops, best of 3: 2.84 µs per loop
Implementation using the primality test implemented above.
def make_prime_list_1(num):
if num < 2:
return []
prime_list = []
for prime in range(2, num + 1):
if is_prime_2(prime):
prime_list.append(prime)
return prime_list
Mounting with an Eratosthenes sieve.
def make_prime_list_2(num):
if num < 2:
return []
#0 is not a prime number
prime_list = [i for i in range(num + 1)]
prime_list[1] = 0 #1 is not a prime number
num_sqrt = math.sqrt(num)
for prime in prime_list:
if prime == 0:
continue
if prime > num_sqrt:
break
for non_prime in range(2 * prime, num, prime):
prime_list[non_prime] = 0
return [prime for prime in prime_list if prime != 0]
Compare execution speeds. The Eratosthenes sieve is faster.
%timeit -n1000 make_prime_list_1(random.randrange(10 ** 5))
%timeit -n1000 make_prime_list_2(random.randrange(10 ** 5))
1000 loops, best of 3: 69.5 ms per loop
1000 loops, best of 3: 10.5 ms per loop
By the way, it is extremely slow to prepare a prime number table first and perform a primality test.
prime_list = make_prime_list_2(10 ** 5)
def is_prime_3(num):
return num in prime_list
%timeit is_prime_3(random.randrange(10 ** 5))
10000 loops, best of 3: 189 µs per loop
You can find the number of divisors using prime factorization.
Implementation using defaultdict.
from collection import defaultdict
def prime_factorization_1(num):
if num <= 1:
return False
else:
num_sqrt = math.floor(math.sqrt(num))
prime_list = make_prime_list_2(num_sqrt)
dict_counter = defaultdict(int)
for prime in prime_list:
while num % prime == 0:
dict_counter[prime] += 1
num //= prime
if num != 1:
dict_counter[num] += 1
dict_counter = dict(dict_counter)
return dict_counter
Implementation without defaultdict.
def prime_factorization_2(num):
if num <= 1:
return False
else:
num_sqrt = math.floor(math.sqrt(num))
prime_list = make_prime_list_2(num_sqrt)
dict_counter = {}
for prime in prime_list:
while num % prime == 0:
if prime in dict_counter:
dict_counter[prime] += 1
else:
dict_counter[prime] = 1
num //= prime
if num != 1:
if num in dict_counter:
dict_counter[num] += 1
else:
dict_counter[num] = 1
return dict_counter
The speed does not change with or without defaultdict.
%timeit prime_factorization_1(random.randrange(10 ** 5))
%timeit prime_factorization_2(random.randrange(10 ** 5))
10000 loops, best of 3: 33.6 µs per loop
10000 loops, best of 3: 33.2 µs per loop
The product of each (number of divisors + 1) of the prime factorization is the number of divisors.
Implementation without prime factorization.
def search_divisor_num_1(num):
if num < 0:
return None
elif num == 1:
return 1
else:
num_sqrt = math.floor(math.sqrt(num))
prime_list = make_prime_list_2(num_sqrt)
divisor_num = 1
for prime in prime_list:
count = 1
while num % prime == 0:
num //= prime
count += 1
divisor_num *= count
if num != 1:
divisor_num *= 2
return divisor_num
Implementation using the prime factorization implemented above.
def search_divisor_num_2(num):
if num < 0:
return None
elif num == 1:
return 1
else:
divisor_num = 1
dict_fact = prime_factorization_1(num)
for value in dict_fact.values():
divisor_num *= (value + 1)
return divisor_num
There is not much difference when comparing the speed.
%timeit search_divisor_num_1(random.randrange(10 ** 5))
%timeit search_divisor_num_2(random.randrange(10 ** 5))
10000 loops, best of 3: 33.5 µs per loop
10000 loops, best of 3: 34.1 µs per loop
To be honest, this is just divided by num / 2, so I think it's faster to do it easily.
def make_divisor_list(num):
if num < 1:
return []
elif num == 1:
return [1]
else:
divisor_list = []
divisor_list.append(1)
for i in range(2, num // 2 + 1):
if num % i == 0:
divisor_list.append(i)
divisor_list.append(num)
return divisor_list
%timeit make_divisor_list(random.randrange(10 ** 5))
1000 loops, best of 3: 1.94 ms per loop
When you factorial, the order becomes quite large and overflows. Implementation using memoization to avoid it.
def search_divisor_num_of_factorial_num(num):
if num <= 0:
return 0
elif num == 1 or num == 2:
return 1
else:
#Enter the prime numbers and their numbers
dict_counter = defaultdict(int)
#Dict for notes
dict_memo = defaultdict(list)
for a_num in range(2, num + 1):
num_sqrt = math.ceil(math.sqrt(a_num))
prime_list = make_prime_list_2(num_sqrt)
#Keep the key to put in the memo dict
now_num = a_num
for prime in prime_list:
while a_num % prime == 0:
#If it's in the memo, move everything from there and exit the loop when you're done
if a_num in dict_memo:
for memo in dict_memo[a_num]:
dict_counter[memo] += 1
dict_memo[now_num].append(memo)
a_num = 1
else:
dict_counter[prime] += 1
dict_memo[now_num].append(prime)
a_num //= prime
if a_num != 1:
dict_counter[a_num] += 1
dict_memo[now_num].append(a_num)
divisor_num = 1
dict_fact = dict(dict_counter)
for value in dict_fact.values():
divisor_num *= (value + 1)
return divisor_num
Comparison with the case where the factorial is thrown as it is to the function to find the number of divisors above.
math.factorial(10) #in the case of
3628800
%timeit search_divisor_num_1(math.factorial(10))
%timeit search_divisor_num_of_factorial_num(10)
1000 loops, best of 3: 331 µs per loop
10000 loops, best of 3: 28.2 µs per loop
math.factorial(50) #in the case of
30414093201713378043612608166064768844377641568960512000000000000
%timeit search_divisor_num_1(math.factorial(50))
%timeit search_divisor_num_of_factorial_num(50)
Overflow
1000 loops, best of 3: 177 µs per loop
I know that SciPy has a module that handles prime numbers, but it's not a standard library, and ... It would be nice to add a library that can handle prime numbers and divisors to the standard library.
Recommended Posts