# Find the largest prime whose length, sum and product is prime

The number 113 is the first prime whose length 3 is prime, digital sum 5 = 1 + 1 + 3 is prime, and digital product 3 = 1 * 1 * 3 is prime.

A prime that has these 3 properties will be called supremely prime. The primes 11117 and 1111151 are other examples.

# Goal

Write a program that can find the largest supremely prime number possible in less than an hour on a decent modern personal computer (such as the preferred spec here).

You should not simply give us a large supreme prime. You need to show us your search process with code that actually works. You can build on your or other people's solutions but be sure to give them credit. We're kind of communally trying to find the largest supreme prime realizable on a normal computer in an hour.

# Scoring

The submission that finds the largest supreme prime wins. If it turns out that there are finitely many supreme primes then the first submission that generates the highest supreme prime wins.

(If you can mathematically prove that either there are or aren't infinitely many supreme primes I'll give you 200 bounty rep just because. :) )

# Details

• You may use any source to generate your primes (e.g. internet).

• You may use probabilistic prime testing methods.

• Everything is in base 10.

• Zero and one are NOT considered prime.

• Primes that contain 0 have a digital product of 0 so obviously they can't be supreme.

• To keep the page less cluttered put large (100+ digit) supreme primes in the form:

  {[number of 1's before the prime digit]}[prime digit]{[number of 1's after the prime digit]}


So 1111151 could be expressed as {5}5{1}.

