Belphegor Primes

October 31, 2014

Belphegor is one of the Seven Princes of Hell, charged with helping people make ingenious inventions and discoveries. Simon Singh gave the name Belphegor’s Prime to the number 1000000000000066600000000000001, which is the Number of the Beast, 666, from the Apocalypse, surrounded on each side by an unlucky 13 zeroes. Generally, Belphegor numbers Bn = (10^(n+3) + 666)*10^(n+1) + 1 have a 1, followed by n zeroes, followed by 666, followed by n zeroes, followed by 1. There are eight known Belphegor primes, with n ∈ {0, 13, 42, 506, 608, 2472, 2623, 28291} (A232448). You can read more about Belphegor primes at Cliff Pickover’s page.

Your task is to write a program that enumerates the Belphegor primes. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

5 Responses to “Belphegor Primes”

  1. Andras said

    Java, with probability chek:

    public static void main(String[] args) {
    int n = -1;
    while (true) {
    n++;
    if (new BigInteger(“1” + StringUtils.repeat(“0”, n) + “666” + StringUtils.repeat(“0”, n) + “1”).isProbablePrime(10)) {
    System.out.println(“n: ” + n + ” is probably (1 – 1/2^” + 10 + “) a prime.”);
    }
    }
    }

  2. Graham said

    In J, though I’m still a rookie:

    ip=:1=#@q: NB. is prime
    bn=:1+(10x^+&1)*666+10^+&3 NB. nth Belphegor number
    bp=:ip@bn#] NB. ns yielding prime Belphegor numbers
    NB. for example:
    pb 11 12 13 NB. prints 13
    pb i.14 NB. prints 0 13
    NB. alternatively, all at once:
    BP=:(1=#@q:)@(1+(10x^+&1)*666+10^+&3)#]

  3. Graham said

    J formatting via code blocks leaves something to be desired; apologies!

    ip=:1=#@q:                      NB. is prime
    bn=:1+(10x^+&1)*666+10^+&3      NB. nth Belphegor number
    bp=:ip@bn#]                     NB. ns yielding prime Belphegor numbers
    NB. for example:
    pb 11 12 13                     NB. prints 13
    pb i.14                         NB. prints 0 13
    NB. alternatively, all at once:
    BP=:(1=#@q:)@(1+(10x^+&1)*666+10^+&3)#]
    
  4. Lucas A. Brown said
    from itertools import count
    from your_math_library import isprime
    
    # Use gmpy if it's available.
    try: from gmpy import mpz; tens = mpz(1)
    except: tens = 1
    
    for n in count():
        belph = (tens * 1000 + 666) * tens * 10 + 1
        if isprime(belph): print n, "prime"
        elif n%100 == 0:   print n
        tens *= 10
    
  5. danaj said

    Perl, with progress indicators every 100. About a minute to get to 2623 so faster than gmpy or Pari/GP for me. It’ll be a long>/i> time getting to 28291. These numbers make easy BLS75-T5 proofs so the proof doesn’t add too much time. Replacing with is_prob_prime (BPSW) or is_prime (BPSW+one or more random base M-R tests) makes it run slightly faster, but between 2623 and 28291 they’d all be the same regardless.

    use ntheory qw/is_prime/;
    for (0 .. 100000) {
      say if is_provable_prime("1"."0"x$_."666"."0"x$_."1");
      say "  - $_" unless $_ % 100;
    }
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: