Find the largest fragile prime

21

3

Consider the function Remove(n, startIndex, count) that removes count digits from the number n starting from the digit at position startIndex. Examples:

Remove(1234, 1, 1) = 234
Remove(123456, 2, 3) = 156
Remove(1507, 1, 2) = 07 = 7
Remove(1234, 1, 4) = 0

We will call the prime number X fragile if every possible Remove operation makes it non-prime. For example, 80651 is a fragile prime because all of the following numbers are not prime:

651, 51, 1, 0, 8651, 851, 81, 8, 8051, 801, 80, 8061, 806, 8065

Goal

Write a program that finds the largest fragile prime. Edit: removed the time limit because there was a relatively fair way to circumvent it.

The score is the fragile prime number found by your program. In the event of a tie, the earlier submission wins.

Rules

  • You may use any language and any third-party libraries.
  • You run the program on you own hardware.
  • You may use probabilistic primality tests.
  • Everything is in base 10.

Leading entries

  • 6629 digits by Qualtagh (Java)
  • 5048 digits by Emil (Python 2)
  • 2268 digits by Jakube (Python 2)

Edit: I've added my own answer.

  • 28164 digits by Suboptimus Prime, based on Qualtagh's algorithm (C#)

Suboptimus Prime

Posted 2014-11-19T16:28:11.047

Reputation: 403

5Even if I don't hard code the answer, I could start searching at a point very close to a large fragile prime. Obviously, nobody wants to start searching at 1. What's stopping me from doing this? Exactly how close can I start my search before I get called out for basically hard-coding the answer? I love the challenge by the way. – Rainbolt – 2014-11-19T16:40:06.240

@Rainbolt You may start searching from any number you want, but you are not allowed to use numbers found by running your program for more than 15 minutes. Obviously, I can't enforce this rule, but hopefully everyone will play fair. – Suboptimus Prime – 2014-11-19T16:47:27.970

2

@SuboptimusPrime You could instead remove the time limit altogether, because I believe at some point that will be so rare that it'll be a feat to find the next one anyway. (Similar to http://codegolf.stackexchange.com/questions/41021/find-highest-scoring-matrix-without-property-x)

– Martin Ender – 2014-11-19T16:49:01.323

1Related – FryAmTheEggman – 2014-11-19T16:49:56.663

1@SuboptimusPrime Fair enough. I'm kind of leaning towards hoping everyone plays fair now, because as some folks pointed out in the chat, its been done at least twice before, and it worked. Leave it to the downvoters to beat up anyone who tries to cheat. – Rainbolt – 2014-11-19T16:53:14.877

@MartinBüttner Not everyone can keep their PC's running for hours or days, so they will be at a disadvantage if I remove the time limit. I want to make this challenge more accessible for everyone. – Suboptimus Prime – 2014-11-19T16:56:13.080

7You're still leaving at a disadvantage those who have slower computers – John Dvorak – 2014-11-19T16:57:14.930

@JanDvorak True, but I don't think I can do anything about it. That's a universal problem for all optimization challenges I guess. – Suboptimus Prime – 2014-11-19T17:00:49.200

1@SuboptimusPrime Well, you could test all submissions yourself (that's what people do for fastest-code challenges). – Martin Ender – 2014-11-19T17:18:21.457

@MartinBüttner I've based my rules on this question. If I'll test all submissions myself than other people would have to optimize their programs for my hardware, which doesn't look like a good idea to me. And I'm still unsure about that time limit rule. I'll remove it if there will be enough demand for it.

– Suboptimus Prime – 2014-11-19T17:31:19.473

1I finally realized that the time limit was pointless. You could run the program repeatedly on the consecutive intervals of numbers, 15 minutes at a time. That would effectively remove the time limit without cheating, because the rules allow you to start searching from any number you want. – Suboptimus Prime – 2014-11-19T21:48:37.630

11It took me an embarrassingly long time to realize that "Write a program that finds the largest fragile prime" did not mean "There exists a largest fragile prime. Write a program that finds it." I guess I've done too much Project Euler. :-P – ruakh – 2014-11-20T06:11:08.630

"You may use probabilistic primality tests." So if I have a number p, and 2**(p-1) % p = 1, I can claim p is prime? – feersum – 2014-11-20T22:08:35.777

@feersum You may use any algorithm you want but it has to produce the correct result, so it's probably better to use something more reliable than 2**(p-1) % p = 1. – Suboptimus Prime – 2014-11-21T06:43:08.333

I don't get it. If the result has to be correct, how is it then OK if it may be incorrect? – feersum – 2014-11-21T06:45:45.280

@feersum I'm veryfing all the answers with a probabilistic test that has a very low possibility of error. You answer has to pass this test to be accepted. – Suboptimus Prime – 2014-11-21T06:49:43.263

1So the result need not be correct, if someone can guess which algorithm you are using ;) – feersum – 2014-11-21T06:52:11.333

