## Belphegor Primes

### October 31, 2014

We use the primality tester from a previous exercise and write this function to calculate Belphegor numbers:

```(define (belphegor x)   (let ((n (+ (expt 10 (+ x x 4))               (* 666 (expt 10 (+ x 1)))               1)))     n))```

Then we can list the Belphegor primes, though be aware that the Belphegor numbers quickly grow large, so the calculation is quite slow:

```> (do ((n 0 (+ n 1))) (#f)     (when (prime? (belphegor n))       (display n) (newline))) 0 13 42 506 ...```

You can run the program at http://programmingpraxis.codepad.org/5hBfAjlj. Happy Halloween!

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;
}
```