felixm/euler
1
0
Fork 0
Code Pull Requests Releases Activity

Figured out problem with Fermat test in problem 58.

Browse Source
This commit is contained in:
Felix Martin
2019-02-02 12:22:43 -05:00
parent fe73c53c46
commit 45c7e862d3
3 changed files with 160 additions and 57 deletions
Download Patch File Download Diff File Expand all files Collapse all files

118
ipython/EulerProblem058.ipynb
View File

@@ -33,7 +33,9 @@
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def get_last_corner_value(side_length):\n",
@@ -81,7 +83,7 @@
" return (base * expmod(base, exp - 1, m) % m)\n",
" \n",
"\n",
"def fermat_test(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",
@@ -145,10 +147,10 @@
" 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))"
"assert(fermat_test_sicp(3))\n",
"assert(fermat_test_sicp(107))\n",
"assert(fermat_test_sicp(108) == False)\n",
"assert(fermat_test_sicp(109))"
]
},
{
@@ -210,44 +212,33 @@
},
{
"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,
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n: 26641 count_total: 53281 count_primes: 5328 ratio: 0.09999812315834913\n"
"n: 26641 count_total: 53281 count_primes: 5328 ratio: 0.09999812315834913\n",
"26641\n",
"Expected n is 26241.\n"
]
}
],
"source": [
"s = get_solution()"
"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.\")"
]
},
{
@@ -261,6 +252,65 @@
"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,