# Spiral primes (Euler Problem 58)¶

Back to overview.

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%?

In [1]:

    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])

In [2]:

    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))

In [3]:

    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()

    n: 26241 count_total: 52481 count_primes: 5248 ratio: 0.09999809454850327

In [4]:

    print(s)
    assert(s == 26241)

    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.

In [12]:

    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.")

    n: 26641 count_total: 53281 count_primes: 5328 ratio: 0.09999812315834913
    26641
    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.

In [10]:

    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))

    n: 26241 count_total: 52481 count_primes: 5248 ratio: 0.09999809454850327
    26241
    For 561 sicp says True and real test says False.
    For 1105 sicp says True and real test says False.
    For 1729 sicp says True and real test says False.
    For 2465 sicp says True and real test says False.
    For 2821 sicp says True and real test says False.
    For 6601 sicp says True and real test says False.

In [ ]: