Prime factorization-related functions

integer_factorization.py

from math import sqrt
from typing import List, Tuple, Union


def gcd(a: int, b: int) -> int:
    """
    Return the greatest common divisor of integers `a` and `b`

    >>> gcd(3, 6)
    3
    >>> gcd(4, 2)
    2
    >>> gcd(121, 2)
    1
    >>> gcd(3, 19)
    1
    >>> gcd(0, 8)
    8
    >>> gcd(0, 0)
    0
    >>> gcd(13, 0)
    13
    >>> gcd(0, -7)
    7
    >>> gcd(-3, -1)
    1
    """
    # Coerce inputs to nonnegative integers
    a *= -1 if a < 0 else 1
    b *= -1 if b < 0 else 1

    if a == 0 or b == 0:
        return max(a, b)
    if a == 1 or b == 1:
        return 1

    return gcd(min(a, b), max(a, b) % min(a, b))


def lcm(a: int, b: int) -> int:
    """
    Return the least common multiple of nonzero integers `a` and `b`

    >>> lcm(3, 6)
    6
    >>> lcm(20, 28)
    140
    >>> lcm(39, 45)
    585
    >>> lcm(-12, 10)
    60
    >>> lcm(-15, -30)
    30
    >>> lcm(0, -5)
    Traceback (most recent call last):
        ...
    AssertionError: `a` and `b` must be nonzero integers
    """
    assert a != 0 and b != 0, "`a` and `b` must be nonzero integers"

    # Coerce inputs to nonnegative integers
    a *= -1 if a < 0 else 1
    b *= -1 if b < 0 else 1

    return a * b // gcd(a, b)


def is_coprime(a: int, b: int) -> bool:
    """
    Return whether nonnegative integers `a` and `b` are relatively prime

    >>> is_coprime(2, 11)
    True
    >>> is_coprime(9, 3)
    False
    >>> is_coprime(0, 1)
    True
    >>> is_coprime(-1, 0)
    True
    >>> is_coprime(-12, 6)
    False
    >>> is_coprime(-10, -5)
    False
    >>> is_coprime(-10, -3)
    True
    """
    return gcd(a, b) == 1


def primality_check(n: int) -> Tuple[bool, int]:
    """
    Return tuple containing info about primality of `n`:
    - bool, whether nonnegative integer `n` is prime
    - integer, depending on primality:
        - if prime or less than 2, `n`
        - if composite, the smallest factor larger than 1 that divides `n`

    >>> primality_check(0)
    (False, 0)
    >>> primality_check(1)
    (False, 1)
    >>> primality_check(2)
    (True, 2)
    >>> primality_check(3)
    (True, 3)
    >>> primality_check(73)
    (True, 73)
    >>> primality_check(1223)
    (True, 1223)
    >>> primality_check(799)
    (False, 17)
    >>> primality_check(-799)
    Traceback (most recent call last):
        ...
    AssertionError: `n` must be a nonnegative integer
    """
    assert n >= 0, "`n` must be a nonnegative integer"

    if n < 2:
        return (False, n)
    if n % 2 == 0 and n > 2:
        return (False, 2)

    k, max_smaller_divisor = 3, (sqrt(n) // 1)
    while k <= max_smaller_divisor:
        if n % k == 0:
            return (False, k)
        k += 1
    return (True, n)


def prime_factorize_tree(n: int) -> Union[int, List[Union[int, List]]]:
    """
    Recursively builds prime factorization tree for nonnegative integer `n`

    >>> prime_factorize_tree(3)
    3
    >>> prime_factorize_tree(2)
    2
    >>> prime_factorize_tree(9)
    [9, [3, 3]]
    >>> prime_factorize_tree(30)
    [30, [2, [15, [3, 5]]]]
    >>> prime_factorize_tree(73)
    73
    >>> prime_factorize_tree(799)
    [799, [17, 47]]
    >>> prime_factorize_tree(1)
    Traceback (most recent call last):
        ...
    AssertionError: `n` must be an integer greater than 1
    >>> prime_factorize_tree(-1)
    Traceback (most recent call last):
        ...
    AssertionError: `n` must be an integer greater than 1
    >>> prime_factorize_tree(0)
    Traceback (most recent call last):
        ...
    AssertionError: `n` must be an integer greater than 1
    """
    assert n > 1, "`n` must be an integer greater than 1"

    is_prime, divisor = primality_check(n)

    if is_prime:
        return n

    return [n, [prime_factorize_tree(divisor), prime_factorize_tree(n // divisor)]]

prime_factors.c

#include "prime_factors.h"
#include <stdint.h>
#include <stdbool.h>
#include <math.h>

static inline uint64_t max_possible_divisor(uint64_t n)
{
    return (uint64_t)floor(sqrt((double)n));
}

static uint64_t next_prime_divisor(uint64_t n)
{
    for (uint64_t k = 2, k_bound = max_possible_divisor(n); k <= k_bound; k++)
    {
        if (n % k == 0)
        {
            return next_prime_divisor(k);
        }
    }
    return n;
}

size_t find_factors(uint64_t n, uint64_t factors[static MAXFACTORS])
{
    size_t len = 0;
    uint64_t k;
    while (n > 1 && len < MAXFACTORS)
    {
        {
            factors[len++] = (k = next_prime_divisor(n));
            n /= k;
        }
    }
    return len;
}

prime_factors.h

#ifndef PRIME_FACTORS_H
#define PRIME_FACTORS_H

#include <stdint.h>
#include <stddef.h>

#define MAXFACTORS 10

size_t find_factors(uint64_t n, uint64_t factors[static MAXFACTORS]);

#endif