# Find all Belphegor primes

A Belphegor number is a number of the form $$\(10^{n+3}+666)*10^{n+1}+1\$$ (1{n zeroes}666{n zeroes}1) where $$\n\$$ is an non-negative integer. A Belphegor prime is a Belphegor number that is also prime.

The $$\n\$$ values of the first few Belphegor primes are 0, 13, 42, 506 (A232448)

Write a program that either:

• takes no input and outputs all Belphegor primes.
• takes a input $$\k\$$ and outputs the first $$\k\$$ Belphegor primes.

A reference python implementation can be found here.

# Rules

• You may output the $$\n\$$ value for the Belphegor prime instead of the prime itself
• You may use probabilistic primality tests as long as there is no known counter case.

# Scoring

This is so shortest bytes wins.

Inspired by The Most Evil Number - Numberphile

• what do you mean by "All"? there might be infinitely many... – J42161217 May 27 at 16:33
• @J42161217 by "All" I mean to write an program that does not stop and will eventually output all Belphegor primes. – Mukundan314 May 27 at 16:37
• Do you mean "$n\text{-th}$ value" by "$n$ value"? – Jonathan Frech May 27 at 22:48
• To be pedantic (and reading the fine print), the only known Belphegor primes are 16661 and 1000000000000066600000000000001. The rest of the numbers in the sequence are only probable primes. – Abigail May 28 at 10:48
• I find the criterion as long as there is no known counter case a bit shaky. It sounds like "it's fine as long as we can't tell". (For example, even if we don't know a counter case for a strong Baillie-PSW primality test to date, it is conjectured that there are infinitely many of them.) – Arnauld May 29 at 8:06

# Wolfram Language (Mathematica), 51 bytes

outputs the n value
"...program that does not stop and will eventually output all..."

PrimeQ[10^c*666+1+100^++c]~If~Print[c-2]~Do~{c,∞}


thanks to @DanTheMan for saving 4 bytes
and also to @mypronoun -7 bytes

• Using Do[...,∞] would be shorter than using While. Additionally, If can use infix syntax. – DanTheMan May 28 at 0:32
• 53 bytes: Try it online! – the default. May 28 at 5:30
• 51 bytes by rearranging 10^(2c+2) to 100^++c: Try it online! – the default. May 28 at 10:53
• @mypronounismonicareinstate good job! – J42161217 May 28 at 13:26

# Pyth, 17 bytes

.fP_sj666_B^TZQ0


Try it online!

Takes k as input and outputs the n corresponding to the first k Belphegor primes.

Explanation:

.fP_sj666_B^TZQ0
.f             Q0    Find the first k values of Z where the following is true,
starting at 0 and counting upwards.
^TZ      Raise 10 to the power of Z
         Convert to a string
_B          Pair with reversal
j666            Join with 666 in the middle
s                Convert to number
P_                 Check for primality.


# 05AB1E, 14 bytes

∞<ε0Xr×66Jû}ʒp


Outputs the infinite sequence.
Extremely slow due to the prime-check on large numbers, so times out before it even reaches the n=13 Belphegor prime on TIO..

Explanation:

∞             # Push an infinite positive list: [1,2,3,...]
<            # Decrease each by one to make it start at 0: [0,1,2,...]
ε           # Map each value to:
0          #  Push a 0
X         #  Push a 1
r        #  Reverse the stack order: [value, 0, 1] to [1, 0, value]
×       #  Repeat the 0 the value amount of times as string
66     #  Push 66
J    #  Join the values on the stack together: "10...066"
û   #  Palindromize it: "10...06660...01"
}ʒ          # After the map: filter the list by:
p         #  Check whether it's a prime number
# (after which the resulting list is output implicitly)


# Pyth, 22 bytes

.V0IP_h*+^T+3b666^Thbb


Try it online!

Implements the formula provided in the question. Prints the n values rather than the primes themselves.

Since this version (not surprisingly) times out on TIO, here is a version that prints all n values lower than the input: Try it online!

• Alternate 22 byte solution by doing string addition – Mukundan314 May 28 at 3:17

# JavaScript (Node.js), 71 bytes

A full program that prints Belphegor primes forever ... and takes forever to print them.

for(k=10n;;)for(d=n=666n*k+(k*=10n)*k+1n;n%--d||d<2n&&console.log(n););


Try it online!

### Commented

for(k = 10n;;)            // outer loop: start with k = 10 and loop forever
for(                    //   inner loop:
666n * k +          //       666 * k +
(k *= 10n) * k +    //       (10 * k)² +
1n;                 //       1
//     and update k to 10 * k
n % --d ||          //     decrement d until it divides n
d < 2n &&         //     if d is less than 2:
console.log(n); //       n is prime --> print it
);                      //


### JavaScript (Node.js), 176 bytes (non-competing)

A much faster version that uses a single iteration of the Miller-Rabin primality test.

for(k=10n;;)(n=666n*k+(k*=10n)*k+1n,~-(x=(g=(d,r,a)=>d?g(d/2n,d&1n?r*a%n:r,a*a%n):r)(d=n/(~-n&1n-n),1n,2n))&&~x+n?(g=d=>~d+n?~-(x=x*x%n)?~x+n&&g(d+d):1:1)(d):0)||console.log(n)


Try it online!

I guess it doesn't comply with the challenge rules since the test is likely to produce false-positives. It does however find the same 5 first terms as other answers.

# Python, 220 164 bytes

def a(k,s=set()):
for i in range(k):
p=1;n=(10**(i+3)+666)*10**-~i+1
for d in range(1,int(n**.5//1/2)):
p*=n%-~(d*2)>0
if~-p:break
return s


Simple prime search by checking modulus below the square root; fastened by skipping every even divisor.

Likely there's room for improvement, as it becomes incredibly slow for k > 10.

Edit: thanks to @JonathanAllan and @mathjunkie for ideas and sources. This update has heavy use of tweaks and bit-operations.

• I doubt that altering the golfiest Python prime identifying program would even print out the first one without stupid amounts of resources :) – Jonathan Allan May 27 at 22:35
• True->1; inline the ifs; use p*=n%((d*2)+1)>0 (maybe even p*=n%-~(d*2)>0?); (i+1)->-~i – Jonathan Allan May 27 at 22:38
• from math import* appears to be shorter. – Jonathan Frech May 27 at 22:46
• @JonathanAllan can you explain -~i? I'm not familiar with bitwise operators that much. How negating the bitwise negation of i` equals +1? – Zoltán Schmidt May 27 at 23:33
• @ZoltánSchmidt Take a look at this Python golfing tip – math junkie May 27 at 23:36