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.
Programming Praxis 
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.”);
}
}
}
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)#]
J formatting via code blocks leaves something to be desired; apologies!
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 *= 10Perl, 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; }One likely, and relatively easy way to exocise Belphegore, is by having the possessed person work nonstop. Belphegore is the demon of sloth and laziness. He will likely find a new host at some point in this process.