@feersum I'm using 40 random miller-rabin tests as a primality test. I don't think it's possible to adjust any algorithm to pass 40 random tests :) – Suboptimus Prime – 2014-11-21T07:00:49.667

Answers

9

Java - 3144 3322 6629 digits

6 0{3314} 8969999

6 0{6623} 49099

This solution is based on FryAmTheEggman's answer.

  1. The last digit is 1 or 9.
  2. If the last digit is 1 then a previous one is 0, 8 or 9.
  3. If the last digit is 9 then a previous one is 0, 4, 6 or 9.
  4. ...

What if we dig deeper?

It becomes a tree structure:

                        S
             -----------------------
             1                     9
    ------------------         ----------------
    0           8    9         0    4    6    9
---------     -----
0   8   9      ...

Let's call number R right composite if R and all its endings are composite.

We're gonna iterate over all right composite numbers in breadth-first manner: 1, 9, 01, 81, 91, 09, 49, 69, 99, 001, 801, 901 etc.

Numbers starting with zero are not checked for primality but are needed to build further numbers.

We will look for a target number N in the form X00...00R, where X is one of 4, 6, 8 or 9 and R is right composite. X cannot be prime. X cannot be 0. And X cannot be 1 because if R ends with 1 or 9 then N would contain 11 or 19.

If XR contains prime numbers after "remove" operation then XYR would contain them too for any Y. So we shouldn't traverse branches starting from R.

Let X be a constant, say 6.

Pseudocode:

X = 6;
for ( String R : breadth-first-traverse-of-all-right-composites ) {
  if ( R ends with 1 or 9 ) {
    if ( remove( X + R, i, j ) is composite for all i and j ) {
      for ( String zeros = ""; zeros.length() < LIMIT; zeros += "0" ) {
        if ( X + zeros + R is prime ) {
          // At this step these conditions hold:
          // 1. X + 0...0 is composite.
          // 2. 0...0 + R = R is composite.
          // 3. X + 0...0 + R is composite if 0...0 is shorter than zeros.
          suits = true;
          for ( E : all R endings )
            if ( X + zeros + E is prime )
              suits = false;
          if ( suits )
            print R + " is fragile prime";
          break; // try another R
                 // because ( X + zeros + 0...0 + R )
                 // would contain prime ( X + zeros + R ).
        }
      }
    }
  }
}

We should limit zeros quantity because it may take too long to find a prime number in form X + zeros + R (or forever if all of them are composite).

The real code is quite verbose and can be found here.

Primality testing for numbers in long int range is performed by deterministic variant of Miller test. For BigInteger numbers a trial division is performed first and then BailliePSW test. It's probabilistic but quite certain. And it's faster than Miller-Rabin test (we should do many iterations for such big numbers in Miller-Rabin to gain enough accuracy).

Edit: the first attempt was incorrect. We should also ignore branches starting with R if X0...0R is prime. Then X0...0YR wouldn't be fragile prime. So an additional check was added. This solution seems to be correct.

Edit 2: added an optimization. If ( X + R ) is divisible by 3 then ( X + zeros + R ) is also divisible by 3. So ( X + zeros + R ) cannot be prime in this case and such R's may be skipped.

Edit 3: it was not necessary to skip prime digits if they're not in the last or first position. So endings like 21 or 51 are ok. But it doesn't change anything much.

Conclusions:

  1. My last answer was checking for being fragile for 100 minutes. The search of the answer (checking all preceding variants) took about 15 minutes. Yes, it makes no sense to restrict search time (we can start searching from the target number, so time would be zero). But it could be meaningful to restrict checking time like in this question.
  2. The answer 60...049099 has digit 4 in the middle. If the "remove" operation touches it, the number becomes divisible by 3. So we should check remove operations in left and right sides. Right side is too short. Left side length is almost n = length( N ).
  3. Primality tests like BPSW and Miller-Rabin use constant quantity of modular exponentiations. Its complexity is O( M( n ) * n ) according to this page, where M( n ) is multiplication complexity. Java uses Toom-Cook and Karatsuba algorithms but we'll take scholar algorithm for simplicity. M( n ) = n2. So primality testing complexity is O( n3 ).
  4. We should check all numbers from length = 6 to 6629. Let's take min = 1 and max = n for commonality. The whole check complexity is O( 13 + 23 + ... + n3 ) = O( ( n * ( n + 1 ) / 2 )2 ) = O( n4 ).
  5. Emil's answer has the same checking asymptotics. But the constant factor is lower. Digit "7" is standing in the middle of the number. Left side and right side can be almost equal. It gives ( n / 2 )4 * 2 = n4 / 8. Speedup: 8X. Numbers in the form 9...9Y9...9 can be 1.7 times longer than in the form X0...0R having the same checking time.

