#!/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