# 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... May 27, 2020 at 16:33
• @J42161217 by "All" I mean to write an program that does not stop and will eventually output all Belphegor primes. May 27, 2020 at 16:37
• Do you mean "$n\text{-th}$ value" by "$n$ value"? May 27, 2020 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. May 28, 2020 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.) May 29, 2020 at 8:06

# 05AB1E, 1410 9 bytes

∞°66+€ûʒ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,4,5,...]
°         # Take 10 to the power each: [10,100,1000,10000,100000,...]
66+      # Add 66 to each: [76,166,1066,10066,100066,...]
€û    # Palindromize each: [767,16661,1066601,100666001,10006660001,...]
ʒ   # Filter the list by:
p  #  Check whether it's a prime number
#  (minor note: 767 is NOT a prime)
# (after which the resulting list is output implicitly)


# 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. May 28, 2020 at 0:32
• 53 bytes: Try it online! May 28, 2020 at 5:30
• 51 bytes by rearranging 10^(2c+2) to 100^++c: Try it online! May 28, 2020 at 10:53
• @mypronounismonicareinstate good job! May 28, 2020 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.


# 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 May 28, 2020 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.

# Vyxal, 108 bitsv2, 13.5 bytes

Þ:'3+↵666+n›↵*›æ


Try it Online!

Using the T flag sets the interpreter to time out after 1 minute :P

• Porting my 05AB1E answer, is 11.375 bytes (91 bits), and possibly even shorter if you know how to prevent casting to/from strings to numbers for the palindromizing: try it online. Oct 2 at 15:05

# 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 :) May 27, 2020 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 May 27, 2020 at 22:38
• from math import* appears to be shorter. May 27, 2020 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? May 27, 2020 at 23:33
• @ZoltánSchmidt Take a look at this Python golfing tip May 27, 2020 at 23:36