https://projecteuler.net/problem=58
Starting with 1 and spiralling anticlockwise in the following way, a square spiral with side length 7 is formed.
37 36 35 34 33 32 31
38 17 16 15 14 13 30
39 18 5 4 3 12 29
40 19 6 1 2 11 28
41 20 7 8 9 10 27
42 21 22 23 24 25 26
43 44 45 46 47 48 49
It is interesting to note that the odd squares lie along the bottom right diagonal, but what is more interesting is that 8 out of the 13 numbers lying along both diagonals are prime; that is, a ratio of 8/13 ≈ 62%.
If one complete new layer is wrapped around the spiral above, a square spiral with side length 9 will be formed. If this process is continued, what is the side length of the square spiral for which the ratio of primes along both diagonals first falls below 10%?
def get_last_corner_value(side_length):
return side_length * side_length
assert(get_last_corner_value(1) == 1)
assert(get_last_corner_value(3) == 9)
assert(get_last_corner_value(5) == 25)
def get_corner_values(side_length):
if side_length == 1:
return [1]
return [get_last_corner_value(side_length) - i * (side_length - 1)
for i in range(0, 4)][::-1]
assert(get_corner_values(1) == [1])
assert(get_corner_values(3) == [3, 5, 7, 9])
assert(get_corner_values(5) == [13, 17, 21, 25])
assert(get_corner_values(7) == [31, 37, 43, 49])
def get_diagonal_values(side_length):
return [corner_value
for length in range(1, side_length + 1, 2)
for corner_value in get_corner_values(length)
]
assert(get_diagonal_values(1) == [1])
assert(get_diagonal_values(3) == [1, 3, 5, 7, 9])
assert(get_diagonal_values(5) == [1, 3, 5, 7, 9, 13, 17, 21, 25])
def expmod(base, exp, m):
if exp == 0:
return 1
if (exp % 2 == 0):
return (expmod(base, exp // 2, m) ** 2 % m)
return (base * expmod(base, exp - 1, m) % m)
def fermat_test_sicp(n):
if n < 40:
return n in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 39]
a = n - 3
if not expmod(a, n, n) == a:
return False
a = n - 5
if not expmod(a, n, n) == a:
return False
a = n - 7
if not expmod(a, n, n) == a:
return False
a = n - 11
if not expmod(a, n, n) == a:
return False
a = n - 13
if not expmod(a, n, n) == a:
return False
a = n - 17
if not expmod(a, n, n) == a:
return False
a = n - 19
if not expmod(a, n, n) == a:
return False
a = n - 23
if not expmod(a, n, n) == a:
return False
a = n - 29
if not expmod(a, n, n) == a:
return False
a = n - 31
if not expmod(a, n, n) == a:
return False
a = n - 37
if not expmod(a, n, n) == a:
return False
a = n - 39
if not expmod(a, n, n) == a:
return False
return True
def trial_division(n):
a = []
if n % 2 == 0:
a.append(2)
while n % 2 == 0:
n //= 2
f = 3
while f * f <= n:
if n % f == 0:
a.append(f)
while n % f == 0:
n //= f
else:
f += 2
if n != 1:
a.append(n)
return a
def is_prime(n):
if n == 1:
return False
return trial_division(n)[0] == n
assert(fermat_test_sicp(3))
assert(fermat_test_sicp(107))
assert(fermat_test_sicp(108) == False)
assert(fermat_test_sicp(109))
def get_solution():
n = 1
count_primes = 0
count_total = 0
while True:
for v in get_corner_values(n):
count_total += 1
if is_prime(v):
count_primes += 1
ratio = count_primes / count_total
if ratio != 0 and ratio < 0.10:
print("n: {} count_total: {} count_primes: {} ratio: {}".format(n, count_total, count_primes, ratio))
return n
n += 2
s = get_solution()
print(s)
assert(s == 26241)
I actually got different results for the Fermat test and for the prime test which relies on reliable computation. I will actually try to solve the problem with ramdonized Fermat test now.
import random
def is_prime(n):
if n == 1:
return False
for _ in range(10):
a = n - random.randint(1, n - 1)
if not expmod(a, n, n) == a:
return False
return True
print(get_solution())
print("Expected n is 26241.")
Something seems to be off with my Fermat test...
Seems like there are systematic errors with the Fermat tests. Certain primes cannot be deteced.
Try this algorithm instead https://en.wikipedia.org/wiki/Miller–Rabin_primality_test.
carmichael_number = [561, 1105, 1729, 2465, 2821, 6601]
def little_fermat(a, p): # SICP uses this
return (a ** p) % p == a % p
def little_fermat_simple(a, p): # Wikipedia uses this
assert(a < p)
return (a ** (p - 1)) % p == 1
def expmod(base, exp, m):
if exp == 0:
return 1
if (exp % 2 == 0):
return (expmod(base, exp // 2, m) ** 2 % m)
return (base * expmod(base, exp - 1, m) % m)
def fermat_test(n):
ps = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 39]
if n <= ps[-1]:
return n in ps
for a in ps:
#if not little_fermat_simple(a, n):
if expmod(a, n - 1, n) != 1:
#print("Fermat witness for {} is {}.".format(n + 1, a))
return False
return True
def is_prime(n):
return fermat_test(n)
print(get_solution())
for n in carmichael_number:
r_sicp = fermat_test_sicp(n)
r = fermat_test(n)
print("For {} sicp says {} and real test says {}.".format(n, r_sicp, r))