Qualtagh

Posted 2014-11-19T16:28:11.047

Reputation: 206

1Thanks for the credit, but your algorithm is vastly more complex than mine! Great work, and welcome to PPCG! :) – FryAmTheEggman – 2014-11-21T18:07:37.687

@FryAmTheEggman: thank you for the idea! It's inspiring. – Qualtagh – 2014-11-24T05:21:32.213

Your analysis of check complexity is very interesting, but search compexity is probably important too. I think that your algorithm requires significantly less primality tests of large numbers (compared to Emil's) to find a large fragile prime. And you can speed up primality tests by using a native library. I'm using Mpir.NET, and checking your number for being a fragile prime takes just a few minutes. – Suboptimus Prime – 2014-11-24T08:05:29.597

13

Python 2 - 126 1221 1337 1719 2268 digits

999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999799999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999

'9' * 1944 + '7' + '9' * 323

There are at about len(n)^2 resulting numbers of Remove(n, startIndex, count). I tried to minimize those numbers. If there are lots of digits next to each other are the same, lots of these resulting numbers can be ignored, because they appear multiple times.

So I took it to the extreme, only 9s and a little prime in the middle. I also took a look at the fragile prime under 1 million, and saw, that there are such fragile prime. Searching for numbers with 2 9s at the end works really good, not sure why. 1 number, 3, or 4 9s at the end results in smaller fragile primes.

It uses the pyprimes module. I'm not sure, if it is any good. It uses the miller_rabin test, so it's probabilistic.

The program find this 126-digit fragile prime in about 1 minute, and for the rest of the time it searches without success.

biggest_found = 80651

n = lambda a,b,c: '9'*a + b + '9'*c

for j in range(1000):
   for digit in '124578':
      for i in range(2000):
         number = int(n(i,digit,j))
         if is_prime(number):
            if (number > biggest_found):
               if all(not is_prime(int(n(i,digit,k))) for k in range(j)):
                  biggest_found = number
                  print(i+j+1, biggest_found)
            break

edit:

Just saw, that you removed the time limit. I will run the program over night, maybe some really big fragile primes appear.

edit 2:

Made my original program faster, so but still no solution with more than 126 digits. So I jumped on the train and searched for x 9s + 1 digit + y 9s. The advantage is, that you have to check O(n) numbers for primality, if you fixe y. It finds a 1221 rather quickly.

edit 3:

For the 2268 digit number I use the same program, only divided the work on multiple cores.

Jakube

Posted 2014-11-19T16:28:11.047

Reputation: 21 462

3"in about 1 minutes" - sorry, have to report a pluralization "bug". :P – hichris123 – 2014-11-19T23:03:56.993

The probabilistic nature of the miller-rabin is what was biting me for my last few entries. You might want to verify with another algorithm as well. – John Meacham – 2014-11-20T14:00:04.010

Why do you only check that the numbers formed from removing digits from the end are composite? Why not check the numbers formed by removing digits from the front? – isaacg – 2014-11-20T20:15:13.247

1Because I checked these before in the 'for i'-loop. Here I append 9s at the beginning, and make a prime check. When I find the first prime number of this form, I know, that all numbers with less 9s at the beginning are not prime. And after checking for removing 9s at the end, I stop (break), because now, every number has a prime number in it and is therefore not prime. – Jakube – 2014-11-20T20:40:29.937

Ah, very clever. – isaacg – 2014-11-20T20:41:51.763

7

Python 2.7 - 429623069 99993799

No optimizations whatsoever, so far. Just using some trivial observations about fragile primes (thanks to Rainbolt in chat):

  1. Fragile primes must end in 1 or 9 (Primes are not even, and the final digit must not be prime)
  2. Fragile primes ending in 1 must start with 8 or 9 (the first number can't be prime, and 11, 41 and 61 and are all primes)
  3. Fragile primes ending in 9 must start with 4,6 or 9 (see reasoning for 1, but only 89 is prime)

Just trying to get the ball rolling :)

This technically runs slightly over 15 minutes, but it only checks a single number in the extra time.

is_prime is taken from here (isaacg used it here) and is probabilistic.

def substrings(a):
    l=len(a)
    out=set()
    for i in range(l):
        for j in range(l-i):
            out.add(a[:i]+a[len(a)-j:])
    return out

import time

n=9
while time.clock()<15*60:
    if is_prime(n):
        if not any(map(lambda n: n!='' and is_prime(int(n)), substrings(`n`))):
            print n
    t=`n`
    if n%10==9 and t[0]=='8':n+=2
    elif n%10==1 and t[0]!='8':n+=8
    elif t[0]=='1' or is_prime(int(t[0])):n+=10**~-len(t)
    else:n+=10

Just a note, when I start this with n=429623069 I get up to 482704669. The extra digit really seems to kill this strategy...

FryAmTheEggman

Posted 2014-11-19T16:28:11.047

Reputation: 16 206

Not bad for a start! Although it seems that is_prime performs a full deteterministic check for 32-bit values, which is a little bit excessive. I think is_prime method may work faster if you'll comment out the full trial division part. – Suboptimus Prime – 2014-11-19T19:50:09.287

@SuboptimusPrime Oh, thanks. I didn't even look at it :P – FryAmTheEggman – 2014-11-19T19:50:43.910

@SuboptimusPrime I think the full deterministic check is faster for small values because the author defined steps to take in between candidate factors. Thanks again for the idea, but it seems much faster when leaving that in :) – FryAmTheEggman – 2014-11-19T20:30:26.583