• Can we start with a list of primes, or fetch a list from the internet, and spend the hour doing supremeness checks? Jul 31, 2014 at 4:07
• If you can start with the highest known supreme prime then this becomes a challenge for who can write a program that spends exactly an hour spanning the largest possible gap between two supreme primes. :( Jul 31, 2014 at 4:24
• Next to not containing a 0, any potential supreme prime must obviously be of the form 1^n[3|5|7]1^m, i.e. some 1s, any prime below 10 and some more 1s. There are more constraints you can put on it right away. Jul 31, 2014 at 4:35
• Ryan has started a related question on MSE as to the existence of infinitely many supreme primes. If you've any insights for that question, please weigh in! Jul 31, 2014 at 15:30
• I can easily show that there is not currently a proof of an infinite number of supreme primes (and that a substantial amount of work has gone towards it). According to michaelnielsen.org/polymath1/…, we know that primes come along infinitely with gaps as small as 246, but for a proof of infinite supreme primes, we need a gap of 2, 4, or 6 (corresponding to the primes with a 3, 5, or 7 somewhere in them). Aug 1, 2014 at 15:36

# Python 2.7 on PyPy, {2404}3{1596} (~10^4000)

11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111113111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111


Found this one about 50 minutes after starting from 4000. Therefore, I would estimate this is the upper limit of this code approach.

Change: I've noticed that some lengths seem to be more fruitful for generating this sort of prime than others, so I've decided to only check 50 random locations of the digit that isn't 1 instead of all possible locations, before moving on. I'm not completely sure this will improve performance, or that 50 is correct, but we'll see.

Possibilities list is generated based on the fact that for the products requirement to be fulfilled, the number must be all ones except for a prime. In addition, the prime can't be 2 because of the sum and length relationship, and the digital sum must not be divisible by three, giving the %3 requirements.

is_prime taken from http://codepad.org/KtXsydxK , written by @primo

Note: this is_prime function is actually a Baillie-PSW pseudoprime test, but there are no known counter-examples, so I'm not going to worry about the distinction.

#http://codepad.org/KtXsydxK
from my_math import is_prime
import time,random
LLIMIT=2748
time.clock()
start_time=time.time()
checked=0
while time.time()-start_time<3600:
small_primes = [a for a in range(LLIMIT,2*LLIMIT) if is_prime(a)]
leng,dig=(0,0)
for a in small_primes:
if a+2 in small_primes:
leng,dig=(a,3)
break
if a+4 in small_primes:
leng,dig=(a,5)
break
if a+6 in small_primes:
leng,dig=(a,7)
break
start=time.clock()
print leng,dig,time.clock(),checked
for loc in random.sample(range(leng),50):
checked+=1
if is_prime(int('1'*loc+str(dig)+'1'*(leng-loc-1))):
print leng-1,loc,dig,time.clock(),time.clock()-start, \
int('1'*loc+str(dig)+'1'*(leng-loc-1))
break
LLIMIT=leng+1

• I know nothing more than the link, unfortunately. I found the link here: codegolf.stackexchange.com/questions/10739/… First answer Jul 31, 2014 at 5:20
• Well, then. I'll credit you. Jul 31, 2014 at 5:35
• It's like an ASCII where's wally... Jul 31, 2014 at 12:57
• Maybe you should rename the function is_very_very_very_very_very_very_very_probably_prime()... Jul 31, 2014 at 12:58
• Mathmatica and Maple both use the same method, so it can't be all that bad. Jul 31, 2014 at 13:14

# PARI/GP, 4127 digits

(104127-1)/9 + 2*10515

This is a fairly straightforward search: check only prime digit lengths, then compute the possible primes to use, then iterate through all possibilities. I special-cased the common cases where there are 0 or 1 suitable prime digits to use.

supreme(lim,startAt=3)={
forprime(d=startAt,lim,
my(N=10^d\9, P=select(p->isprime(d+p),[1,2,4,6]), D, n=1);
if(#P==0, next);
if(#P==1,
for(i=0,d-1,
if (ispseudoprime(D=N+n*P[1]), print(D));
n*=10
);
next
);
D=vector(#P-1,i,P[i+1]-P[i]);
for(i=0,d-1,
forstep(k=N+n*P[1],N+n*P[#P],n*D,
if (ispseudoprime(k), print(k))
);
n*=10
)
)
};
supreme(4200, 4100)


This took 36 minutes to compute on one core of a fairly old machine. It would have no trouble finding such a prime over 5000 digits in an hour, I'm sure, but I'm also impatient.

A better solution would be to use any reasonable language to do everything but the innermost loop, then construct an abc file for primeform which is optimized for that particular sort of calculation. This should be able to push the calculation up to at least 10,000 digits.

Edit: I implemented the hybrid solution described above, but on my old machine I can't find the first term with >= 10,000 digits in less than an hour. Unless I run it on something faster I'll have to change to a less lofty target.

• How did you know to start at 4100? Jul 31, 2014 at 20:10
• @isaacg: I was just trying to be larger than the (incorrect) Mathematica solution, which was just over 4000. I just went to the next multiple of 100 as a 'nothing-up-my-sleeve' number. Actually it seems that this was an unfortunate starting place, since I had to go longer than I expected (and longer than Mathematica!) to find a prime. Jul 31, 2014 at 20:15
• No, actually, you were incredibly lucky. (4127,3) is the first pair after 4100, and by pure chance it happens to have a prime. A lot of pairs have no prime at all. Jul 31, 2014 at 20:23
• @isaacg: Maybe so. My heuristics are clearly off, since I would have expected probability ~80% to find a prime in a given pair: 1 - exp(-15/(4*log 10)), but they seem to be rarer than that, so they don't act like random {2, 3, 5}-smooth numbers of their size (unless I goofed the calculation). Jul 31, 2014 at 20:28
• @isaacg: In any case I'm working on the "better solution" I mentioned now: pushing the hard work to pfgw. It's searched the first 20 pairs above 10^10000 without finding anything but that only took ~15 minutes. Jul 31, 2014 at 20:31

# Perl, 15101 digits, {83}7{15017}, 8 minutes. Max found: 72227 digits

Using my module Math::Prime::Util and its GMP back end. It has a number of compositeness tests including is_prob_prime() with an ES BPSW test (slightly more stringent than Pari's ispseudoprime), is_prime() which adds one random-base M-R, and is_provable_prime() which will run BLS75 T5 or ECPP. At these sizes and types, doing a proof is going to take a long time. I threw in another M-R test in the pedantic verifier sub. Times on a Core2 E7500 which is definitely not my fastest computer (it takes 2.5 minutes on my i7-4770K).

As Tim S. points out, we could keep searching for larger values, up to the point where a single test takes an hour. At ~15000 digits on this E7500 it takes about 26s for an M-R test and 2 minutes for the full is_prime (trial division plus base-2 M-R plus ES Lucas plus one random-base M-R). My i7-4770K is over 3x faster. I tried a few sizes, mainly seeing how it did on other people's results. I tried 8k, 20k, and 16k, killing each after ~5 minutes. I then tried 15k in progression for ~10m each and got lucky on the 4th one.

OpenPFGW's PRP tests are certainly faster once past 4000 or so digits, and much faster indeed in the 50k+ range. Its test has a fair amount of false positives however, which makes it a great pre-test but one would still like to verify the results with something else.

This could be parallelized with perl threads or using MCE similar to the parallel Fibonacci prime finder examples in the module.

Timing and results on an idle i7-4770K using single core:

• input 3000, 16 seconds, 3019 digits, {318}5{2700}
• input 4000, 47 seconds, 4001 digits, {393}7{3607}
• input 4100, 5 seconds, 4127 digits, {29}7{4097}
• input 6217, 5 seconds, 6217 digits, {23}5{6193}
• input 6500, 5 minutes, 6547 digits, {598}5{5948}
• input 7000, 15 minutes, 7013 digits, {2411}7{4601}
• input 9000, 11 minutes, 9001 digits, {952}7{8048}
• input 12000, 10 minutes, 12007 digits, {652}5{11354}
• input 15100, 2.5 minutes, 15101 digits, {83}7{15017}
• input 24600, 47 minutes, 24671 digits, {621}7{24049}
• input 32060, 18 minutes, 32063 digits, {83}7{31979}
• input 57000, 39 minutes, 57037 digits, {112}5{56924}
• input 72225, 42 minutes, 72227 digits, {16}3{72210}

For the 32k digit result, I started 6 scripts running at the same time each with successive arguments starting at 32000. After 26.5 minutes one finished with the 32063 digit result shown. For 57k I let successive scripts run 6 at a time for an hour in input increments of 500 until the 57k result returned in 57 minutes. The 72k digit result was found by doing successive primes from 70k up, so definitely not found in an hour (though once you know where to start, it is).

The script:

#!/usr/bin/env perl
use warnings;
use strict;
use Math::Prime::Util qw/:all/;
use Math::Prime::Util::GMP;  # Just to ensure it is used.

my $l = shift || 1000;$l--;

while (1) {
$l = next_prime($l);
my @D = grep { is_prime($l-1 +$_) } (3,5,7);
next unless scalar @D > 0;
for my $s (0 ..$l-1) {
my $e =$l-$s-1; warn " checking$l $s\n" unless$s % 100;
for my $d (@D) { my$n = "1"x$s .$d . "1"x$e; die unless length($n) == $l; verify_supreme($n,$s,$d,$e) if is_prime($n);  # ES BPSW + 1 rand-base M-R
}
}
}
sub verify_supreme {  # Be pedantic and verify the result
my($n,$s,$d,$e) = @_;
die "Incorrect length" unless is_prime(length($n)); die "Incorrect sum" unless is_prime(vecsum(split(//,$n)));
my $prod = 1;$prod *= $_ for split(//,$n);
die "Incorrect product" unless is_prime($prod); die "n is not a prime!" unless miller_rabin_random($n,1);  # One more M-R test
die "{$s}$d {$e}\n"; }  • +1 for introducing me to this library! Timings on my machine to iterate primes less than 10^7 (compared to CPython with gmpy2, and PyPy with my_math): codepad.org/aSzc0esT Aug 2, 2014 at 6:49 • Glad you like it! There are other ways including forprimes { ...do stuff... } 1e7; which is 10x or more faster (kudos to Pari/GP for lots of great ideas). I always appreciate feedback so let me know if something isn't working the way you'd like. Aug 2, 2014 at 17:44 # Mathematica 3181 digits Update: There were some serious mistakes in my first submission. I was able to devote some time to checking the results for this one. The output is formatted as a list of digits. That makes for easy checking of each of the conditions. f[primeDigitLength_]:= Module[{id=ConstantArray[1,primeDigitLength-1]}, permutations=Reverse@Sort@Flatten[Table[Insert[id,#,pos],{pos,primeDigitLength}]&/@{3,5,7},1]; Flatten[Select[permutations,PrimeQ[FromDigits[#]]\[And]PrimeQ[Plus@@#]&,1],1]]  Example This was my first test, a search for a solution with 3181 digits. It found the first case in 26 seconds. Let's go through the reasoning. Then we'll step into the program itself. Let's start, as I did, "What is the 450th prime?" Can we find a solution with that many digits (3181)? primeDigits = Prime[450]  3181 The number is found by joining the digits. number = FromDigits[digits];  But rather than display it, we can ask instead what the digits are and where they lie. DigitCount[number]  {3180, 0, 0, 0, 0, 0, 1, 0, 0, 0} This means that there were 3180 instances of the digit 1, and a single instance of the digit 7. At what position is the digit 7? Position[digits, 7][[1, 1]]  142 So the digit 7 is the 142nd digit. All the others are 1's. Of course, the product of the digits must be a prime, namely 7. digitProduct = Times @@ digits  7 And the sum of the digits is also a prime. digitSum = Plus @@ digits PrimeQ[digitSum]  3187 True And we know that the number of digits is a prime. Remember, we selected the 450th prime, namely 3118. So all the conditions have been met. • If I am not mistaken, its sum is 4009 which is not prime. – gerw Jul 31, 2014 at 13:20 • One thing : shouldn't it be the sum of all the digits that is prime and not the number of digits? In your case you would have to test 4002 * 1 + 7 = 4009 and not 4003 according to the spec. Jul 31, 2014 at 13:20 • @Johnride: Both should be prime. – gerw Jul 31, 2014 at 13:25 • @gerw Is right. The number of digits AND the sum of the digits AND the product of digits all have to be prime. Jul 31, 2014 at 16:04 • You all were correct. In my earlier submission I forgot to check the digit sum for primality. This is now done by seeing whether one of the following (it doesn't matter which) is prime: digit length+2, digit length _4, or digit Length +6. Aug 1, 2014 at 7:57 # Python 2.7, 6217 digits: {23}5{6193} 6 mins 51 secs I was working on my own version and was disappointed to see that @issacg had beaten me to the punch with a very similar approach, albeit with is_(very_probably)_prime(). However, I see that I have some significant differences that result in a better answer in less time (when I also use is_prime). To make this clear, when also starting from 4000, I arrive at a better 4001 digit answer ({393}7{3607}) in only 26 mins, 37 seconds using the standard Python interpreter (also at version 2.7), not the PyPy version. Also, I'm not 'spot' checking the numbers. All of the candidates are checked. Here's the primary improvements: 1. Use a prime generator (https://stackoverflow.com/questions/567222/simple-prime-generator-in-python/568618#568618) to create a list of primes to check against and (his version of "small primes") and for generating eligible number lengths. 2. We want to spend our time looking for the largest supreme prime of a given length, not the smallest, so, I construct the largest possible numbers first for checking, not the smallest. Then, once one is found, we can immediately move on to the next length. # EDIT: Now with multiprocessing This is a significant change to previous versions. Before, I noticed that my 8-core machine was hardly working, so, I decided to try my hand at multiprocessing in Python (first time). The results are very nice! In this version, 7 children processes are spawned, which grab a 'task' off a queue of potential possibilities (num_length + eligible digits). They churn through trying different [7,5,3] positions until it finds one. If it does, informs the master process of the new longest length that has been found. If children are working on a num_length that is shorter, they just bail, and go get the next length. I started this run with 6000, and it is still running, but so far, I'm very pleased with the results. The program doesn't yet stop correctly, but not a huge deal to me. Now the code: #!/usr/bin/env python from __future__ import print_function import sys from multiprocessing import Pool, cpu_count, Value from datetime import datetime, timedelta # is_prime() et al from: http://codepad.org/KtXsydxK - omitted for brevity # gen_primes() from: https://stackoverflow.com/questions/567222/simple-prime-generator-in-python/568618#568618 - ommitted for brevity from external_sources import is_prime, gen_primes def gen_tasks(start_length): """ A generator that produces a stream of eligible number lengths and digits """ for num_length in gen_primes(): if num_length < start_length: continue ns = [ n for n in [7,5,3] if num_length + n - 1 in prime_list ] if ns: yield (num_length, ns) def hunt(num_length, ns): """ Given the num_length and list of eligible digits to try, build combinations to try, and try them. """ if datetime.now() > end_time or num_length <= largest.value: return print('Tasked with: {0}, {1}'.format(num_length, ns)) sys.stdout.flush() template = list('1' * num_length) for offset in range(num_length): for n in ns: if datetime.now() > end_time or num_length <= largest.value: return num_list = template[:] num_list[offset] = str(n) num = int(''.join(num_list)) if is_prime(num): elapsed = datetime.now() - start_time largest.value = num_length print('\n{0} - "{1}"\a'.format(elapsed, num)) if __name__ == '__main__': start_time = datetime.now() end_time = start_time + timedelta(seconds=3600) print('Starting @ {0}, will stop @ {1}'.format(start_time, end_time)) start_length = int(sys.argv[1]) # # Just create a list of primes for checking. Up to 20006 will cover the first # 20,000 digits of solutions # prime_list = [] for prime in gen_primes(): prime_list.append(prime) if prime > 20006: break; print('prime_list is primed.') largest = Value('d', 0) task_generator = gen_tasks(start_length) cores = cpu_count() print('Number of cores: {0}'.format(cores)) # # We reduce the number of cores by 1 because __main__ is another process # pool = Pool(processes=cores - 1) while datetime.now() < end_time: pool.apply_async(hunt, next(task_generator))  • it would read more cleanly if you represent the codepad link as a [broken, if necessary] import Aug 1, 2014 at 3:05 • I think that would be confusing, as the code at the other end isn't really importable like that. Aug 1, 2014 at 3:08 • use isaacg's syntax. comment the URL, then import from a nonexistent package (my_math, in his case) Aug 1, 2014 at 3:14 • Actually, I generate the numbers from largest to smallest, too. I don't think our code differences are very significant. I'd expect the difference to lie more in our computers than our code. Nevertheless, well done, and welcome to the site. Aug 1, 2014 at 3:37 • my_math can also be used to generate a list of primes, à la while prime < 20006: prime = next_prime(prime). Seems to be about 3 times as fast as gen_primes, and far more memory efficient. Aug 1, 2014 at 5:51 ## C, GMP - {7224}5{564} = 7789 Kudos to @issacg and all of you guys for the inspirations and tricks. And also the masterful question asker @Calvin's Hobbies for this question. Compile: gcc -I/usr/local/include -o p_out p.c -pthread -L/usr/local/lib -lgmp  If you feel like donating your computation power or curious of the performance, feel free copy the code and compile. ;) You will need GMP installed. #include<stdio.h> #include<stdlib.h> #include<sys/time.h> #include<gmp.h> #include<pthread.h> #define THREAD_COUNT 1 #define MAX_DIGITS 7800 #define MIN_DIGITS 1000 static void huntSupremePrime(int startIndex) { char digits[MAX_DIGITS + 1]; for (int i = 0; i < MAX_DIGITS; digits[i++] = '1'); digits[MAX_DIGITS] = '\0'; mpz_t testPrime, digitSum, digitCount, increment; for (int j = 0; j < MAX_DIGITS - startIndex; digits[j++] = '0'); int step = THREAD_COUNT * 2; for (int i = startIndex, l = MAX_DIGITS - startIndex; i > MIN_DIGITS - 1; i -= step, l += step) { fprintf(stderr, "Testing for %d digits.\n", i); mpz_init_set_ui(digitCount, i); if (mpz_millerrabin(digitCount, 10)) { for (int j = 3; j < 8; j += 2) { mpz_init_set_ui(digitSum, i - 1 + j); if (mpz_millerrabin(digitSum, 10)) { gmp_printf("%Zd \n", digitSum); digits[MAX_DIGITS - 1] = j + 48; mpz_init_set_str(testPrime, digits, 10); mpz_init_set_ui(increment, (j - 1) * 99); for (int k = 0; k < i/20; ++k) { if (mpz_millerrabin(testPrime, 25)) { i = 0; j = 9; k = l; gmp_printf("%Zd\n", testPrime); break; } mpz_add(testPrime, testPrime, increment); mpz_mul_ui(increment, increment, 100); fprintf(stderr, "TICK %d\n", k); } } } } for (int j = 0; j < step; digits[l + j++] = '0'); } } static void *huntSupremePrimeThread(void *p) { int* startIndex = (int*) p; huntSupremePrime(*startIndex); pthread_exit(NULL); } int main(int argc, char *argv[]) { int startIndexes[THREAD_COUNT]; pthread_t threads[THREAD_COUNT]; int startIndex = MAX_DIGITS; for (int i = 0; i < THREAD_COUNT; ++i) { for (;startIndex % 2 == 0; --startIndex); startIndexes[i] = startIndex; int rc = pthread_create(&threads[i], NULL, huntSupremePrimeThread, (void*)&startIndexes[i]); if (rc) { printf("ERROR; return code from pthread_create() is %d\n", rc); exit(-1); } --startIndex; } for (int i = 0; i < THREAD_COUNT; ++i) { void * status; int rc = pthread_join(threads[i], &status); if (rc) { printf("ERROR: return code from pthread_join() is %d\n", rc); exit(-1); } } pthread_exit(NULL); return 0; }  ## PFGW, 6067 digits, {5956}7{110} Run PFGW with the following input file and -f100 to prefactor the numbers. In about 2-3 CPU minutes on my computer (i5 Haswell), it finds the PRP (10^(6073-6)-1)/9+6*10^110, or {5956}7{110}. I chose 6000 digits as the starting point as a nothing-up-my-sleeve number that's a little higher than all previous submissions. ABC2$a-$b & (10^($a-$b)-1)/9+$b*10^\$c
a: primes from 6000 to 6200
b: in { 2 4 6 }
c: from 0 to 5990


Based on how quickly I was able to find this one, I could easily bump up the number of digits and still find a PRP within an hour. With how the rules are written, I might even just find the size where my CPU, running on all 4 cores, can finish one PRP test in an hour, take a long time to find a PRP, and have my "search" consist solely of the one PRP.

P.S. In some senses, this isn't a "code" solution because I didn't write anything besides the input file...but then, many one-line Mathematica solutions to mathematical problems could be described in the same way, as could using a library that does the hard work for you. In reality, I think it's hard to draw a good line between the two. If you like, I could write a script that creates the PFGW input file and calls PFGW. The script could even search in parallel, to use all 4 cores and speed up the search by ~4 times (on my CPU).

P.P.S. I think LLR can do the PRP tests for these numbers, and I'd expect it to be far faster than PFGW. A dedicated sieving program could also be better at factoring these numbers than PFGW's one-at-a-time. If you combined these, I'm sure you could push the bounds much higher than current solutions.

# Python 2.7, 17-19 digits

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
if add in primes and mul in primes:
# check for primality next
if checkprime(candidate):
# supreme prime!
print candidate


# Mathematica 4211 4259 digits

With the number: {1042}7{3168} {388}3{3870}

Which was generated by the following code:

TimeConstrained[
Do[
p = Prime[n];
curlargest = Catch[
If[PrimeQ[p + 6],
list = ConstantArray[1, p];
Do[
temp = FromDigits[ReplacePart[list, i -> 7]];
If[PrimeQ[temp],
Throw[temp]
], {i, p}]
];

If[PrimeQ[p + 4],
list = ConstantArray[1, p];
Do[
temp = FromDigits[ReplacePart[list, i -> 5]];
If[PrimeQ[temp],
Throw[temp]
], {i, p}]
];
If[PrimeQ[p + 2],
list = ConstantArray[1, p];
Do[
temp = FromDigits[ReplacePart[list, i -> 3]];
If[PrimeQ[temp],
Throw[temp]
], {i, p}]
];
Throw[curlargest];
]

, {n, 565, 10000}]
, 60*60]


The throws cause it to stop testing for other numbers with the same digits as the one currently found. since it begins testing at the most significant digit this will mean that it always returns the largest number unless the number of digits is a member of a prime triplet.

Simply started testing just below the value of one of the preceding answers :)

Once it finishes, the number is stored in the variable curlargest

# JavaScript, 3019 digits, {2,273}5{745}

This uses the MillerRabin test included in BigInteger.js by Tom Wu.

Starting from 0 => 2,046 digits = {1799}7{263} in one hour.

Starting from 3000 => 3,019 digits = {2,273}5{745} in one hour, less 3 seconds.

When it started from 0, the program skipped ahead and began searching again at a length of 1.5X the length of the last s-prime found. Then when I saw how fast it was running I guessed it would find one starting at 3000 in one hour - which it did with only 3 seconds to spare.

You can try it here: http://goo.gl/t3TmTk
(Set to calculate all s-primes, or skip ahead.)

The program works by constructing strings of all "1"s, but with one "3", "5", or "7". I added a quick check in the IsStrPrime function to reject numbers ending in "5".

if (IsIntPrime(length)) {

var do3s = IsIntPrime(length + 2);
var do5s = IsIntPrime(length + 4);
var do7s = IsIntPrime(length + 6);

if (do3s || do5s || do7s) {

// loop through length of number
var before, digit, after;

for (var after = 0; after <= length - 1; after++) {

before = length - after - 1;
beforeStr = Ones(before);
afterStr = Ones(after);

if (do3s && IsStrPrime(beforeStr + (digit = "3") + afterStr)) { RecordAnswer(); if (brk) break; }
if (AreDone()) break;

if (do5s && IsStrPrime(beforeStr + (digit = "5") + afterStr)) { RecordAnswer(); if (brk) break; }
if (AreDone()) break;

if (do7s && IsStrPrime(beforeStr + (digit = "7") + afterStr)) { RecordAnswer(); if (brk) break; }
if (AreDone()) break;

if (after % 10 == 0) document.title = "b=" + bestLength + ", testing=" + length + "-" + after;
}
}
}


This was pretty fun. Reminds me of a puzzle I did many years back to calculate what's called a digit removed prime. This is a prime number that if you remove any digit, then the remaining number is still prime. For example 1037 is a digit removed prime because 1037, 037, 137, 107, and 103 are prime. I found one 84 digits long, and the longest I know of is 332 digits long. I'm sure we could find one much longer with the techniques used for this puzzle. (But choosing the trial numbers is a little bit trickier - maybe?)

• RE: digit removed prime, we've had that one here. 332 digits would have won, too. Aug 28, 2014 at 10:23

# Axiom, 3019 digits {318}5{2700}

)time on

-- Return true if n is pseudo prime else return false
-- **Can Fail**
prime1(n:PI):Boolean==
n=1=>false
n<4=>true
i:=5;sq:=sqrt(1.*n)
repeat
if i>sq or i>50000 then break
if n rem i=0       then return false
i:=i+2
if i~=50001        then return true
--output("i")
if powmod(3,n,n)=3 then return true
--output("e")
false

-- input  'n': must be n>1 prime
-- output [0] if not find any number, else return
-- [n,a,b,c,z] 'n' digits of solution,
-- 'a' number of '1', 'b' central digit, 'b' number of final digit '1'
-- 'z' the number found
g(n:PI):List NNI==
x:=b:=z:=1
for i in 1..n-1 repeat
z:=z*10+1
b:=b*10
repeat
--output b
k:=0    -- 3 5 7 <-> 2 4 6
for i in [2,4,6] repeat
~prime?(n+i)=>0 --somma
k:=k+1
t:=z+b*i
if prime1(t) then return [n,x-1,i+1,n-x,t]
--if x=1 then output ["k=", k]
if k=0  then break
x:=x+1
b:=b quo 10
if b<=0 then break
[0]

-- start from number of digits n
-- and return g(i) with i prime i>=n
find(n:PI):List NNI==
i:=n
if i rem 2=0 then i:=i+1
repeat
if prime?(i) then --solo le lunghezze prime sono accettate
output i
a:=g(i)
if a~=[0] then return a
i:=i+2


result from the start value 3000 in 529 sec

(4) -> find(3000)
3001
3011
3019

(4)
[3019, 318, 5, 2700, Omissis]
Type: List NonNegativeInteger
Time: 0.02 (IN) + 525.50 (EV) + 0.02 (OT) + 3.53 (GC) = 529.07 sec