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!

Advertisements

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 )

Google photo

You are commenting using your Google 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 )

Connecting to %s

%d bloggers like this: