From 45c7e862d366b9efd63e53dbe3b6f6604d766fee Mon Sep 17 00:00:00 2001 From: Felix Martin Date: Sat, 2 Feb 2019 12:22:43 -0500 Subject: [PATCH] Figured out problem with Fermat test in problem 58. --- ipython/EulerProblem058.ipynb | 118 +++++++++++++++++++++--------- ipython/html/EulerProblem058.html | 97 ++++++++++++++++++------ ipython/html/EulerProblem068.html | 2 +- 3 files changed, 160 insertions(+), 57 deletions(-) diff --git a/ipython/EulerProblem058.ipynb b/ipython/EulerProblem058.ipynb index c68c6d0..322c1b5 100644 --- a/ipython/EulerProblem058.ipynb +++ b/ipython/EulerProblem058.ipynb @@ -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, diff --git a/ipython/html/EulerProblem058.html b/ipython/html/EulerProblem058.html index 063afb0..9b96e8b 100644 --- a/ipython/html/EulerProblem058.html +++ b/ipython/html/EulerProblem058.html @@ -11846,7 +11846,7 @@ div#notebook { return (base * expmod(base, exp - 1, m) % m) -def fermat_test(n): +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 @@ -11910,10 +11910,10 @@ div#notebook { return False return trial_division(n)[0] == n -assert(fermat_test(3)) -assert(fermat_test(107)) -assert(fermat_test(108) == False) -assert(fermat_test(109)) +assert(fermat_test_sicp(3)) +assert(fermat_test_sicp(107)) +assert(fermat_test_sicp(108) == False) +assert(fermat_test_sicp(109)) @@ -11989,7 +11989,7 @@ div#notebook {
-
In [5]:
+
In [12]:
import random
@@ -11997,27 +11997,14 @@ div#notebook {
 def is_prime(n):
     if n == 1:
         return False
-    for _ in range(100):
+    for _ in range(10):
         a = n - random.randint(1, n - 1)
         if not expmod(a, n, n) == a:
             return False
     return True
 
-assert(fermat_test(3))
-assert(fermat_test(107))
-assert(fermat_test(108) == False)
-assert(fermat_test(109))
-
-
-
-
-
-
-
-
In [6]:
-
-
-
s = get_solution()
+print(get_solution())
+print("Expected n is 26241.")
@@ -12028,6 +12015,8 @@ div#notebook { 
n: 26641 count_total: 53281 count_primes: 5328 ratio: 0.09999812315834913
+26641
+Expected n is 26241.
@@ -12046,6 +12035,70 @@ div#notebook {
+
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 [ ]:
diff --git a/ipython/html/EulerProblem068.html b/ipython/html/EulerProblem068.html index 3387d58..13da685 100644 --- a/ipython/html/EulerProblem068.html +++ b/ipython/html/EulerProblem068.html @@ -11770,7 +11770,7 @@ div#notebook {
-

Problem Name (Euler Problem xxx)

Back to overview.

+

Magic 5-gon ring (Euler Problem 68)

Back to overview.