21
\$\begingroup\$

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#)
\$\endgroup\$
18
  • 5
    \$\begingroup\$ Even 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. \$\endgroup\$
    – Rainbolt
    Nov 19, 2014 at 16:40
  • 2
    \$\begingroup\$ @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 codegolf.stackexchange.com/questions/41021/…) \$\endgroup\$ Nov 19, 2014 at 16:49
  • 1
    \$\begingroup\$ Related \$\endgroup\$ Nov 19, 2014 at 16:49
  • 7
    \$\begingroup\$ You're still leaving at a disadvantage those who have slower computers \$\endgroup\$ Nov 19, 2014 at 16:57
  • 11
    \$\begingroup\$ It 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 \$\endgroup\$
    – ruakh
    Nov 20, 2014 at 6:11

6 Answers 6

13
\$\begingroup\$

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.

\$\endgroup\$
5
  • 3
    \$\begingroup\$ "in about 1 minutes" - sorry, have to report a pluralization "bug". :P \$\endgroup\$
    – hichris123
    Nov 19, 2014 at 23:03
  • \$\begingroup\$ 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. \$\endgroup\$ Nov 20, 2014 at 14:00
  • \$\begingroup\$ 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? \$\endgroup\$
    – isaacg
    Nov 20, 2014 at 20:15
  • 1
    \$\begingroup\$ Because 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. \$\endgroup\$
    – Jakube
    Nov 20, 2014 at 20:40
  • \$\begingroup\$ Ah, very clever. \$\endgroup\$
    – isaacg
    Nov 20, 2014 at 20:41
9
\$\begingroup\$

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.
\$\endgroup\$
3
  • 1
    \$\begingroup\$ Thanks for the credit, but your algorithm is vastly more complex than mine! Great work, and welcome to PPCG! :) \$\endgroup\$ Nov 21, 2014 at 18:07
  • \$\begingroup\$ @FryAmTheEggman: thank you for the idea! It's inspiring. \$\endgroup\$
    – Qualtagh
    Nov 24, 2014 at 5:21
  • \$\begingroup\$ 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. \$\endgroup\$ Nov 24, 2014 at 8:05
7
\$\begingroup\$

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

\$\endgroup\$
5
  • \$\begingroup\$ 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. \$\endgroup\$ Nov 19, 2014 at 19:50
  • \$\begingroup\$ @SuboptimusPrime Oh, thanks. I didn't even look at it :P \$\endgroup\$ Nov 19, 2014 at 19:50
  • \$\begingroup\$ @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 :) \$\endgroup\$ Nov 19, 2014 at 20:30
  • \$\begingroup\$ Small correction to your answer: 91 = 13x7, so 91 is composite, and fragile primes ending in 1 can actually start with 9. \$\endgroup\$ Nov 19, 2014 at 21:59
  • \$\begingroup\$ @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. \$\endgroup\$ Nov 19, 2014 at 22:36
7
\$\begingroup\$

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

\$\endgroup\$
5
  • 2
    \$\begingroup\$ Your 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. \$\endgroup\$
    – Jakube
    Nov 20, 2014 at 0:06
  • 1
    \$\begingroup\$ Ah, 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. \$\endgroup\$
    – Emil
    Nov 20, 2014 at 6:55
  • 1
    \$\begingroup\$ Checked it again, now your results are correct. \$\endgroup\$
    – Jakube
    Nov 20, 2014 at 7:47
  • \$\begingroup\$ Your 5048 digit number is, indeed, a fragile prime according to my program. \$\endgroup\$ Nov 22, 2014 at 8:12
  • \$\begingroup\$ @SuboptimusPrime: Great, thanks for checking! \$\endgroup\$
    – Emil
    Nov 22, 2014 at 8:29
4
\$\begingroup\$

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);
    }
}
\$\endgroup\$
6
  • \$\begingroup\$ 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. \$\endgroup\$
    – Qualtagh
    Nov 26, 2014 at 11:46
  • \$\begingroup\$ @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. \$\endgroup\$ Nov 26, 2014 at 16:35
  • \$\begingroup\$ 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. \$\endgroup\$
    – Qualtagh
    Nov 28, 2014 at 4:26
  • 1
    \$\begingroup\$ @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). \$\endgroup\$ Nov 28, 2014 at 6:57
  • 2
    \$\begingroup\$ MPIR 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. \$\endgroup\$ Nov 28, 2014 at 20:12
2
\$\begingroup\$

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)
\$\endgroup\$
16
  • \$\begingroup\$ I think you can remove everything except the last 3... \$\endgroup\$
    – Sp3000
    Nov 20, 2014 at 7:42
  • \$\begingroup\$ it ends in 5 if i remove the last 3 so is divisible by 5 \$\endgroup\$ Nov 20, 2014 at 7:44
  • 2
    \$\begingroup\$ No I mean removing everything until you have just the last 3 left, which is prime. \$\endgroup\$
    – Sp3000
    Nov 20, 2014 at 7:45
  • 1
    \$\begingroup\$ @JohnMeacham My program suggests that this number turns into prime if you remove 386 digits from the left. \$\endgroup\$ Nov 20, 2014 at 12:20
  • 1
    \$\begingroup\$ Please, verify your numbers before posting. If you remove the left 1256 digits from your 1276-digit number, you'll get 99999994999999999999, which is prime. \$\endgroup\$ Nov 20, 2014 at 14:23

Your Answer

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

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