Small correction to your answer: 91 = 13x7, so 91 is composite, and fragile primes ending in 1 can actually start with 9. – Suboptimus Prime – 2014-11-19T21:59:40.743

@SuboptimusPrime Quite right, dunno how I messed that up. The value I posted should still be valid, as I was just skipping some possible values. – FryAmTheEggman – 2014-11-19T22:36:05.777

7

Python 2, 828 digits 5048 digits

99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999799999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999
155*'9'+'7'+4892*'9'

As @Jakube pointed out, the first prime I submitted was not actually fragile due to a bug in my code. Fixing the bug was easy but it also made the algorithm significantly slower.

I limited myself to an easily searchable subset of the fragile primes, namely those that only consist of the digit 9 and exactly one digit 7.

def fragile_prime_generator(x, b_max):
  bs, cs = set(), set()
  prime = dict()

  def test_prime(b,c):
    if (b,c) not in prime:
      prime[(b,c)] = is_prime(int('9'*b+`x`+'9'*c))
    return prime[(b,c)]

  def test_frag(b,c):
    for b2 in xrange(b):
      if test_prime(b2,c):
        bs.add(b2)
        return False
    for c2 in xrange(c):
      if test_prime(b,c2):
        cs.add(c2)
        return False
    return True

  a = 1
  while len(bs)<b_max:
    for b in xrange(min(a, b_max)):
      c = a-b
      if b not in bs and c not in cs and test_prime(b,c):
        bs.add(b)
        cs.add(c)
        if test_frag(b,c): yield b,c
    a += 1
  print "no more fragile primes of this form"

for b,c in fragile_prime_generator(7, 222):
  print ("%d digit fragile prime found: %d*'9'+'%d'+%d*'9'"
          % (b+c+1, b, x, c))

I used the same is_prime function (from here) as @FryAmTheEggman.

Edit:

I made two changes to make the algorithm faster:

  • I try to skip as many primality checks as possible and only go back when a potential fragile prime is found to make sure it's really fragile. There's a small number of duplicate checks, so I crudely memoized the prime checking function.

  • For the numbers of the form b*'9' + '7' + c*'9' I limited the size of b. The lower the limit, the less numbers have to be checked, but chances increase to not find any large fragile prime at all. I kind of arbitrarily chose 222 as the limit.

At a few thousand digits a single prime check already can take my program a few seconds. So, I probably can't do much better with this approach.

Please feel free to check the correctness of my submission. Due to the probabilistic primality check my number could theoretically not be prime, but if it is, it should be fragile. Or I've done something wrong. :-)

Emil

Posted 2014-11-19T16:28:11.047

Reputation: 1 438

2Your found prime is not fragile. If you call Remove(n,83,838) [Remove everything except the first 82 digits], you'll end up with a prime. – Jakube – 2014-11-20T00:06:39.630

1Ah, thanks @Jakube. I was trying to be too clever. Turns out I was skipping more primality checks then I should have. I'm on my way to fix it. – Emil – 2014-11-20T06:55:48.900

1Checked it again, now your results are correct. – Jakube – 2014-11-20T07:47:54.647

Your 5048 digit number is, indeed, a fragile prime according to my program. – Suboptimus Prime – 2014-11-22T08:12:01.800

@SuboptimusPrime: Great, thanks for checking! – Emil – 2014-11-22T08:29:58.820

4

C#, 10039 28164 digits

6 0{28157} 169669

Edit: I've made another program based on Qualtagh's algorithm with some minor modifications:

  • I'm searching for the numbers of the form L000...000R, where L is left composite, R is right composite. I allowed the left composite number L to have multiple digits, although this is mostly a stylistic change, and it probably doesn't affect the efficiency of the algorithm.
  • I've added multithreading to speed up the search.
using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Threading;
using System.Threading.Tasks;
using Mpir.NET;

class Program
{
    const int PrimeNotFound = int.MaxValue;

    private static BitArray _primeSieve;
    private static HashSet<Tuple<int, int>> _templatesToSkip = new HashSet<Tuple<int, int>>();

