Demonstrate the "prime conspiracy"

prime_conspiracy.py

'''Demonstrate the "prime number conspiracy" as described at
   https://www.quantamagazine.org/20160313-mathematicians-discover-prime-conspiracy/

"Among the first billion prime numbers, for instance, a prime ending in
9 is almost 65 percent more likely to be followed by a prime ending in 1
than another prime ending in 9. In a paper posted online today, Kannan
Soundararajan and Robert Lemke Oliver of Stanford University present
both numerical and theoretical evidence that prime numbers repel other
would-be primes that end in the same digit, and have varied
predilections for being followed by primes ending in the other possibl
final digits. "

'''

from collections import Counter
from random import randrange
from pprint import pprint

def is_prime(n, k=30):
    # http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
    if n <= 3:
        return n == 2 or n == 3
    neg_one = n - 1

    # write n-1 as 2^s*d where d is odd
    s, d = 0, neg_one
    while not d & 1:
        s, d = s+1, d>>1
    assert 2 ** s * d == neg_one and d & 1

    for i in range(k):
        a = randrange(2, neg_one)
        x = pow(a, d, n)
        if x in (1, neg_one):
            continue
        for r in range(1, s):
            x = x ** 2 % n
            if x == 1:
                return False
            if x == neg_one:
                break
        else:
            return False
    return True

def rand_prime(n):
    p = 1
    while not is_prime(p):
        p = randrange(n)
    return p

def next_prime(p):
    while True:
        p += 2
        if is_prime(p):
            return p

def check_conspiracy(trials, digits):
    limit = 10 ** digits
    pair_counts = Counter()
    for i in range(trials):
        p = rand_prime(limit)
        np = next_prime(p)
        pair = (p % 10, np % 10)
        pair_counts[pair] += 1
    return pair_counts

if __name__ == '__main__':
    pprint(check_conspiracy(trials=10000, digits=9).most_common())