2 reword

1

# Python 2.7, ~10^17

11111111171111111


Found 5111111111111 (13 digits) in 3 seconds and this 17 digit supreme prime in 3 minutes. I'll guess that the target machine could run this and get a 19 digit supreme prime in less than an hour. This approach does not scale well because it keeps primes up to half the number of target digits in memory. 17 digit search requires storing an array of 100M booleans. 19 digits would require a 1B element array, and memory would be exhausted before getting to 23 digits. Runtime probably would be, too.

Primality test approaches that don't involve a ridiculously large array of divisor primes will fare much better.

#!/usr/bin/env python
import math
import numpy as np
import sys

max_digits = int(sys.argv[1])
max_num = 10**max_digits

print "largest supreme prime of " + str(max_digits) + " or fewer digits"

def sum_product_digits(num):
mul = 1
while num:
add, mul, num = add + num % 10, mul * (num % 10), num / 10

def primesfrom2to(n):
# http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Input n>=6, Returns a array of primes, 2 <= p < n """
sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
sieve[0] = False
for i in xrange(int(n**0.5)/3+1):
if sieve[i]:
k=3*i+1|1
sieve[      ((k*k)/3)      ::2*k] = False
sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]

def checkprime(n):
for divisor in primes:
if (divisor>math.sqrt(n)):
break
if n%divisor==0:
return False
return True

# make an array of all primes we need to check as divisors of our max_num
primes = primesfrom2to(math.sqrt(max_num))
# only consider digit counts that are prime
for num_digits in primes:
if num_digits > max_digits:
break
for ones_on_right in range(0,num_digits):
for mid_prime in ['3','5','7']:
# assemble a number of the form /1*[357]1*/
candidate = int('1'*(num_digits-ones_on_right-1)+mid_prime+'1'*ones_on_right)
# check for primeness of digit sum first digit product first
add, mul = sum_product_digits(candidate)
if add in primes and mul in primes:
# check for primality next
if checkprime(candidate):
# supreme prime!
print candidate