    static void Main(string[] args)
    {
        int bestDigitCount = 0;
        foreach (Tuple<int, int> template in GetTemplates())
        {
            int left = template.Item1;
            int right = template.Item2;
            if (SkipTemplate(left, right))
                continue;

            int zeroCount = GetZeroCountOfPrime(left, right);
            if (zeroCount != PrimeNotFound)
            {
                int digitCount = left.ToString().Length + right.ToString().Length + zeroCount;
                if (digitCount >= bestDigitCount)
                {
                    string primeStr = left + " 0{" + zeroCount + "} " + right;
                    Console.WriteLine("testing " + primeStr);
                    bool isFragile = IsFragile(left, right, zeroCount);
                    Console.WriteLine(primeStr + " is fragile: " + isFragile);

                    if (isFragile)
                        bestDigitCount = digitCount;
                }

                _templatesToSkip.Add(template);
            }
        }
    }

    private static int GetZeroCountOfPrime(int left, int right)
    {
        _zeroCount = 0;

        int threadCount = Environment.ProcessorCount;
        Task<int>[] tasks = new Task<int>[threadCount];
        for (int i = 0; i < threadCount; i++)
            tasks[i] = Task.Run(() => InternalGetZeroCountOfPrime(left, right));
        Task.WaitAll(tasks);

        return tasks.Min(task => task.Result);
    }

    private static int _zeroCount;

    private static int InternalGetZeroCountOfPrime(int left, int right)
    {
        const int maxZeroCount = 40000;
        int zeroCount = Interlocked.Increment(ref _zeroCount);
        while (zeroCount <= maxZeroCount)
        {
            if (zeroCount % 1000 == 0)
                Console.WriteLine("testing " + left + " 0{" + zeroCount + "} " + right);

            if (IsPrime(left, right, zeroCount))
            {
                Interlocked.Add(ref _zeroCount, maxZeroCount);
                return zeroCount;
            }
            else
                zeroCount = Interlocked.Increment(ref _zeroCount);
        }

        return PrimeNotFound;
    }

    private static bool SkipTemplate(int left, int right)
    {
        for (int leftDiv = 1; leftDiv <= left; leftDiv *= 10)
            for (int rightDiv = 1; rightDiv <= right; rightDiv *= 10)
                if (_templatesToSkip.Contains(Tuple.Create(left / leftDiv, right % (rightDiv * 10))))
                    return true;

        return false;
    }

    private static bool IsPrime(int left, int right, int zeroCount)
    {
        return IsPrime(left.ToString() + new string('0', zeroCount) + right.ToString());
    }

    private static bool IsPrime(string left, string right, int zeroCount)
    {
        return IsPrime(left + new string('0', zeroCount) + right);
    }

    private static bool IsPrime(string s)
    {
        using (mpz_t n = new mpz_t(s))
        {
            return n.IsProbablyPrimeRabinMiller(20);
        }
    }

    private static bool IsFragile(int left, int right, int zeroCount)
    {
        string leftStr = left.ToString();
        string rightStr = right.ToString();

        for (int startIndex = 0; startIndex < leftStr.Length - 1; startIndex++)
            for (int count = 1; count < leftStr.Length - startIndex; count++)
                if (IsPrime(leftStr.Remove(startIndex, count), rightStr, zeroCount))
                    return false;

        for (int startIndex = 1; startIndex < rightStr.Length; startIndex++)
            for (int count = 1; count <= rightStr.Length - startIndex; count++)
                if (IsPrime(leftStr, rightStr.Remove(startIndex, count), zeroCount))
                    return false;

        return true;
    }

    private static IEnumerable<Tuple<int, int>> GetTemplates()
    {
        const int maxDigitCount = 8;
        PreparePrimeSieve((int)BigInteger.Pow(10, maxDigitCount));
        for (int digitCount = 2; digitCount <= maxDigitCount; digitCount++)
        {
            for (int leftCount = 1; leftCount < digitCount; leftCount++)
            {
                int rightCount = digitCount - leftCount;
                int maxLeft = (int)BigInteger.Pow(10, leftCount);
                int maxRight = (int)BigInteger.Pow(10, rightCount);

                for (int left = maxLeft / 10; left < maxLeft; left++)
                    for (int right = maxRight / 10; right < maxRight; right++)
                        if (IsValidTemplate(left, right, leftCount, rightCount))
                            yield return Tuple.Create(left, right);
            }

        }
    }

    private static void PreparePrimeSieve(int limit)
    {
        _primeSieve = new BitArray(limit + 1, true);
        _primeSieve[0] = false;
        _primeSieve[1] = false;

        for (int i = 2; i * i <= limit; i++)
            if (_primeSieve[i])
                for (int j = i * i; j <= limit; j += i)
                    _primeSieve[j] = false;
    }

