Proth Primes And Sierpinski Numbers

April 6, 2018

We begin with a function proth? that is #t if any of 3, 5, 7 or 17 identifies a prime, and #f otherwise:

(define (proth? n)
  (let loop ((a '(3 5 7 17)))
    (if (null? a) #f
      (if (not (= (jacobi (car a) n) -1)) (loop (cdr a))
        (if (= (expm (car a) (/ (- n 1) 2) n) (- n 1)) #t
          (loop (cdr a)))))))
(define (proth k)
  (do ((k2n (* k 2) (* k2n 2)) (n 1 (+ n 1))) (#f)
    (when (proth? (+ k2n 1)) (display n) (newline))))

And here’s the beginning of the output for k = 12909:

> (proth 12909)
1
2
5
7
11
14
17
18
22
30
39
58
62
77
85
91
99
109
125
131
149
171
185
215
223
253
267
285
287
299
310
337
347
391
639
641
685
767
771
781
887
CTRL-C

Here we see a Sierpinski number:

> (proth 78557)
... wait forever ...

You can run the program at https://ideone.com/e6FnQS.

Advertisements

Pages: 1 2

3 Responses to “Proth Primes And Sierpinski Numbers”

  1. Daniel said

    @programmingpraxis, the problem states that 2^n > k, but the output of your solution includes n for which the inequality does not hold when k = 12909. For example, the first output is n = 1. In that case 2^n <= k and N = (53118 * 2^1) + 1 = 106237, which is not prime.

  2. Daniel said

    Here’s a solution in Python.

    from itertools import count
    import sys
    
    def jacobi(a, m):
      if m % 2 == 0: raise RuntimeError()
      t = 1
      a %= m
      while a != 0:
        while a % 2 == 0:
          a //= 2
          if m % 8 in (3,5): t = -t
        a, m = m, a
        if a % 4 == m % 4 == 3: t = -t
        a %= m
      return t if m == 1 else 0
    
    def proth(k, n):
      n2 = 2 ** n
      if n2 <= k: raise RuntimeError()
      N = k * n2 + 1
      for a in (3,5,7,17):
        j = jacobi(a, N)
        if j == -1: break
      return pow(a, (N-1) // 2, N) == N - 1
    
    if __name__ == "__main__":
      if len(sys.argv) != 2:
        sys.stderr.write("Usage: {} <k>\n".format(sys.argv[0]))
        sys.exit(1)
      k = int(sys.argv[1])
      idx = 0
      for n in count(k.bit_length()):
        is_prime = proth(k, n)
        if is_prime:
          idx += 1
          print(idx, n)
    

    Examples:

    $ python proth.py 12909
    1 14
    2 17
    3 18
    4 22
    5 30
    6 39
    7 58
    8 62
    9 77
    10 85
    11 91
    12 99
    13 109
    14 125
    15 131
    16 149
    ...
    
    $ python proth.py 78557
    ...
    

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 )

w

Connecting to %s

%d bloggers like this: