## 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.

Pages: 1 2

### 3 Responses to “Proth Primes And Sierpinski Numbers”

2. 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.

3. 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))
sys.exit(1)
k = int(sys.argv)
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
...
```