Proth Primes And Sierpinski Numbers

April 6, 2018

One of the neat things about mathematics is that, throughout history, many contributions to mathematics have been made by self-taught amateur mathematicians. About 150 years ago, a French farmer named François Proth (1852–1879), who lived near Verdun, proved this theorem:

Let N = k 2n + 1 with k odd and 2n > k. Choose an integer a so that the Jacobi symbol (a, N) is -1. Then N is prime if and only if a(N−1)/2 ≡ -1 (mod N).

Primes of that form are known as Proth primes. Testing for a Proth prime is easy: choose a ∈ {3, 5, 7, 17} so that the Jacobi symbol is -1 and perform the modular multiplication. The game that recreational mathematicians play is to fix k and iterate n and see which k are fertile in producing primes; for instance, k = 12909 produces 81 primes before n = 53118.

Sierpinski numbers have the same form as Proth numbers, but “opposite” in the sense that there are no primes for a given k. The smallest known Sierpinski k is 78557; the only smaller k for which the Sierpinski character is not known are 10223, 21181, 22699, 24737, 55459 and 67607.

Your task is to write a program that identifies Proth primes and use it to find fertile k and Sierpinski k. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

Posted by programmingpraxis
Filed in Exercises
3 Comments »

3 Responses to “Proth Primes And Sierpinski Numbers”

  2. Daniel said
    May 24, 2018 at 9:30 PM

    @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
    May 24, 2018 at 9:33 PM

    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 comment