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

Question

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

Answer 1

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";
}

Answer 2

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

Answer 3

PARI/GP, 4127 digits

(10⁴¹²⁷-1)/9 + 2*10⁵¹⁵

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.

Answer 4

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.

Answer 5

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

Answer 6

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

Answer 7

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.

Answer 8

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

Answer 9

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

Answer 10

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

Answer 11

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