    private static bool IsValidTemplate(int left, int right, int leftCount, int rightCount)
    {
        int rightDigit = right % 10;
        if ((rightDigit != 1) && (rightDigit != 9))
            return false;

        if (left % 10 == 0)
            return false;

        if ((left + right) % 3 == 0)
            return false;

        if (!Coprime(left, right))
            return false;

        int leftDiv = 1;
        for (int i = 0; i <= leftCount; i++)
        {
            int rightDiv = 1;
            for (int j = 0; j <= rightCount; j++)
            {
                int combination = left / leftDiv * rightDiv + right % rightDiv;
                if (_primeSieve[combination])
                    return false;

                rightDiv *= 10;
            }

            leftDiv *= 10;
        }

        return true;
    }

    private static bool Coprime(int a, int b)
    {
        while (b != 0)
        {
            int t = b;
            b = a % b;
            a = t;
        }
        return a == 1;
    }
}

Old answer:

8 0{5436} 4 0{4600} 1

The are a few notable patterns for fragile primes:

600..00X00..009
900..00X00..009
800..00X00..001
999..99X99..999

where X can be 1, 2, 4, 5, 7 or 8.

For such numbers we only have to consider (length - 1) possible Remove operations. The other Remove operations produce either duplicates or obviously composite numbers. I tried to search for all such numbers with up to 800 digits and noticed that 4 patterns come up more frequently than the rest: 8007001, 8004001, 9997999 and 6004009. Since Emil and Jakube are using the 999X999 pattern, I decided to use 8004001 just to add some variety.

I've added the following optimizations to the algorithm:

  • I start searching from numbers with 7000 digits and then increment length by 1500 every time a fragile prime is found. If there is no fragile prime with a given length then I increment it by 1. 7000 and 1500 are just arbitrary numbers that seemed appropriate.
  • I'm using multithreading to search for the numbers with different length at the same time.
  • The result of each prime check is stored in a hash table to prevent duplicate checks.
  • I'm using Miller-Rabin implementation from Mpir.NET, which is very fast (MPIR is a fork of GMP).
using System;
using System.Collections.Concurrent;
using System.Threading.Tasks;
using Mpir.NET;

class Program
{
    const string _template = "8041";

    private static ConcurrentDictionary<Tuple<int, int>, byte> _compositeNumbers = new ConcurrentDictionary<Tuple<int, int>, byte>();
    private static ConcurrentDictionary<int, int> _leftPrimes = new ConcurrentDictionary<int, int>();
    private static ConcurrentDictionary<int, int> _rightPrimes = new ConcurrentDictionary<int, int>();

    static void Main(string[] args)
    {
        int threadCount = Environment.ProcessorCount;
        Task[] tasks = new Task[threadCount];
        for (int i = 0; i < threadCount; i++)
        {
            int index = i;
            tasks[index] = Task.Run(() => SearchFragilePrimes());
        }
        Task.WaitAll(tasks);
    }

    private const int _lengthIncrement = 1500;
    private static int _length = 7000;
    private static object _lengthLock = new object();
    private static object _consoleLock = new object();

    private static void SearchFragilePrimes()
    {
        int length;
        lock (_lengthLock)
        {
            _length++;
            length = _length;
        }

        while (true)
        {
            lock (_consoleLock)
            {
                Console.WriteLine("{0:T}: length = {1}", DateTime.Now, length);
            }

            bool found = false;
            for (int rightCount = 1; rightCount <= length - 2; rightCount++)
            {
                int leftCount = length - rightCount - 1;
                if (IsFragilePrime(leftCount, rightCount))
                {
                    lock (_consoleLock)
                    {
                        Console.WriteLine("{0:T}: {1} {2}{{{3}}} {4} {2}{{{5}}} {6}",
                            DateTime.Now, _template[0], _template[1], leftCount - 1,
                            _template[2], rightCount - 1, _template[3]);
                    }
                    found = true;
                    break;
                }
            }

            lock (_lengthLock)
            {
                if (found && (_length < length + _lengthIncrement / 2))
                    _length += _lengthIncrement;
                else
                    _length++;
                length = _length;
            }
        }
    }

    private static bool IsFragilePrime(int leftCount, int rightCount)
    {
        int count;
        if (_leftPrimes.TryGetValue(leftCount, out count))
            if (count < rightCount)
                return false;

        if (_rightPrimes.TryGetValue(rightCount, out count))
            if (count < leftCount)
                return false;

        if (!IsPrime(leftCount, rightCount))
            return false;

        for (int i = 0; i < leftCount; i++)
            if (IsPrime(i, rightCount))
                return false;

        for (int i = 0; i < rightCount; i++)
            if (IsPrime(leftCount, i))
                return false;

        return true;
    }

