2018-12-24 16:57:03 +01:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Spiral primes (Euler Problem 58)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"[https://projecteuler.net/problem=58](https://projecteuler.net/problem=58)\n",
"\n",
"Starting with 1 and spiralling anticlockwise in the following way, a square spiral with side length 7 is formed.\n",
"\n",
2018-12-29 10:17:35 +01:00
" 37 36 35 34 33 32 31\n",
" 38 17 16 15 14 13 30\n",
" 39 18 5 4 3 12 29\n",
" 40 19 6 1 2 11 28\n",
" 41 20 7 8 9 10 27\n",
" 42 21 22 23 24 25 26\n",
" 43 44 45 46 47 48 49\n",
2018-12-24 16:57:03 +01:00
"\n",
"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%.\n",
"\n",
"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%?"
]
},
2018-12-29 10:17:35 +01:00
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"def get_last_corner_value(side_length):\n",
" return side_length * side_length\n",
"\n",
"assert(get_last_corner_value(1) == 1)\n",
"assert(get_last_corner_value(3) == 9)\n",
"assert(get_last_corner_value(5) == 25)\n",
"\n",
"\n",
"def get_corner_values(side_length):\n",
" if side_length == 1:\n",
" return [1]\n",
" return [get_last_corner_value(side_length) - i * (side_length - 1)\n",
" for i in range(0, 4)][::-1]\n",
"\n",
"assert(get_corner_values(1) == [1])\n",
"assert(get_corner_values(3) == [3, 5, 7, 9])\n",
"assert(get_corner_values(5) == [13, 17, 21, 25])\n",
"assert(get_corner_values(7) == [31, 37, 43, 49])\n",
"\n",
"\n",
"def get_diagonal_values(side_length):\n",
" return [corner_value\n",
" for length in range(1, side_length + 1, 2)\n",
" for corner_value in get_corner_values(length)\n",
" ]\n",
"\n",
"assert(get_diagonal_values(1) == [1])\n",
"assert(get_diagonal_values(3) == [1, 3, 5, 7, 9])\n",
"assert(get_diagonal_values(5) == [1, 3, 5, 7, 9, 13, 17, 21, 25])"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def expmod(base, exp, m):\n",
" if exp == 0:\n",
" return 1\n",
" if (exp % 2 == 0):\n",
" return (expmod(base, exp // 2, m) ** 2 % m)\n",
" return (base * expmod(base, exp - 1, m) % m)\n",
" \n",
"\n",
"def fermat_test(n):\n",
" if n < 40:\n",
" return n in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 39]\n",
" a = n - 3\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" a = n - 5\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" a = n - 7\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" a = n - 11\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" a = n - 13\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" a = n - 17\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" a = n - 19\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" a = n - 23\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" a = n - 29\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" a = n - 31\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" a = n - 37\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" a = n - 39\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" return True\n",
"\n",
"def trial_division(n):\n",
" a = [] \n",
" if n % 2 == 0:\n",
" a.append(2)\n",
" while n % 2 == 0:\n",
" n //= 2\n",
" f = 3\n",
" while f * f <= n:\n",
" if n % f == 0:\n",
" a.append(f)\n",
" while n % f == 0:\n",
" n //= f\n",
" else:\n",
" f += 2 \n",
" if n != 1:\n",
" a.append(n)\n",
" return a\n",
"\n",
"def is_prime(n):\n",
" if n == 1:\n",
" return False\n",
" return trial_division(n)[0] == n\n",
"\n",
"assert(fermat_test(3))\n",
"assert(fermat_test(107))\n",
"assert(fermat_test(108) == False)\n",
"assert(fermat_test(109))"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n: 26241 count_total: 52481 count_primes: 5248 ratio: 0.09999809454850327\n"
]
}
],
"source": [
"def get_solution():\n",
" n = 1\n",
" count_primes = 0\n",
" count_total = 0\n",
" while True:\n",
" for v in get_corner_values(n):\n",
" count_total += 1\n",
" if is_prime(v):\n",
" count_primes += 1\n",
" ratio = count_primes / count_total\n",
" if ratio != 0 and ratio < 0.10:\n",
" print(\"n: {} count_total: {} count_primes: {} ratio: {}\".format(n, count_total, count_primes, ratio))\n",
" return n\n",
" n += 2\n",
"\n",
"s = get_solution()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"26241\n"
]
}
],
"source": [
"print(s)\n",
"assert(s == 26241)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import random\n",
"\n",
"def is_prime(n):\n",
" if n == 1:\n",
" return False\n",
" for _ in range(100):\n",
" a = n - random.randint(1, n - 1)\n",
" if not expmod(a, n, n) == a:\n",
" return False\n",
" return True\n",
"\n",
"assert(fermat_test(3))\n",
"assert(fermat_test(107))\n",
"assert(fermat_test(108) == False)\n",
"assert(fermat_test(109))"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n: 26641 count_total: 53281 count_primes: 5328 ratio: 0.09999812315834913\n"
]
}
],
"source": [
"s = get_solution()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Something seems to be off with my Fermat test... \n",
"\n",
"Seems like there are systematic errors with the Fermat tests. Certain primes cannot be deteced.\n",
"\n",
"Try this algorithm instead [https://en.wikipedia.org/wiki/Miller– –
]
},
2018-12-24 16:57:03 +01:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
2018-12-29 10:17:35 +01:00
"completion_date": "Sat, 29 Dec 2018, 09:00",
2018-12-24 16:57:03 +01:00
"kernelspec": {
"display_name": "Python 3",
"language": "python3.6",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
},
2018-12-29 10:17:35 +01:00
"tags": [
"prime",
"fermat",
"spiral",
"diagonal",
"todo"
]
2018-12-24 16:57:03 +01:00
},
"nbformat": 4,
"nbformat_minor": 2
}
