| Bytes | Lang | Time | Link |
|---|---|---|---|
| 3019 | Python 3.9.2 | 230722T064055Z | 138 Aspe |
| nan | Perl | 140801T210432Z | DanaJ |
| nan | Axiom | 170920T153448Z | user5898 |
| nan | JavaScript | 140827T003237Z | JeffSB |
| nan | C | 140801T072335Z | Vectoriz |
| nan | Python 2.7 | 140801T025205Z | mkoistin |
| nan | PFGW | 140801T122213Z | Tim S. |
| nan | 140801T024207Z | Tally | |
| nan | Python 2.7 on PyPy | 140731T045409Z | isaacg |
| 4127 | PARI/GP | 140731T184736Z | Charles |
| nan | 140731T121116Z | DavidC | |
| 1719 | Python 2.7 | 140731T054118Z | Sparr |
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:
continue
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
continue
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):
print(k)
n *= 10
start_time=time.time()
supreme(3019, 3000)
print(f"The search ends in {time.time()-start_time} seconds")
1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111151111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111115111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
The search ends in 606.180704832077 seconds
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";
}
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
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?)
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;
}
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:
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.
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))
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.
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
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
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.
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.
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