    private static bool IsPrime(int leftCount, int rightCount)
    {
        Tuple<int, int> tuple = Tuple.Create(leftCount, rightCount);
        if (_compositeNumbers.ContainsKey(tuple))
            return false;

        using (mpz_t n = new mpz_t(BuildStr(leftCount, rightCount)))
        {
            bool result = n.IsProbablyPrimeRabinMiller(20);

            if (result)
            {
                _leftPrimes.TryAdd(leftCount, rightCount);
                _rightPrimes.TryAdd(rightCount, leftCount);
            }
            else
                _compositeNumbers.TryAdd(tuple, 0);

            return result;
        }
    }

    private static string BuildStr(int leftCount, int rightCount)
    {
        char[] chars = new char[leftCount + rightCount + 1];
        for (int i = 0; i < chars.Length; i++)
            chars[i] = _template[1];
        chars[0] = _template[0];
        chars[leftCount + rightCount] = _template[3];
        chars[leftCount] = _template[2];
        return new string(chars);
    }
}

Suboptimus Prime

Posted 2014-11-19T16:28:11.047

Reputation: 403

While I was trying to verify your first answer, you have already posted a new one)). Checking already took 24 hours. The answer seems to be correct. I can't believe that Java's BigInteger is so much slower than native implementations. I thought about 2, 3 or even 10 times slower. But 24 hours against several minutes is too much. – Qualtagh – 2014-11-26T11:46:08.477

@Qualtagh To be fair, the 10039 digit number took 35 hours to find due to the inferior algorithm:) My current program takes about 3 minutes to find your 6629 digit number, and 6 hours to find the 28164 digit one. – Suboptimus Prime – 2014-11-26T16:35:59.680

Your first answer is correct. Verified! Verification took 48 hours. And I won't even try to verify the second answer)). I'm wondering why is BigInteger so slow comparing to MPIR. Is it only JVM/native difference? I set a "-server" flag, so expect the code to be JIT-compiled. Algorithms of modular exponentiation differ: both Java and MPIR use 2<sup>k</sup>-ary sliding window, but k=3 is fixed in Java and MPIR chooses k according to the size of the exponent. Is MPIR using parallel computations on several cores or probably GPU capabilities? Java's BigInteger doesn't. – Qualtagh – 2014-11-28T04:26:25.270

1

@Qualtagh I'm pretty sure MPIR is only using one CPU core. I've added multithreading myself which resulted in nearly 4 times faster search on a quad-core CPU. I didn't compare the internal implementation of MPIR and Java BigInteger, but I guess that MPIR is using better algorithms for multiplication and modular division. Also, it's probably better optimized for 64-bit CPUs (see the benchmark in this blog post).

– Suboptimus Prime – 2014-11-28T06:57:22.333

Wow, that's an interesting point. I didn't think about x64. One more idea: Java's BigInteger is immutable though it uses MutableBigInteger for some operations inside. Probably immutability leads to performance loss. But it shouldn't matter that much. Yes, MPIR has FFT and other fast multiplication algorithms that Java lacks of. Though I don't know if they are used in modular exponentiation. – Qualtagh – 2014-11-28T10:52:46.643

2MPIR is indeed single core and it doesn't use GPU. It's a highly optimized and fine tuned mix of C and assembler code. There's an MPIR version that uses only C (for portability reasons), but the C+ASM version is notably faster. The MPIR version I'm using for MPIR.Net is C+ASM using the K8 (1st gen x64) instruction set, because I wanted MPIR.Net to run on all x64 PCs. The versions for later instruction sets weren't noticeably faster in my crypto benchmark, but that may of course differ for other operations. – John Reynolds – 2014-11-28T20:12:44.350

2

Haskell - 1220 1277 digits fixed for really reals

99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999997999999999999999999999999999999999999999999999999999999999999999999999

9{1150} 7 9{69}

Better one - 1277 digits

9{871} 8 9{405}

Haskell code

downADigit :: Integer -> [Integer]
downADigit n = f [] 1 where
     f xs a | nma /= n = f (((n `div` a10)*a + nma):xs) a10
            | otherwise = xs where
        a10 = a * 10
        nma = n `mod` a

