```#!/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):
# 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

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
for i in range(len(p_list)):
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
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
```