# This Python script will find all permutable primes having two or more distinct digits.

import sys
from numtheory import is_probable_prime    # Miller-Rabin probabilistic primality test
from itertools import permutations, combinations_with_replacement

# Read a list of digits in a given base (default = 10) into an integer.

def from_base(v, base = 10):
    n = 0
    for d in v:
        n = n*base + d
    return n

# Prove that a permutable prime cannot contain a given multiset of digits
# by checking all congruence classes modulo m. The proof is indented to
# the specified level.

def prove_case(digits, m, level=1):
    mod_count = 0
    witness = [0] * m
    indent = "    " * level
    for n in map(from_base, permutations(digits)):
        r = n % m
        if witness[r] == 0:
            witness[r] = n
            mod_count += 1
        if mod_count == m:
            for r in range(m):
                print "%s%d = %d (mod %d)" % (indent, witness[r], r, m)
            print
            break
    else:
        sys.exit("Proof failed.")

# Determine if a number with a given list of digits is a permutable prime.

def is_permutable_prime(digits):
    for n in map(from_base, permutations(digits)):
        if not(is_probable_prime(n)):
            return False
    return True

print "We begin with a brute force search for permutable primes less than 10^6."
print "We only list the permutable primes whose digits are in ascending order."
print "The one-digit primes are omitted.\n"

print "    ",
for num_digits in range(2,7):
    for digits in combinations_with_replacement([1,3,7,9], num_digits):
        if is_permutable_prime(digits):
            print from_base(digits),            
        
print "\n\nWe will prove that no other permutable primes exist, except for repunits."
print "The proof uses a lemma: If N and m are integers greater than 1 such that"
print "the digit permutations of N cover all congruence classes modulo m, then"
print "no permutable prime can contain all of the digits of N (as a multiset)."
print "By choosing a finite number of pairs (N,m) we can rule out all numbers"
print "having six or more digits, except for repunits.\n"

print "Case 1: N has four distinct digits (1,3,7,9). We use m = 7. \n"

prove_case((1,3,7,9), 7)

print "Case 2: N has three distinct digits.\n"

print "    Case 2A: N contains aabbc.\n"

for a in (1,3,7):
    for b in (3,7,9):
        if b > a:
            for c in (1,3,7,9):
                if c != a and c != b:
                    prove_case((a,a,b,b,c), 7, 2)

print "    Case 2B: N contains aaaaabc.\n"

for a in (1,3,7,9):
    for b in (1,3,7,9):
        if b != a:
            for c in (3,7,9):
                if c > b and c != a:
                    prove_case((a,a,a,a,a,b,c), 7, 2)

print "Case 3: N has two distinct digits.\n"

print "    Case 3A: N contains aaabb.\n"

for a in (1,3,7,9):
    for b in (1,3,7,9):
        if b != a:
            prove_case((a,a,a,b,b), 7, 2)

print "    Case 3B: N is a near-repdigit (one digit is different from the others).\n"

print "         We show that a near-repunit cannot be prime if it has 17 or more digits."
print "         Suppose that N = aaa...ab, where a is repeated 16 or more times."
print "         Then the numbers N - (b-a) + (b-a)*10^k are permutations of N"
print "         for k = 0, 1, 2, ..., 16. These numbers cover all congruence classes"
print "         modulo 17, since 10 is a primitive root. Therefore, 17 divides one of"
print "         these numbers.\n"

print "         We check (using the Miller-Rabin primality test) that the following"
print "         numbers are NOT prime.\n"

indent = " "*11
repunit = 111111
for digits in range(7, 17):
    repunit = repunit * 10 + 1
    for a in (1,3,7,9):
        for b in (1,3,7,9):
            if b != a:
                for k in range(digits):
                    n = a*repunit + (b-a)*10**k
                    if not(is_probable_prime(n)):
                        print indent, n
                        break
                else:
                    sys.exit("Proof failed.")
    print
print "QED."