isFragile = all (not . isPrime') . downADigit
findNextPrime :: Integer -> Integer
findNextPrime n | even n = f (n + 1)
                | otherwise = f n where
    f n | isPrime' n  = n
        | otherwise = f (n + 2)

primesFrom n = f (findNextPrime n) where
    f n = n:f (findNextPrime $ n + 1)

primeLimit = 10000

isPrime' n | n < primeLimit = isPrime n
isPrime' n = all (millerRabinPrimality n) [2,3,5,7,11,13,17,19,984,7283,6628,8398,2983,9849,2739]

-- (eq. to) find2km (2^k * n) = (k,n)
find2km :: Integer -> (Integer,Integer)
find2km n = f 0 n
    where 
        f k m
            | r == 1 = (k,m)
            | otherwise = f (k+1) q
            where (q,r) = quotRem m 2        

-- n is the number to test; a is the (presumably randomly chosen) witness
millerRabinPrimality :: Integer -> Integer -> Bool
millerRabinPrimality n a
    | a <= 1 || a >= n-1 = 
        error $ "millerRabinPrimality: a out of range (" 
              ++ show a ++ " for "++ show n ++ ")" 
    | n < 2 = False
    | even n = False
    | b0 == 1 || b0 == n' = True
    | otherwise = iter (tail b)
    where
        n' = n-1
        (k,m) = find2km n'
        b0 = powMod n a m
        b = take (fromIntegral k) $ iterate (squareMod n) b0
        iter [] = False
        iter (x:xs)
            | x == 1 = False
            | x == n' = True
            | otherwise = iter xs

-- (eq. to) pow' (*) (^2) n k = n^k
pow' :: (Num a, Integral b) => (a->a->a) -> (a->a) -> a -> b -> a
pow' _ _ _ 0 = 1
pow' mul sq x' n' = f x' n' 1
    where 
        f x n y
            | n == 1 = x `mul` y
            | r == 0 = f x2 q y
            | otherwise = f x2 q (x `mul` y)
            where
                (q,r) = quotRem n 2
                x2 = sq x

mulMod :: Integral a => a -> a -> a -> a
mulMod a b c = (b * c) `mod` a
squareMod :: Integral a => a -> a -> a
squareMod a b = (b * b) `rem` a

-- (eq. to) powMod m n k = n^k `mod` m
powMod :: Integral a => a -> a -> a -> a
powMod m = pow' (mulMod m) (squareMod m)

-- simple for small primes
primes :: [Integer]
primes = 2:3:primes' where
    1:p:candidates = [6*k+r | k <- [0..], r <- [1,5]]
    primes'        = p : filter isPrime candidates
    isPrime n      = all (not . divides n)
                                   $ takeWhile (\p -> p*p <= n) primes'
    divides n p    = n `mod` p == 0
isPrime :: Integer -> Bool
isPrime n | n < 2 = False
          | otherwise = f primes where
            f (p:ps) | p*p <= n = if n `rem` p == 0 then False else f ps
                     | otherwise = True

main = do
    print . head $ filter isFragile (primesFrom $ 10^1000)

John Meacham

Posted 2014-11-19T16:28:11.047

Reputation: 129

I think you can remove everything except the last 3... – Sp3000 – 2014-11-20T07:42:43.293

it ends in 5 if i remove the last 3 so is divisible by 5 – John Meacham – 2014-11-20T07:44:54.260

2No I mean removing everything until you have just the last 3 left, which is prime. – Sp3000 – 2014-11-20T07:45:28.763

doh. you are right. just a sec. – John Meacham – 2014-11-20T07:47:14.593

Still not fragile. 61 (removing the first 1000 digits) is prime – Jakube – 2014-11-20T08:20:45.313

@Optimizer No, it don't match the lengths of 9s in the front or the back so a single remove could not make them match. – John Meacham – 2014-11-20T12:15:33.830

1@JohnMeacham My program suggests that this number turns into prime if you remove 386 digits from the left. – Suboptimus Prime – 2014-11-20T12:20:47.710

@JohnMeacham Oh, forgot that removal is in range .. – Optimizer – 2014-11-20T12:21:24.670

@SuboptimusPrime Hmm.. pari/gp seems to think it is composite. – John Meacham – 2014-11-20T12:28:07.457

@JohnMeacham My program also think's it's prime. And so does WolframAlpha: http://goo.gl/YmXW4q

– Jakube – 2014-11-20T12:33:39.340

strange, pari/gp is supposed to be quite good. I'll dig up some other CAS software and try to figure out what is going on. – John Meacham – 2014-11-20T12:52:43.100

Sure you are checking the same number? Because I just checked it again with pari/gp and it also says prime. We are talking about the number with 410 9s, than a 7 and than 253 9s = 10^664 - 1 - 2*10^253. – Jakube – 2014-11-20T12:59:22.460

Okay, I wrote a C program using openssl's crypto library to verify the primes. It looks like there were some spurious ones misidentified, I should have pumped up my probabalistic prime tester for how huge of numbers i am testing. A chance in a few hundred million of misidentifying a prime suddenly becomes relevant when working with a thousand digit numbers. :) I guess it makes me feel better that Maple had the same bug and ran into it on a measely 350 digit number. – John Meacham – 2014-11-20T13:39:23.857

Verified, the new number is a fragile prime. – Suboptimus Prime – 2014-11-20T13:53:15.700

1Please, verify your numbers before posting. If you remove the left 1256 digits from your 1276-digit number, you'll get 99999994999999999999, which is prime. – Suboptimus Prime – 2014-11-20T14:23:06.067

erf. cut-n-paste wrong line by one, 1277 was the verified one. – John Meacham – 2014-11-20T14:44:07.090