354 lines
10 KiB
Plaintext
354 lines
10 KiB
Plaintext
{
|
||
"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",
|
||
" 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",
|
||
"\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%?"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 1,
|
||
"metadata": {
|
||
"collapsed": true
|
||
},
|
||
"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_sicp(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_sicp(3))\n",
|
||
"assert(fermat_test_sicp(107))\n",
|
||
"assert(fermat_test_sicp(108) == False)\n",
|
||
"assert(fermat_test_sicp(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": 12,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"n: 26641 count_total: 53281 count_primes: 5328 ratio: 0.09999812315834913\n",
|
||
"26641\n",
|
||
"Expected n is 26241.\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"import random\n",
|
||
"\n",
|
||
"def is_prime(n):\n",
|
||
" if n == 1:\n",
|
||
" return False\n",
|
||
" for _ in range(10):\n",
|
||
" a = n - random.randint(1, n - 1)\n",
|
||
" if not expmod(a, n, n) == a:\n",
|
||
" return False\n",
|
||
" return True\n",
|
||
"\n",
|
||
"print(get_solution())\n",
|
||
"print(\"Expected n is 26241.\")"
|
||
]
|
||
},
|
||
{
|
||
"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–Rabin_primality_test](https://en.wikipedia.org/wiki/Miller–Rabin_primality_test)."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 10,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"n: 26241 count_total: 52481 count_primes: 5248 ratio: 0.09999809454850327\n",
|
||
"26241\n",
|
||
"For 561 sicp says True and real test says False.\n",
|
||
"For 1105 sicp says True and real test says False.\n",
|
||
"For 1729 sicp says True and real test says False.\n",
|
||
"For 2465 sicp says True and real test says False.\n",
|
||
"For 2821 sicp says True and real test says False.\n",
|
||
"For 6601 sicp says True and real test says False.\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"carmichael_number = [561, 1105, 1729, 2465, 2821, 6601]\n",
|
||
"\n",
|
||
"def little_fermat(a, p): # SICP uses this\n",
|
||
" return (a ** p) % p == a % p\n",
|
||
"\n",
|
||
"def little_fermat_simple(a, p): # Wikipedia uses this\n",
|
||
" assert(a < p)\n",
|
||
" return (a ** (p - 1)) % p == 1\n",
|
||
"\n",
|
||
"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",
|
||
"def fermat_test(n):\n",
|
||
" ps = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 39]\n",
|
||
" if n <= ps[-1]:\n",
|
||
" return n in ps\n",
|
||
" for a in ps:\n",
|
||
" #if not little_fermat_simple(a, n):\n",
|
||
" if expmod(a, n - 1, n) != 1:\n",
|
||
" #print(\"Fermat witness for {} is {}.\".format(n + 1, a))\n",
|
||
" return False\n",
|
||
" return True\n",
|
||
"\n",
|
||
"def is_prime(n):\n",
|
||
" return fermat_test(n)\n",
|
||
"\n",
|
||
"print(get_solution())\n",
|
||
"\n",
|
||
"for n in carmichael_number:\n",
|
||
" r_sicp = fermat_test_sicp(n)\n",
|
||
" r = fermat_test(n)\n",
|
||
" print(\"For {} sicp says {} and real test says {}.\".format(n, r_sicp, r))\n"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {
|
||
"collapsed": true
|
||
},
|
||
"outputs": [],
|
||
"source": []
|
||
}
|
||
],
|
||
"metadata": {
|
||
"completion_date": "Sat, 29 Dec 2018, 09: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"
|
||
},
|
||
"tags": [
|
||
"prime",
|
||
"fermat",
|
||
"spiral",
|
||
"diagonal",
|
||
"todo"
|
||
]
|
||
},
|
||
"nbformat": 4,
|
||
"nbformat_minor": 2
|
||
}
|