Prime numbers and divisors

I often use prime numbers and divisors in python, but I wrote it because there was no library that handles prime numbers and divisors.

Primality test and prime number enumeration

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.

Primality test

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

Prime number enumeration

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

Prime factorization and divisors

You can find the number of divisors using prime factorization.

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

Number of divisors

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

Divisor enumeration

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

Number of divisors of factorial

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

Prime numbers and divisors
Discrimination of prime numbers
Recursively find prime numbers
Prime numbers in Python
This year's prime numbers
Determine prime numbers with python
Fibonacci and prime implementations (python)
Sorting with mixed numbers and letters
Project Euler 10 "Sum of Prime Numbers"
It's a prime number ... Counting prime numbers ...
[Python] nCr mod Compute prime numbers
Aggregation and visualization of accumulated numbers
Algorithm learned with Python 4th: Prime numbers
I searched for prime numbers in python
Solve with Ruby and Python AtCoder ABC084 D Cumulative sum of prime numbers