monster_divisors.py

# Source - https://stackoverflow.com/a/70638665
# Posted by Jérôme Richard
# Retrieved 2026-04-30, License - CC BY-SA 4.0

import numba as nb

@nb.njit('List(int_)(int_)')
def get_prime_divisors(n):
    divisors = []
    while n % 2 == 0:
        divisors.append(2)
        n //= 2
    while n % 3 == 0:
        divisors.append(3)
        n //= 3
    i = 5
    while i*i <= n:
        for k in (i, i+2):
            while n % k == 0:
                divisors.append(k)
                n //= k
        i += 6
    if n > 1:
        divisors.append(n)
    return divisors

@nb.njit('List(int_)(int_)')
def get_divisors(n):
    divisors = []
    if n == 1:
        divisors.append(1)
    elif n > 1:
        prime_factors = get_prime_divisors(n)
        divisors = [1]
        last_prime = 0
        factor = 0
        slice_len = 0
        # Find all the products that are divisors of n
        for prime in prime_factors:
            if last_prime != prime:
                slice_len = len(divisors)
                factor = prime
            else:
                factor *= prime
            for i in range(slice_len):
                divisors.append(divisors[i] * factor)
            last_prime = prime
        divisors.sort()
    return divisors

# Source - https://stackoverflow.com/a/70635601
# Posted by Alain T., modified by community. See post 'Timeline' for change history
# Retrieved 2026-04-30, License - CC BY-SA 4.0

import numpy as np

def npDivs(N):
    divs = np.arange(1,int(N**0.5)+1)  # potential divisors up to √N
    divs = divs[N % divs==0]           # divisors
    comp = N//divs[::-1]               # complement quotients
    return np.concatenate((divs,comp[divs[-1]==comp[0]:])) # combined
    

MONSTER = 808017424794512875886459904961710757005754368000000000
print( npDivs(MONSTER))