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.


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.


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. :) )


  • 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}.

  • \$\begingroup\$ Can we start with a list of primes, or fetch a list from the internet, and spend the hour doing supremeness checks? \$\endgroup\$
    – Sparr
    Jul 31, 2014 at 4:07
  • 3
    \$\begingroup\$ 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. :( \$\endgroup\$
    – Sparr
    Jul 31, 2014 at 4:24
  • 8
    \$\begingroup\$ 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. \$\endgroup\$
    – Ingo Bürk
    Jul 31, 2014 at 4:35
  • 3
    \$\begingroup\$ 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! \$\endgroup\$ Jul 31, 2014 at 15:30
  • 1
    \$\begingroup\$ 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). \$\endgroup\$
    – Tim S.
    Aug 1, 2014 at 15:36

12 Answers 12


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


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.

from my_math import is_prime
import time,random
while time.time()-start_time<3600:
    small_primes = [a for a in range(LLIMIT,2*LLIMIT) if is_prime(a)]
    for a in small_primes:
        if a+2 in small_primes:
        if a+4 in small_primes:
        if a+6 in small_primes:
    print leng,dig,time.clock(),checked
    for loc in random.sample(range(leng),50):
        if is_prime(int('1'*loc+str(dig)+'1'*(leng-loc-1))):
            print leng-1,loc,dig,time.clock(),time.clock()-start, \
  • \$\begingroup\$ I know nothing more than the link, unfortunately. I found the link here: codegolf.stackexchange.com/questions/10739/… First answer \$\endgroup\$
    – isaacg
    Jul 31, 2014 at 5:20
  • \$\begingroup\$ Well, then. I'll credit you. \$\endgroup\$
    – isaacg
    Jul 31, 2014 at 5:35
  • 10
    \$\begingroup\$ It's like an ASCII where's wally... \$\endgroup\$ Jul 31, 2014 at 12:57
  • 5
    \$\begingroup\$ Maybe you should rename the function is_very_very_very_very_very_very_very_probably_prime()... \$\endgroup\$ Jul 31, 2014 at 12:58
  • 2
    \$\begingroup\$ Mathmatica and Maple both use the same method, so it can't be all that bad. \$\endgroup\$
    – primo
    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.

        my(N=10^d\9, P=select(p->isprime(d+p),[1,2,4,6]), D, n=1);
        if(#P==0, next);
                if (ispseudoprime(D=N+n*P[1]), print(D));
                if (ispseudoprime(k), print(k))
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.

  • \$\begingroup\$ How did you know to start at 4100? \$\endgroup\$
    – isaacg
    Jul 31, 2014 at 20:10
  • \$\begingroup\$ @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. \$\endgroup\$
    – Charles
    Jul 31, 2014 at 20:15
  • \$\begingroup\$ 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. \$\endgroup\$
    – isaacg
    Jul 31, 2014 at 20:23
  • \$\begingroup\$ @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). \$\endgroup\$
    – Charles
    Jul 31, 2014 at 20:28
  • \$\begingroup\$ @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. \$\endgroup\$
    – Charles
    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";
  • \$\begingroup\$ +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 \$\endgroup\$
    – primo
    Aug 2, 2014 at 6:49
  • \$\begingroup\$ 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. \$\endgroup\$
    – DanaJ
    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.



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]


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.


{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]]


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


And the sum of the digits is also a prime.

digitSum = Plus @@ digits


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.

  • 3
    \$\begingroup\$ If I am not mistaken, its sum is 4009 which is not prime. \$\endgroup\$
    – gerw
    Jul 31, 2014 at 13:20
  • \$\begingroup\$ 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. \$\endgroup\$
    – Johnride
    Jul 31, 2014 at 13:20
  • 2
    \$\begingroup\$ @Johnride: Both should be prime. \$\endgroup\$
    – gerw
    Jul 31, 2014 at 13:25
  • \$\begingroup\$ @gerw Is right. The number of digits AND the sum of the digits AND the product of digits all have to be prime. \$\endgroup\$ Jul 31, 2014 at 16:04
  • \$\begingroup\$ 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. \$\endgroup\$
    – DavidC
    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:

        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:

    print('Tasked with: {0}, {1}'.format(num_length, ns))
    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:

            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():
        if prime > 20006:
    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))
  • \$\begingroup\$ it would read more cleanly if you represent the codepad link as a [broken, if necessary] import \$\endgroup\$
    – Sparr
    Aug 1, 2014 at 3:05
  • \$\begingroup\$ I think that would be confusing, as the code at the other end isn't really importable like that. \$\endgroup\$
    – mkoistinen
    Aug 1, 2014 at 3:08
  • \$\begingroup\$ use isaacg's syntax. comment the URL, then import from a nonexistent package (my_math, in his case) \$\endgroup\$
    – Sparr
    Aug 1, 2014 at 3:14
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – isaacg
    Aug 1, 2014 at 3:37
  • \$\begingroup\$ 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. \$\endgroup\$
    – primo
    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.


#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);
                        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;

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);

    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);

    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


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):
    add = 0
    mul = 1
    while num:
         add, mul, num = add + num % 10, mul * (num % 10), num / 10
    return add, mul

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]:
            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)):
        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:
    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

Mathematica 4211 4259 digits

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

Which was generated by the following code:

  p = Prime[n];
  curlargest = Catch[
    If[PrimeQ[p + 6],
     list = ConstantArray[1, p];
      temp = FromDigits[ReplacePart[list, i -> 7]];
       ], {i, p}]
    If[PrimeQ[p + 4],
     list = ConstantArray[1, p];
      temp = FromDigits[ReplacePart[list, i -> 5]];
       ], {i, p}]
    If[PrimeQ[p + 2],
     list = ConstantArray[1, p];
      temp = FromDigits[ReplacePart[list, i -> 3]];
       ], {i, p}]
  , {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.)

enter image description here enter image description here
enter image description here

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?)

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

Python 3.9.2, with sympy and gmpy2, 3019 digits

Completed the search from 3000 digits to 3019 digits in approximately 10 minutes on my laptop.

Port of @Charles's Pari/GP answer in Python.

import sympy
import gmpy2
import time

def supreme(lim, startAt=3):
    for d in sympy.primerange(startAt, lim + 1):
        N = 10**d // 9
        P = [p for p in [1,2,4,6] if sympy.isprime(d + p)]
        if len(P) == 0:
        if len(P) == 1:
            n = 1
            for i in range(0, d):
                if gmpy2.is_strong_bpsw_prp(N + n * P[0]):
                    print(N + n * P[0])
                n *= 10
        D = [P[i+1] - P[i] for i in range(len(P) - 1)]
        n = 1
        for i in range(0, d):
            for k in range(N + n * P[0], N + n * P[-1] + 1, n * D[0]):
                if gmpy2.is_strong_bpsw_prp(k):
            n *= 10

supreme(3019, 3000)
print(f"The search ends in {time.time()-start_time} seconds")
The search ends in 606.180704832077 seconds

Axiom, 3019 digits {318}5{2700}

)time on

-- Return true if n is pseudo prime else return false
-- **Can Fail**
       if i>sq or i>50000 then break
       if n rem i=0       then return false
     if i~=50001        then return true
     if powmod(3,n,n)=3 then return true

-- 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==
    for i in 1..n-1 repeat
       --output b
       k:=0    -- 3 5 7 <-> 2 4 6
       for i in [2,4,6] repeat
           ~prime?(n+i)=>0 --somma
           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
       b:=b quo 10
       if b<=0 then break

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

result from the start value 3000 in 529 sec

(4) -> find(3000)

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.