# nayuki/Project-Euler-solutions

Switch branches/tags
Nothing to show
Fetching contributors…
Cannot retrieve contributors at this time
89 lines (78 sloc) 3.47 KB
 # # Solution to Project Euler problem 451 # Copyright (c) Project Nayuki. All rights reserved. # # https://www.nayuki.io/page/project-euler-solutions # https://github.com/nayuki/Project-Euler-solutions # import array, eulerlib, itertools # Let n be an arbitrary integer such that n >= 3. # When we say that the modular inverse of m modulo n equals m itself, # the formula is m^-1 = m mod n, which is equivalent to m^2 = 1 mod n. # # We know that if n is prime, then m^2 = 1 mod n has exactly two solutions: # m = 1, n-1. It is easy to verify that these two numbers are solutions. # The equation factorizes as (m - 1)(m + 1) = 0 mod n. Because n is prime, # the numbers form a field, and there are no zero divisors (two arbitrary # non-zero numbers x and y such that xy = 0). Hence 1 and -1 mod n are # the only possible solutions to the equation. (Note that for the excluded # special prime case where n = 2, the solutions 1 and -1 are the same number.) # # Suppose we can find the smallest prime factor of n quickly. (Note that if n is # prime, then the smallest prime factor is n itself.) This can be achieved by # building a table ahead of time, using a modification of the sieve of Eratosthenes. # # Suppose that for every n' < n, we know the set of solutions to m^2 = 1 mod n'. # This means whenever we solve the equation for the number n, we save its solutions # in an ever-growing list, so that when we work on the next value of n we can access # all possible smaller solutions. This is essentially an argument by strong induction. # # Let p be the smallest prime factor of n. If p = n, then the set of # solutions is {1, n - 1}, and we are finished with this value of n. # # Otherwise p < n, and obviously n is an integer multiple of p. Because we are looking # for values of m such that m^2 = 1 mod n, these candidate m values also must satisfy # m^2 = 1 mod k for any k that divides n (i.e. k is a factor of n). We look at the set # of solutions for the modulus k = n/p, which has already been solved because k < n. # We know that any solution modulo n must be congruent to these solutions modulo k. # Hence we can try to extend and check these old solutions by brute force. Namely, suppose # m' is a solution modulo k. Then we check the sequence m = m' + 0k, m' + 1k, m' + 2k, ..., # m' + (p-1)k modulo n. Because p is usually a small number, this isn't a lot of work to do. def compute(): LIMIT = 20000000 # Build table of smallest prime factors smallestprimefactor = array.array("L", itertools.repeat(0, LIMIT + 1)) end = eulerlib.sqrt(len(smallestprimefactor) - 1) for i in range(2, len(smallestprimefactor)): if smallestprimefactor[i] == 0: smallestprimefactor[i] = i if i <= end: for j in range(i * i, len(smallestprimefactor), i): if smallestprimefactor[j] == 0: smallestprimefactor[j] = i # Returns all the solutions (in ascending order) such that # for each k, 1 <= k < n and k^2 = 1 mod n. def get_solutions(n): if smallestprimefactor[n] == n: # n is prime return (1, n - 1) else: temp = [] p = smallestprimefactor[n] sols = solutions[n // p] for i in range(0, n, n // p): for j in sols: k = i + j if k * k % n == 1: temp.append(k) return tuple(temp) # Process every integer in range solutions = [(), (), (1,)] ans = 0 for i in range(3, LIMIT + 1): sols = get_solutions(i) if i <= LIMIT // 2: solutions.append(sols) ans += sols[-2] # Second-largest solution return str(ans) if __name__ == "__main__": print(compute())