Consecutive primes last digit blog article

scipython · Mar 14, 2018 · a4f08d7 · a4f08d7
1 parent b2fd58b
commit a4f08d7
Show file tree

Hide file tree

Showing 5 changed files with 173 additions and 0 deletions.
diff --git a/primes_last_digits/last_digits_hmap.py b/primes_last_digits/last_digits_hmap.py
@@ -0,0 +1,40 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.ticker import FuncFormatter
+
+# Create a heatmap summarizing the probabilities associated with
+# the last digits of consecutive prime numbers. More details
+# at https://scipython.com/blog/do-consecutive-primes-avoid-sharing-the-same-last-digit/
+# Christian Hill, March 2016.
+
+# First 10,000,000 primes
+digit_count = {1: {1: 446808, 3: 756072, 9: 526953, 7: 769924},
+               3: {1: 593196, 3: 422302, 9: 769915, 7: 714795},
+               9: {1: 820369, 3: 640076, 9: 446032, 7: 593275},
+               7: {1: 639384, 3: 681759, 9: 756852, 7: 422289}}
+last_digits = [1,3,7,9]
+
+hmap = np.empty((4,4))
+for i, d1 in enumerate(last_digits):
+    total = sum(digit_count[d1].values())
+    for j, d2 in enumerate(last_digits):
+        hmap[i,j] = digit_count[d1][d2] / total * 100
+
+fig = plt.figure()
+ax = fig.add_axes([0.1,0.3,0.8,0.6])
+im = ax.imshow(hmap, interpolation='nearest', cmap=plt.cm.YlOrRd, origin='lower')
+tick_labels = [str(d) for d in last_digits]
+ax.set_xticks(range(4))
+ax.set_xticklabels(tick_labels)
+ax.set_xlabel('Last digit of second prime')
+ax.set_yticks(range(4))
+ax.set_yticklabels(tick_labels)
+ax.set_ylabel('Last digit of first prime')
+
+cbar_axes = fig.add_axes([0.1,0.1,0.8,0.05])
+
+cbar = plt.colorbar(im, orientation='horizontal', cax=cbar_axes)
+cbar.ax.set_xlabel('Probability /%')
+
+plt.savefig('prime_digits_hmap.png')
+plt.show()
diff --git a/primes_last_digits/plot_prime_digits.py b/primes_last_digits/plot_prime_digits.py
@@ -0,0 +1,45 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Create a set of bar charts summarizing the probabilities associated with
+# the last digits of consecutive prime numbers. More details
+# at https://scipython.com/blog/do-consecutive-primes-avoid-sharing-the-same-last-digit/
+# Christian Hill, March 2016.
+
+# Dictionary of consecutive digit counts for the first 10,000,000 primes
+digit_count = {1: {1: 446808, 3: 756072, 9: 526953, 7: 769924},
+               3: {1: 593196, 3: 422302, 9: 769915, 7: 714795},
+               9: {1: 820369, 3: 640076, 9: 446032, 7: 593275},
+               7: {1: 639384, 3: 681759, 9: 756852, 7: 422289}}
+fig, ax = plt.subplots(nrows=2, ncols=2, facecolor='#dddddd')
+
+xticks = [0,1,2,3]
+last_digits = [1,3,7,9]
+for i, d in enumerate(last_digits):
+    ir, ic = i // 2, i % 2 
+    this_ax = ax[ir,ic]
+    this_ax.patch.set_alpha(1)
+    count = np.array([digit_count[d][j] for j in last_digits])
+    total = sum(count)
+    prob = count / total * 100
+    this_ax.bar(xticks, prob, align='center', color='maroon', ec='maroon',
+                alpha=0.7)
+    this_ax.set_title('Last digit of prime: {:d}'.format(d), fontsize=14)
+    this_ax.set_xticklabels(['{:d}'.format(j) for j in last_digits])
+    this_ax.set_xticks(xticks)
+    this_ax.set_yticks([0,10,20,30,40])
+    this_ax.set_ylim(0,35)
+    this_ax.set_yticks([])
+    for j, pr in enumerate(prob):
+        this_ax.annotate('{:.1f}%'.format(pr), xy=(j, pr-2), ha='center',
+                         va='top', color='w', fontsize=12)
+    this_ax.set_xlabel('Next prime ends in')
+    this_ax.set_frame_on(False)
+    this_ax.tick_params(axis='x', length=0)
+    this_ax.tick_params(axis='y', length=0)
+fig.subplots_adjust(wspace=0.2, hspace=0.7, left=0, bottom=0.1, right=1,
+                    top=0.95)
+plt.savefig('prime_digits.png')
+plt.show()
+
+
diff --git a/primes_last_digits/prime-last-digits.py b/primes_last_digits/prime-last-digits.py
@@ -0,0 +1,88 @@
+import time
+import math
+
+# Prime numbers greater than 5 end in 1, 3, 7 or 9. It seems that a prime
+# ending in one digit is less likely to be followed by another ending in the
+# same digit. Return a dictionary summarizing this distribution. More details
+# at https://scipython.com/blog/do-consecutive-primes-avoid-sharing-the-same-last-digit/
+# Christian Hill, March 2016.
+
+def approx_nth_prime(n):
+    """Return an upper bound for the value of the nth prime"""
+
+    return n * (math.log(n) + math.log(math.log(n)))
+
+nmax = 10000000
+pmax = approx_nth_prime(nmax)
+print('The {:d}th prime is approximately {:d}'.format(nmax,int(pmax)))
+N = int(math.sqrt(pmax)) + 1
+print('Our sieve will therefore contain primes up to', N)
+
+def primes_up_to(N):
+    """A generator yielding all primes less than N."""
+
+    yield 2
+    # Only consider odd numbers up to N, starting at 3
+    bsieve = [True] * ((N-1)//2)
+    for i,bp in enumerate(bsieve):
+        p = 2*i + 3
+        if bp:
+            yield p
+            # Mark off all multiples of p as composite
+            for m in range(i, (N-1)//2, p):
+                bsieve[m] = False
+
+gen_primes = primes_up_to(N)
+sieve = list(gen_primes)
+
+def is_prime(n, imax):
+    """Return True if n is prime, else return False.
+
+    imax is the maximum index in the sieve of potential prime factors that
+    needs to be considered; this should be the index of the first prime number
+    larger than the square root of n.
+
+    """
+    return not any(n % p == 0 for p in sieve[:imax])
+
+
+digit_count = {1: {1: 0, 3: 0, 7: 0, 9: 0},
+               3: {1: 0, 3: 0, 7: 0, 9: 0},
+               7: {1: 0, 3: 0, 7: 0, 9: 0},
+               9: {1: 0, 3: 0, 7: 0, 9: 0}}
+
+# nprimes is the number of prime numbers encountered
+nprimes = 0
+# the most recent prime number considered (we start with the first prime number
+# which ends with 1,3,7 or 9 and is followed by a number ending with one of
+# these digits, 7 since 2, 3 and 5 are somewhat special cases.
+last_prime = 7
+# The current prime number to consider, initially the one after 7 which is 11
+n = 11
+# The index of the maximum prime in our sieve we need to consider when testing
+# for primality: initially 2, since sieve[2] = 5 is the nearest prime larger
+# than sqrt(11). plim is this largest prime from the sieve.
+imax = 2
+plim = sieve[imax]
+start_time = time.time()
+
+while nprimes <= nmax:
+    # Output a progress indicator
+    if not nprimes % 1000:
+        print(nprimes)
+
+    if is_prime(n, imax):
+        # n is prime: update the dictionary of last digits
+        digit_count[last_prime % 10][n % 10] += 1
+        last_prime = n
+        nprimes += 1
+    # Move on to the next candidate (skip even numbers)
+    n += 2
+
+    # Update imax and plim if necessary
+    if math.sqrt(n) >= plim:
+        imax += 1
+        plim = sieve[imax]
+end_time = time.time()
+print(digit_count)
+print('Time taken: {:.2f} s'.format(end_time - start_time))
diff --git a/primes_last_digits/prime_digits.png b/primes_last_digits/prime_digits.png
diff --git a/primes_last_digits/prime_digits_hmap.png b/primes_last_digits/prime_digits_hmap.png