This is the Python code module I have developed for number theoretical problems:

#!/usr/local/bin/python3

from fractions import gcd
from math import log, ceil

class DataError(Exception):
    def __init__(self, value):
        self.value = value

    def __str__(self):
        return self.value

# I picked this one up from ActiveState. Thanks, guys.
# The classic Sieve of Eratosthenes
# Input n -- the largest number in the sieve
# Return -- the sieved primes up to n as an array

def eratosthenes_sieve(n):
    sq_root, half = n ** 0.5, (n + 1) // 2 - 1
    i, m = 0, 3
    s = [j for j in range(3, n + 1, 2)]
    while m <= sq_root:
        if s[i]:
            j = (m * m - 3) // 2
            while j < half:
                s[j] = 0
                j += m
        i += 1
        m = 2 * i + 3

    return [2] + [j for j in s if j]

# Input n -- the number to find the prime factors of
# Return -- the prime factors (repeated, as necessary) as an array

def prime_factors(n):
    while True:
        if rabin_miller(n, 1.0e-7):
            result = [n]
            break
        if n & 0x01:
            p_list = eratosthenes_sieve(n // 3)
        else:
            p_list = eratosthenes_sieve(n // 2)
        result, done = [], False
        for a_p in p_list:
            while n % a_p == 0:
                result.append(a_p)
                n //= a_p
                if n == 1:
                    done = True
                    break
            if done:
                break
        break

    return result

# Input n -- the number to find the prime factors of
# Return -- the prime factors (not repeated) as an array

def unique_prime_factors(n):
    result, previous_factor = [], 0
    for factor in prime_factors(n):
        if factor != previous_factor:
            result.append(factor)
            previous_factor = factor

    return result

# Input n -- the number to find the Euler totient of
# Return -- the totient of n

def totient(n):
    if rabin_miller(n, 1.0e-7):
        result = n - 1
    else:
        result = n
        for a_p in unique_prime_factors(n):
            result -= result // a_p

    return result

# From Rosen's Number Theory book
# Multiplicative inverse via extended Euclid's algorithm
# Input a -- The integer to find the multiplicative inverse of
# Input m -- The modulus of the inverse operation
# Return -- The multiplicative inverse of a (mod m)

def m_i_euclid(a, m):
    if gcd(a, m) > 1:
        raise DataError(repr(a) + ' and ' + repr(m) + ' are not relatively prime!')

    m_in, x, x_back = m, 0, 1
    while m != 0:
        q, r = divmod(a, m)
        a, m = m, r
        x, x_back = x_back - q * x, x
    
    return x_back % m_in

# Find s and t in a * s + b * t = gcd(a, b)
# Input a -- First integer
# Input b -- Second integer
# Return -- s and t

def extended_euclid(a, b):
    x, x_back, y, y_back = 0, 1, 1, 0
    while b != 0:
        q, r = divmod(a, b)
        a, b = b, r
        x, x_back = x_back - q * x, x
        y, y_back = y_back - q * y, y

    return x_back, y_back

# Fast and efficient algorithm for computing b^e (mod m)
# Input b -- the base, or number to be exponentiated
# Input e -- the exponent to raise the base, b, to
# Input m -- the modulus of the exponentiation problem, b^e (mod m)
# Return -- b^e (mod m)
def pow_mod(b, e, m):
    result = 1
    while e > 0:
        if e & 1:
            result = (result * b) % m
        e >>= 1
        b = (b * b) % m

    return result

# Check for pairwise primality of a set of integers
# Input ps -- array of integers to test
# Return -- True if integers are pairwise prime, False otherwise

def pairwise_prime(ps):
    l, result = len(ps), True
    
    for i in range(l - 1):
        p1 = ps[i]
        for j in range(i + 1, l):
            p2 = ps[j]
            if gcd(p1, p2) != 1:
                result = False
                break
        if not result:
            break
            
    return result
    
# Chinese Remainder Theorem for a single variable group of congruences
# Input a_of_t -- array of tuples of congruences, with each subarray of the form (xi, modi)
# Return -- the simultaneous solution of the congruences
    
def crt_1_variable(a_of_t):
    result = 0

    ps = [tup[1] for tup in a_of_t]
    if not pairwise_prime(ps):
        raise DataError(repr(ps) + ' not pairwise prime!')
    m = 1
    for p in ps:
        m *= p
    ms = [m // p for p in ps]

    for i in range(len(ps)):
        md = ps[i]
        m_i = m_i_euclid(ms[i], md)
        result += a_of_t[i][0] * ms[i] * m_i
    return result % m
    
# Chinese Remainder Theorem for two variables with the same modulus
# Inputs a, b, c, d, e, f, m from
# a * x + b * y = e (mod m) and c * x + d * y = f (mod m)
# Return -- the array x and y

def crt_2_variables(a, b, c, d, e, f, m):
    delta = a * d - b * c
    if gcd(delta, m) != 1:
        raise DataError("Discriminant is not relatively prime to modulus!")
    delta_i = m_i_euclid(delta, m)
    x = (delta_i * (d * e - b * f)) % m
    y = (delta_i * (a * f - c * e)) % m

    return x, y

# Pull the factors of 2 out of a number, n - 1
# Input n -- the number to split into n - 1 = 2^s * t
# Return -- s and t
def miller_s_and_t(n):
    s, t = 0, n - 1
    while t % 2 == 0:
        t //= 2
        s += 1

    return s, t

# Miller's test for primality with respect to a given base
# Input n -- the number to test for primality
# Input base -- the base for the test
# Input s -- n - 1 = 2**s * t
# Input t -- n - 1 = 2**s * t
# Return -- True if n is prime to base, False otherwise

def millers_test(n, base, s, t):
    rv, n_minus_1 = False, n - 1

    base_to_power = pow_mod(base, t, n)
    if base_to_power == 1 or base_to_power == n_minus_1:
        # Obeys Fermat's Little Theorem, so probably prime. 
        rv = True
    else:
        for i in range(1, s):
            base_to_power = pow_mod(base_to_power, 2, n)
            if base_to_power == 1:
                # base-to_power's square root is not + or - 1.
                break
            elif base_to_power == n_minus_1:
                # Obeys Fermat's Little Theorem, so probably prime.
                rv = True
                break

        # Composite.
    return rv

# Rabin-Miller probabilistic test for primality
# Input n -- the number to test for primality
# Input uncertainty -- the level of uncertainty of the prime test
# Return -- True if n is believed to be prime, False otherwise

def rabin_miller(n, uncertainty):
    rv = True

    if n != 2:
        count, base = int(log(uncertainty) / ceil(log(0.25))), 2
        if count < 1:
            count = 1

        s, t = miller_s_and_t(n)
        i = 1
        while i <= count:
            if base % n != 0:
                if not millers_test(n, base, s, t):
                    rv = False
                    break
                
                i += 1
            base += 1

    return rv

# Rabin-Miller deterministic test for primality of integers < 25.0e9
# Input n -- the number to test for primality
# Return -- True if prime, False otherwise

def rabin_miller_certain(n):
    if n >= 25.0e9:
        raise DataError(n, 'must be <', 25.0e9)

    strong_pseudoprime, rv = 3215031751, True
    if n == strong_pseudoprime:
        rv = False
    elif n != 2:
        s, t = miller_s_and_t(n)
        for base in (2, 3, 5, 7):
            if not millers_test(n, base, s, t):
                rv = False
                break

    return rv

# Second order Newton method (Halley's method) for square roots of large numbers
# Input n -- the number to take the root of
# There are three return values:
#    the best approximation to the square root (integer)
#    either 'Exact' or 'Floor' (string)
#    the number of iterations (integer)
# A return of 'Exact' means that n is a perfect square and that the root is the true root
# A return of 'Floor' means the root is within 1 of the square root and is greater than it

APPROX_SQRT_MAX_ITERS = 1000
def approx_sqrt(n):
    msg, n_3, iters = 'Not found', 3 * n, 0
    # base 2 logarithm approximation:
    # bin prepends '0b' to its string output, so must subtract 2 
    x = n >> (len(bin(n)) - 2) // 2
    for i in range(1, APPROX_SQRT_MAX_ITERS + 1):
        iters += 1
        x_next = (x * (x * x + n_3)) // (3 * x * x + n)
        if abs(x_next - x) < 1:
            if x_next * x_next == n:
                msg = 'Exact'
            else:
                msg = 'Floor'
            break
        x = x_next

    return x_next, msg, iters

# The Baby Step Giant Step method for the discrete logarithm using a hash for storage. Runs in O(sqrt(n)) time
# Input beta -- the number to take the discrete logarithm of
# Input alpha -- the base of the logarithm
# Input mod -- the modulus of the logarithm
# Return y where a^y = x (mod modulus)

def baby_step_giant_step(beta, alpha, mod):
    table = {}
    m, msg, iters = approx_sqrt(mod)
    if msg == 'Floor':
        m += 1
    for j in range(m):
        table[pow_mod(alpha, j, mod)] = j

    a_to_minus_m = m_i_euclid (alpha**m, mod)
    gamma, found = beta, False
    for i in range(m):
        for key, j in table.items():
            if key == gamma:
                result = i * m + j
                found = True
                break
        if found:
            break
        gamma = (gamma * a_to_minus_m) % mod
    return result

# Pollard's Rho factorization method
# Input n -- the number to attempt to factor
# Input starting_value -- the starting values for x and y
# Return -- either 'Failure' a factor of n

def pollard_factorization(n, starting_value):
    x = y = starting_value
    d = 1
    while d == 1:
        x = (x * x + 1) % n
        temp_y = (y * y + 1) % n
        y = (temp_y * temp_y + 1) % n
        d = gcd(abs(x - y), n)
    if d == n:
        result = 'Failure'
    else:
        result = d
    return d

# A function for the tortoise and hare and Brent algorithms.

def f(x, n):
    return (x * x + 1) % n

# Floyd's cycle-finding algorithm
# Input f -- a function of two variables that cycles at some point
# Input x0 -- a starting value for f
# Input n -- a constant second parameter of f
# Return -- lambda and mu where mu is the starting index of f's cycle and lambda is the period

def tortoise_and_hare(f, x0, n):
    mu, lam = 0, 1
    tortoise = f(x0, n)
    hare = f(tortoise, n)
    while tortoise != hare:
        tortoise = f(tortoise, n)
        hare = f(f(hare, n), n)

# At this point the hare and tortoise are separated by some multiple of the period, lambda, call it nu.
# But the hare is twice as many steps from the origin as the hare, at 2 * nu.
# Reset the tortoise to the starting point and let both the tortoise and hare step along until
# they are at the same value.

    tortoise = x0
    while tortoise != hare:
        tortoise = f(tortoise, n)
        hare = f(hare, n)
        mu += 1

# That position is mu, the start of the loop.
# Now reset just the hare to the tortoise and let the hare step along until the values are equal.
# That number of steps is lambda, the period.

    hare = f(tortoise, n)
    while tortoise != hare:
        hare = f(hare, n)
        lam += 1
        
    return lam, mu

# Brent's modification of Floyd's cycle-finding algorithm
# Input f -- a function of two variables that cycles at some point
# Input x0 -- a starting value for f
# Input n -- a constant second parameter of f
# Return -- lambda and mu where mu is the starting index of f's cycle and lambda is the period

def t_and_h_brent(f, x0, n):
    power = lam = 1
    tortoise, hare = x0, f(x0, n)
    while tortoise != hare:
        if power == lam:
            tortoise = hare
            power *= 2
            lam = 0
        hare = f(hare, n)
        lam += 1

# We now have determined the period, lambda.
# Next, position the tortoise at x0 and step the hare a period ahead, which won't
# initially be the same value.

    mu = 0
    tortoise = hare = x0
    for i in range(lam):
        hare = f(hare, n)

    while tortoise != hare:
        tortoise = f(tortoise, n)
        hare = f(hare, n)
        mu += 1

# Now we've stepped the hare to mu because the values agree.

    return lam, mu

# Pollard's algorithm for finding the discrete logarithm
# Input big_n -- the modulus of the cyclic group
# Input n -- the group's size
# Input alpha -- the base of the logarithm
# Input beta -- the number to take the logarithm of
# Return -- either the Left String "No solution" or the Right Integer logarithm

def pollard_log(N, n, alpha, beta):
    x, a, b = 1, 0, 0
    X, A, B = x, a, b
    done = False
    for i in range(n):
        x, a, b = new_xab(x, a, b, N, n, alpha, beta)
        X, A, B = new_xab(X, A, B, N, n, alpha, beta)
        X, A, B = new_xab(X, A, B, N, n, alpha, beta)
        # print('{0:3d}  {1:4d} {2:3d} {3:3d}  {4:4d} {5:3d} {6:3d}'.format((i + 1), x, a, b, X, A, B))
        if x == X:
            done = True
            break

    if not done:
        result = 'No solution'
    else:
        lhs, rhs = B - b, a - A
        if lhs < 0:
            lhs += n
        if rhs < 0:
            rhs += n
        gcd_lhs_rhs = gcd(lhs, rhs)
        gcd_lhs_n = gcd(lhs, n)
        lhs //= gcd_lhs_rhs
        rhs //= gcd_lhs_rhs
        n //= gcd_lhs_n
        inv = m_i_euclid(lhs, n)
        log_try = (inv * rhs) % n
        beta_test = pow_mod(alpha,log_try, N)
        if beta_test == beta:
            result = log_try
        else:
            result = 'Alternate solution to beta = ' + str(beta_test) + ', gamma = ' + str(log_try)
    return result

def new_xab(x, a, b, big_n, n, alpha, beta):
    xp3 = x % 3
    if xp3 == 0:
        new_x = x * x % N
        new_a = a * 2 % n
        new_b = b * 2 % n
    elif xp3 == 1:
        new_x = x * alpha % N
        new_a = (a + 1) % n
        new_b = b
    elif xp3 == 2:
        new_x = x * beta % N
        new_a = a
        new_b = (b + 1) % n
    return new_x, new_a, new_b

# The Dixon algorithm for factoring integers
# Input -- n, The number to factor
# Input -- starting_num, The start of the range to check for possible congruence candidates
# Input -- count, The number of candidate numbers in the list to pair-up
# Return -- Either "No solution" or f1, f2 as factors of n

def dixon_algorithm(n, starting_number, count):
    p_list, possibles, test_num, i = (2, 3, 5, 7), [], starting_number, 0
    
    while i < count:
        temp = (test_num * test_num) % n
        if is_possible(temp, p_list):
            possibles.append(test_num)
            i += 1
        test_num += 1
    
    pairs = [(possibles[i], possibles[j]) for i in range(1, count) for j in range(i)]
    goods = [a_pair for a_pair in pairs if test_pair(a_pair, n, p_list)]
    goods_keepers =[a_good for a_good in goods if get_factors(a_good, n, p_list)[0] != 1 and get_factors(a_good, n, p_list)[0] != n]

    if len(goods_keepers) == 0:
        result = 'No solution'
    else:
        result = get_factors(goods_keepers[0], n, p_list)
    
    return result

# Given an integer, see if its only prime factors are 2, 3, 5, and 7
# Input n -- the integer to check for prime factors
# Input p_list -- the tuple (2, 3, 5, 7)
# Return -- True if all the prime factors are <= 7, False otherwise

def is_possible(n, p_list):
    result = False
    for a_prime in p_list:
        d, r = divmod(n, a_prime)
        while r == 0:
            if d == 1:
                result = True
                break
            n = d
            d, r = divmod(n, a_prime)
        if result == True:
            break
        
    return result

# Collect the prime factors 2, 3, 5, and 7, repeated, as necessary, returning them in an array of tuples
# Input n -- the integer to check for prime factors
# Input p_list -- the tuple (2, 3, 5, 7)
# Return -- An array of 4 tuples including missing primes

def padded_possible(n, p_list):
    result = []
    for a_prime in p_list:
        count = 0
        d, r = divmod(n, a_prime)      
        while r == 0:
            count += 1
            n = d
            d, r = divmod(n, a_prime)
        result.append((a_prime, count))
        
    return result
  
# Test a pair of integers to see if their squared modulus is a perfect square
# Input a_pair -- A tuple of two integers to test
# Input n -- the number to factor
# Input p_list -- the tuple (2, 3, 5, 7)
# Return -- True if the input pair satisfies the Dixon criterion, False otherwise
   
def test_pair(a_pair, n, p_list):
    n1, n2 = a_pair
    result = True
    n1_padded = padded_possible(n1 * n1 % n, p_list)
    n2_padded = padded_possible(n2 * n2 % n, p_list)
    for i in range(len(p_list)):
        if (n1_padded[i][1] + n2_padded[i][1]) & 1:
            result = False
            break

    return result

# Given a pair of integers satisfying the Dixon criterion, find the factors of the number n
# Input a_pair -- A tuple of two integers satisfying the Dixon criterion
# Input n -- the number to factor
# Input p_list -- the tuple (2, 3, 5, 7)
# Return -- The integer factors of n

def get_factors(a_pair, n, p_list):
    n1, n2 = a_pair
    d1 = n1 * n2 % n
    n1_padded = padded_possible(n1 * n1 % n, p_list)
    n2_padded = padded_possible(n2 * n2 % n, p_list)
    comb = [(n1_padded[i][0], (n1_padded[i][1] + n2_padded[i][1]) // 2) for i in range(len(p_list))]
    d2 = 1
    for a_pair_2 in comb:
        d2 *= a_pair_2[0]**a_pair_2[1]
    f1 = gcd(d1 + d2, n)
    f2 = n // f1
    
    return f1, f2

# Solve the quadratic residue equation x^2 = n (mod p) for x via the Tonelli-Shanks algorithm
# Input -- p, the prime for the residue calculations
# Input -- n, the quadratic residue
# Return -- either "n is not a quadratic residue modulo p" or s1 and s2 where s1 and s2 are solutions 

def shanks_tonelli(p, n):
    if n ** ((p - 1) // 2) % p == 1:
        s, q = miller_s_and_t(p)
        pm1 = p - 1
        z = 0
        for z in range(2, p):
            if z**(pm1 // 2) % p == pm1:
                break
        c = z**q % p
        r = n**((q + 1) // 2) % p
        t = n**q % p
        m = s
        
        while t != 1:
            i = 0
            for i in range(1, m):
                if t**(2**i) % p == 1:
                    break
            b = c**(2**(m - i - 1)) % p
            r = r * b % p
            t = t * b * b % p
            m = i
            c = b * b % p
        return r, p - r
    else:
        result = str(n) + ' is not a quadratic residue modulo ' + str(p)

    return result

if __name__ == "__main__":
    import sys
    n = 252
    result = eratosthenes_sieve(n)
    print('Primes <=', n, ':', result)
    result = prime_factors(n)
    print('\nPrime factors of', n, 'are', result)
    result = unique_prime_factors(n)
    print('Unique prime factors of', n, 'are', result)
    result = totient(n)
    print('\nThe totient of', n, '=', result)
    m, a, m_i = 7, 11, 0
    try:
        m_i = m_i_euclid(a, m)
        print('\nThe multiplicative inverse of', a, 'modulus', m, '=', m_i)
    except DataError as e:
        print('\nDataError exception occurred:', e.value)
    print('Check:', a, '*', m_i, '=', a * m_i % 7, '(mod', m, ')')
        
    a, b = 9, 11
    s, t = extended_euclid(a, b)
    print('\nFor the integers', a, 'and', b, 'with gcd =', gcd(a, b))
    print(a, '*', s, '+', b, '*', t, '=', gcd(a, b))
    
    b, e, m = 3, 10, 17
    result = pow_mod(b, e, m)
    print('\npow_mod:', b, '**', e, '=', result, '(mod', m, ')')

    ps = [10, 9, 7, 25]
    result = pairwise_prime(ps)
    print('\nPairwise prime test on:', ps, 'is', result)
    ps = [4, 9, 7, 25]
    result = pairwise_prime(ps)
    print('Pairwise prime test on:', ps, 'is', result)

    t_of_t = ((1, 3), (2, 5), (3, 7))
    try:
        result = crt_1_variable(t_of_t)
        print('\nChinese Remainder Theorem on', t_of_t, '=', result)
        print('Check:')
        print(result, '(mod', t_of_t[0][1], ') =', result % t_of_t[0][1])
        print(result, '(mod', t_of_t[1][1], ') =', result % t_of_t[1][1])
        print(result, '(mod', t_of_t[2][1], ') =', result % t_of_t[2][1])
    except DataError as e:
        print('\nDataError exception occurred:', e.value)

    a, b, c, d, e, f, m = 3, 4, 2, 5, 5, 7, 13
    x, y = crt_2_variables(a, b, c, d, e, f, m)
    
    print('\nChinese Remainder Theorem in two variables:')
    print(a, '* x +', b, '* y =', e, '(mod', m, ')')
    print(c, '* x +', d, '* y =', f, '(mod', m, ')')
    print('x =', x, '(mod', m, ')', 'y =', y, '(mod', m, ')')
    print('Check:')
    print(a, '*', x, '+', b, '*', y, '=', (a * x + b * y) % m, '(mod', m, ')')
    print(c, '*', x, '+', d, '*', y, '=', (c * x + d * y) % m, '(mod', m, ')')

    n = 997
    s, t = miller_s_and_t(n)
    print('\nMiller s and t on', 997)
    print('s =', s, 'and t =', t)
    print('Check:', '(2 **', s, ') *', t, '=', 2**s * t)
    
    base = 2
    result = millers_test(n, base, s, t)
    print("\nMiller's test on", n, 'is', result )
    
    uncertainty = 1.0e-5
    result = rabin_miller(n, uncertainty)
    print('\nRabin-Miller test on', n, 'is', result, 'with uncertainty', uncertainty )
    
    result = rabin_miller_certain(n)
    print('\nRabin-Miller certain test on', n, 'is', result)
    
    n = 50000000
    result, msg, iters = approx_sqrt(n)
    print('\nInteger (', msg, ') square root of', n, '=', result, 'after', iters, 'iterations')
    print('Check:')
    print(result, '** 2 = ', result**2)
    print(result + 1, '** 2 = ', (result + 1)**2)
    
    beta, alpha, mod = 57, 3, 113
    result = baby_step_giant_step(beta, alpha, mod)
    print('\nLog to base', alpha, 'of', beta, 'mod(', mod, ') =', result)
    print('Check:', alpha, '**', result, '(mod', mod, ') =', alpha**result % mod)
    
    n, starting_value = 999, 3
    result = pollard_factorization(n, starting_value)
    print('\nPollard factorization of', n, 'with starting value of', starting_value, '=', result)
    print('Check:', n, 'divided by', result, '=', n // result)

#    x0, n = 2, 8031
#    This should work. Must be a bug in Python:
#    lam, mu = tortoise_and_hare(f, x0, n)
#    print('\nFor tortoise_and_hare with f = (x**2 + 1) % n,', 'x0 = ', x0, 'and', 'n =', n)
#    print('mu =', mu, 'and lambda =', lam)

#    x0, n = 2, 1011
#    lam, mu = t_and_h_brent(f, x0, n)
#    print('\nFor t_and_h_brent with f = (x**2 + 1) % n,', 'x0 = ', x0, 'and', 'n =', n)
#    print('mu =', mu, 'and lambda =', lam)
    print('\nPollard integer logarithm:')
    N, n, alpha, beta = 383, 191, 2, 228
    result = pollard_log(N, n, alpha, beta)
    print('log(', beta, ') to base', alpha, 'mod(', N, ') =', result)
    print('Check:', alpha, '**', result, '(mod', N, ') ==', 2**result % N)
    N, n, alpha, beta = 1019, 1018, 2, 1024
    result = pollard_log(N, n, alpha, beta)
    print('log(', beta, ') to base', alpha, 'mod(', N, ') =', result)
    print('Check:', alpha, '**', 10, '(mod', N, ') ==', 2**10 % N)

    n, starting_number, count = 100160063, 1200000, 4
    print('\nDixon factoring algorithm')
    print('With n =', n, 'and starting from', starting_number, 'and collecting', count, 'possible numbers:')
    result = dixon_algorithm(n, starting_number, count)
    if result == 'No solution':
        print('the Dixon factoring algorithm yields no solution')
    else:
        print('the factors are', result[0], 'and', result[1])
        print('Check:', result[0], '*', result[1], '=', result[0] * result[1])
   
    p, n = 1000033, 5
    print('\nShanks Tonelli')
    result = shanks_tonelli(p, n)
    print('shanks_tonelli(', p, ',', n, ') fails with message:', result)
    
    n = 6
    r1, r2 = shanks_tonelli(p, n)
    print('\nshanks_tonelli(', p, ',', n, ') =', r1, ',', r2)
    print('Check:', r1, '** 2 (mod', p, ')=', r1**2 % p)
    print('Check:', r2, '** 2 (mod', p, ')=', r2**2 % p)

End of Code:

Primes <= 252 : [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251]

Prime factors of 252 are [2, 2, 3, 3, 7]
Unique prime factors of 252 are [2, 3, 7]

The totient of 252 = 72

The multiplicative inverse of 11 modulus 7 = 2
Check: 11 * 2 = 1 (mod 7 )

For the integers 9 and 11 with gcd = 1
9 * 5 + 11 * -4 = 1

pow_mod: 3 ** 10 = 8 (mod 17 )

Pairwise prime test on: [10, 9, 7, 25] is False
Pairwise prime test on: [4, 9, 7, 25] is True

Chinese Remainder Theorem on ((1, 3), (2, 5), (3, 7)) = 52
Check:
52 (mod 3 ) = 1
52 (mod 5 ) = 2
52 (mod 7 ) = 3

Chinese Remainder Theorem in two variables:
3 * x + 4 * y = 5 (mod 13 )
2 * x + 5 * y = 7 (mod 13 )
x = 7 (mod 13 ) y = 9 (mod 13 )
Check:
3 * 7 + 4 * 9 = 5 (mod 13 )
2 * 7 + 5 * 9 = 7 (mod 13 )

Miller s and t on 997
s = 2 and t = 249
Check: (2 ** 2 ) * 249 = 996

Miller's test on 997 is True

Rabin-Miller test on 997 is True with uncertainty 1e-05

Rabin-Miller certain test on 997 is True

Integer ( Floor ) square root of 50000000 = 7071 after 3 iterations
Check:
7071 ** 2 =  49999041
7072 ** 2 =  50013184

Log to base 3 of 57 mod( 113 ) = 100
Check: 3 ** 100 (mod 113 ) = 57

Pollard factorization of 999 with starting value of 3 = 111
Check: 999 divided by 111 = 9

Pollard integer logarithm:
log( 228 ) to base 2 mod( 383 ) = 110
Check: 2 ** 110 (mod 383 ) == 228
log( 1024 ) to base 2 mod( 1019 ) = Alternate solution to beta = 5, gamma = 10
Check: 2 ** 10 (mod 1019 ) == 5

Dixon factoring algorithm
With n = 100160063 and starting from 1200000 and collecting 4 possible numbers:
the factors are 10009 and 10007
Check: 10009 * 10007 = 100160063

Shanks Tonelli
shanks_tonelli( 1000033 , 5 ) fails with message: 5 is not a quadratic residue modulo 1000033

shanks_tonelli( 1000033 , 6 ) = 348481 , 651552
Check: 348481 ** 2 (mod 1000033 )= 6
Check: 651552 ** 2 (mod 1000033 )= 6