Следуя хорошей традиции таких вопросов, как « Найти наибольшее простое число, длина, сумма и произведение которого простое» , это вариант самой большой простой задачи.

вход

Ваш код не должен принимать никаких данных.

Определение

Мы говорим, что премьер pявляется , goodесли p-1имеет ровно 2различные простые множители.

Выход

Ваш код должен выводить абсолютную разницу между последовательными хорошими простыми числами qи, pтаким образом, |q-p|быть как можно большей и qнаименьшим хорошим простым, большим, чем p. Вы можете вывести любое количество хороших пар, и ваш последний результат будет принят за оценку.

пример

Последовательность первых 55 хороших простых чисел https://oeis.org/A067466 .

Гол

Ваша оценка просто |q-p| для пары хороших простых чисел, которые вы выводите.

Языки и библиотеки

Вы можете использовать любой язык или библиотеку, которая вам нравится (которая не была предназначена для этой задачи), кроме любых библиотечных функций для тестирования простоты или факторизации целых чисел. Тем не менее, для целей скоринга я буду запускать ваш код на моей машине, поэтому, пожалуйста, предоставьте четкие инструкции о том, как запустить его в Ubuntu.

Моя машина Время будет запущено на моей машине. Это стандартная установка Ubuntu на 8-гигабайтный 8-ядерный процессор AMD FX-8350. Это также означает, что мне нужно иметь возможность запускать ваш код.

Детали

• Я убью твой код через 2 минуты, если до этого ему не хватит памяти. Поэтому следует обязательно вывести что-то перед обрезкой.
• Вы не можете использовать какой-либо внешний источник простых чисел.
• Вы можете использовать вероятностные методы первичного тестирования, хотя Mego сказал мне, что с хорошими таблицами Миллер-Рабин может детерминистически тестировать до 341,550,071,728,321 (или даже выше). Смотрите также http://miller-rabin.appspot.com/ .

Лучшие записи, которые проверяют все целые числа из 1

• 756 котом в го
• 756 от El'endia Starman в Python
• 1932 Аднаном в C # (используя моно 3.2.8)
• 2640 по йети в Python (с использованием Pypy 4.01)
• 2754 Рето Коради в C ++
• 3486 Питер Тейлор на Яве
• 3900 примо в RPython (используя pypy 4.01)
• 4176 кодером на Яве

Лучшие записи, которые могут пропустить большое количество целых чисел, чтобы найти большой пробел

• 14226 Рето Коради в C ++
• 22596 по primo в RPython (используя pypy 4.01). Запись достигнута через 5 секунд!

RPython (PyPy 4.0.1), 4032

RPython - это ограниченное подмножество Python, которое можно преобразовать в C и затем скомпилировать с помощью RPython Toolchain. Его выраженная цель - помочь в создании языковых переводчиков, но его также можно использовать для компиляции простых программ.

Для компиляции загрузите текущий исходный код PyPy (PyPy 4.0.1) и выполните следующее:

$ pypy /pypy-4.0.1-src/rpython/bin/rpython --opt=3 good-primes.py

Полученный исполняемый файл будет назван good-primes-cили похож в текущем рабочем каталоге.

Замечания по реализации

Генератор простых чисел primes- это неограниченное сито Эратосфена, которое использует колесо, чтобы избежать любых кратных 2 , 3 , 5 или 7 . Он также вызывает себя рекурсивно для генерации следующего значения, которое будет использоваться для маркировки. Я вполне доволен этим генератором. Профилирование линий показывает, что две самые медленные линии:

37>      n += o
38>      if n not in sieve:

поэтому я не думаю, что есть много возможностей для улучшения, кроме, возможно, использования большего колеса.

Для проверки «добродетели» сначала из n-1 удаляются все факторы из двух , используя хитроумный хак, чтобы найти наибольшую степень двух, которая является делителем (n-1 & 1-n). Поскольку p-1 обязательно даже для любого простого числа p> 2 , отсюда следует, что 2 должно быть одним из различных простых факторов. То, что осталось, отправляется is_prime_powerфункции, которая делает то, что подразумевает ее имя. Проверка, является ли значение простой степенью, «почти свободна», так как это делается одновременно с проверкой первичности, с не более чем O (log _p n) операциями, где p - наименьший простой коэффициент n, Деление проб может показаться немного наивным, но по моим тестам это самый быстрый метод для значений менее 2 ³², Я немного экономлю, повторно используя колесо из сита. Особенно:

59>      while p*p < n:
60>        for o in offsets:

перебирая колесо длиной 48, p*p < n проверка будет пропущена тысячи раз по низкой низкой цене, не превышающей 48 дополнительных операций по модулю. Это также пропускает более 77% всех кандидатов, а не 50%, принимая только шансы.

Последние несколько выходов:

3588 (987417437 - 987413849) 60.469000s
3900 (1123404923 - 1123401023) 70.828000s
3942 (1196634239 - 1196630297) 76.594000s
4032 (1247118179 - 1247114147) 80.625000s
4176 (1964330609 - 1964326433) 143.047000s
4224 (2055062753 - 2055058529) 151.562000s

Код также является допустимым Python и должен достигать 3588 ~ 3900 при запуске с последним интерпретатором PyPy.

# primes less than 212
small_primes = [
    2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 37,
   41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89,
   97,101,103,107,109,113,127,131,137,139,149,151,
  157,163,167,173,179,181,191,193,197,199,211]

# pre-calced sieve of eratosthenes for n = 2, 3, 5, 7
# distances between sieve values, starting from 211
offsets = [
  10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6,
   6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4,
   2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6,
   4, 2, 4, 6, 2, 6, 4, 2, 4, 2,10, 2]

# tabulated, mod 105
dindices =[
  0,10, 2, 0, 4, 0, 0, 0, 8, 0, 0, 2, 0, 4, 0,
  0, 6, 2, 0, 4, 0, 0, 4, 6, 0, 0, 6, 0, 0, 2,
  0, 6, 2, 0, 4, 0, 0, 4, 6, 0, 0, 2, 0, 4, 2,
  0, 6, 6, 0, 0, 0, 0, 6, 6, 0, 0, 0, 0, 4, 2,
  0, 6, 2, 0, 4, 0, 0, 4, 6, 0, 0, 2, 0, 6, 2,
  0, 6, 0, 0, 4, 0, 0, 4, 6, 0, 0, 2, 0, 4, 8,
  0, 0, 2, 0,10, 0, 0, 4, 0, 0, 0, 2, 0, 4, 2]

def primes(start = 0):
  for n in small_primes[start:]: yield n
  pg = primes(6)
  p = pg.next()
  q = p*p
  sieve = {221: 13, 253: 11}
  n = 211
  while True:
    for o in offsets:
      n += o
      stp = sieve.pop(n, 0)
      if stp:
        nxt = n/stp
        nxt += dindices[nxt%105]
        while nxt*stp in sieve: nxt += dindices[nxt%105]
        sieve[nxt*stp] = stp
      else:
        if n < q:
          yield n
        else:
          sieve[q + dindices[p%105]*p] = p
          p = pg.next()
          q = p*p

def is_prime_power(n):
  for p in small_primes:
    if n%p == 0:
      n /= p
      while n%p == 0: n /= p
      return n == 1
  p = 211
  while p*p < n:
    for o in offsets:
      p += o
      if n%p == 0:
        n /= p
        while n%p == 0: n /= p
        return n == 1
  return n > 1

def main(argv):
  from time import time
  t0 = time()
  m = 0
  p = q = 7
  pgen = primes(3)

  for n in pgen:
    d = (n-1 & 1-n)
    if is_prime_power(n/d):
      p, q = q, n
      if q-p > m:
        m = q-p
        print m, "(%d - %d) %fs"%(q, p, time()-t0)

  return 0

def target(*args):
  return main, None

if __name__ == '__main__':
  from sys import argv
  main(argv)

RPython (PyPy 4.0.1), 22596

Это представление немного отличается от других, опубликованных до сих пор, в том, что оно не проверяет все хорошие простые числа, но вместо этого делает относительно большие скачки. Одним из недостатков этого является то, что сита не могут быть использованы ^{[Я исправлен?]} , Поэтому нужно полностью полагаться на тестирование на простоту, которое на практике немного медленнее. Существует также радостная среда между скоростью роста и количеством проверяемых значений каждый раз. Меньшие значения гораздо быстрее проверяются, но большие значения с большей вероятностью будут иметь большие пропуски.

Чтобы успокоить математических богов, я решил следовать последовательности, подобной Фибоначчи, имея следующую отправную точку как сумму двух предыдущих. Если после проверки 10 пар новые записи не найдены, скрипт переходит к следующей.

Последние несколько выходов:

6420 (12519586667324027 - 12519586667317607) 0.364000s
6720 (707871808582625903 - 707871808582619183) 0.721000s
8880 (626872872579606869 - 626872872579597989) 0.995000s
10146 (1206929709956703809 - 1206929709956693663) 4.858000s
22596 (918415168400717543 - 918415168400694947) 8.797000s

При компиляции используются 64-битные целые числа, хотя в некоторых местах предполагается, что два целых числа могут быть добавлены без переполнения, поэтому на практике можно использовать только 63. По достижении 62 значащих битов текущее значение уменьшается вдвое, чтобы избежать переполнения в вычислениях. В результате сценарий перебирает значения в диапазоне 2 ⁶⁰ - 2 ⁶² . Непревзойденная точность целых чисел также делает скрипт быстрее при интерпретации.

Следующий скрипт PARI / GP может использоваться для подтверждения этого результата:

isgoodprime(n) = isprime(n) && omega(n-1)==2

for(n = 918415168400694947, 918415168400717543, {
  if(isgoodprime(n), print(n" is a good prime"))
})

try:
  from rpython.rlib.rarithmetic import r_int64

  from rpython.rtyper.lltypesystem.lltype import SignedLongLongLong
  from rpython.translator.c.primitive import PrimitiveType

  # check if the compiler supports long long longs
  if SignedLongLongLong in PrimitiveType:

    from rpython.rlib.rarithmetic import r_longlonglong

    def mul_mod(a, b, m):
      return r_int64(r_longlonglong(a)*b%m)

  else:

    from rpython.rlib.rbigint import rbigint

    def mul_mod(a, b, m):
      biga = rbigint.fromrarith_int(a)
      bigb = rbigint.fromrarith_int(b)
      bigm = rbigint.fromrarith_int(m)

      return biga.mul(bigb).mod(bigm).tolonglong()


  # modular exponentiation b**e (mod m)
  def pow_mod(b, e, m):
    r = 1
    while e:
      if e&1: r = mul_mod(b, r, m)
      e >>= 1
      b = mul_mod(b, b, m)
    return r

except:

  import sys

  r_int64 = int
  if sys.maxint == 2147483647:
    mul_mod = lambda a, b, m: a*b%m
  else:
    mul_mod = lambda a, b, m: int(a*b%m)
  pow_mod = pow


# legendre symbol (a|m)
# note: returns m-1 if a is a non-residue, instead of -1
def legendre(a, m):
  return pow_mod(a, (m-1) >> 1, m)


# strong probable prime
def is_sprp(n, b=2):
  if n < 2: return False
  d = n-1
  s = 0
  while d&1 == 0:
    s += 1
    d >>= 1

  x = pow_mod(b, d, n)
  if x == 1 or x == n-1:
    return True

  for r in xrange(1, s):
    x = mul_mod(x, x, n)
    if x == 1:
      return False
    elif x == n-1:
      return True

  return False


# lucas probable prime
# assumes D = 1 (mod 4), (D|n) = -1
def is_lucas_prp(n, D):
  Q = (1-D) >> 2

  # n+1 = 2**r*s where s is odd
  s = n+1
  r = 0
  while s&1 == 0:
    r += 1
    s >>= 1

  # calculate the bit reversal of (odd) s
  # e.g. 19 (10011) <=> 25 (11001)
  t = r_int64(0)
  while s:
    if s&1:
      t += 1
      s -= 1
    else:
      t <<= 1
      s >>= 1

  # use the same bit reversal process to calculate the sth Lucas number
  # keep track of q = Q**n as we go
  U = 0
  V = 2
  q = 1
  # mod_inv(2, n)
  inv_2 = (n+1) >> 1
  while t:
    if t&1:
      # U, V of n+1
      U, V = mul_mod(inv_2, U + V, n), mul_mod(inv_2, V + mul_mod(D, U, n), n)
      q = mul_mod(q, Q, n)
      t -= 1
    else:
      # U, V of n*2
      U, V = mul_mod(U, V, n), (mul_mod(V, V, n) - 2 * q) % n
      q = mul_mod(q, q, n)
      t >>= 1

  # double s until we have the 2**r*sth Lucas number
  while r:
    U, V = mul_mod(U, V, n), (mul_mod(V, V, n) - 2 * q) % n
    q = mul_mod(q, q, n)
    r -= 1

  # primality check
  # if n is prime, n divides the n+1st Lucas number, given the assumptions
  return U == 0


# primes less than 212
small_primes = [
    2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 37,
   41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89,
   97,101,103,107,109,113,127,131,137,139,149,151,
  157,163,167,173,179,181,191,193,197,199,211]

# pre-calced sieve of eratosthenes for n = 2, 3, 5, 7
indices = [
    1, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
   53, 59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,
  107,109,113,121,127,131,137,139,143,149,151,157,
  163,167,169,173,179,181,187,191,193,197,199,209]

# distances between sieve values
offsets = [
  10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6,
   6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4,
   2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6,
   4, 2, 4, 6, 2, 6, 4, 2, 4, 2,10, 2]

bit_lengths = [
  0x00000000, 0x00000001, 0x00000003, 0x00000007,
  0x0000000F, 0x0000001F, 0x0000003F, 0x0000007F,
  0x000000FF, 0x000001FF, 0x000003FF, 0x000007FF,
  0x00000FFF, 0x00001FFF, 0x00003FFF, 0x00007FFF,
  0x0000FFFF, 0x0001FFFF, 0x0003FFFF, 0x0007FFFF,
  0x000FFFFF, 0x001FFFFF, 0x003FFFFF, 0x007FFFFF,
  0x00FFFFFF, 0x01FFFFFF, 0x03FFFFFF, 0x07FFFFFF,
  0x0FFFFFFF, 0x1FFFFFFF, 0x3FFFFFFF, 0x7FFFFFFF]

max_int = 2147483647


# returns the index of x in a sorted list a
# or the index of the next larger item if x is not present
# i.e. the proper insertion point for x in a
def binary_search(a, x):
  s = 0
  e = len(a)
  m = e >> 1
  while m != e:
    if a[m] < x:
      s = m
      m = (s + e + 1) >> 1
    else:
      e = m
      m = (s + e) >> 1
  return m


def log2(n):
  hi = n >> 32
  if hi:
    return binary_search(bit_lengths, hi) + 32
  return binary_search(bit_lengths, n)


# integer sqrt of n
def isqrt(n):
  c = n*4/3
  d = log2(c)

  a = d>>1
  if d&1:
    x = r_int64(1) << a
    y = (x + (n >> a)) >> 1
  else:
    x = (r_int64(3) << a) >> 2
    y = (x + (c >> a)) >> 1

  if x != y:
    x = y
    y = (x + n/x) >> 1
    while y < x:
      x = y
      y = (x + n/x) >> 1
  return x

# integer cbrt of n
def icbrt(n):
  d = log2(n)

  if d%3 == 2:
    x = r_int64(3) << d/3-1
  else:
    x = r_int64(1) << d/3

  y = (2*x + n/(x*x))/3
  if x != y:
    x = y
    y = (2*x + n/(x*x))/3
    while y < x:
      x = y
      y = (2*x + n/(x*x))/3
  return x


## Baillie-PSW ##
# this is technically a probabalistic test, but there are no known pseudoprimes
def is_bpsw(n):
  if not is_sprp(n, 2): return False

  # idea shamelessly stolen from Mathmatica's PrimeQ
  # if n is a 2-sprp and a 3-sprp, n is necessarily square-free
  if not is_sprp(n, 3): return False

  a = 5
  s = 2
  # if n is a perfect square, this will never terminate
  while legendre(a, n) != n-1:
    s = -s
    a = s-a
  return is_lucas_prp(n, a)


# an 'almost certain' primality check
def is_prime(n):
  if n < 212:
    m = binary_search(small_primes, n)
    return n == small_primes[m]

  for p in small_primes:
    if n%p == 0:
      return False

  # if n is a 32-bit integer, perform full trial division
  if n <= max_int:
    p = 211
    while p*p < n:
      for o in offsets:
        p += o
        if n%p == 0:
          return False
    return True

  return is_bpsw(n)


# next prime strictly larger than n
def next_prime(n):
  if n < 2:
    return 2

  # first odd larger than n
  n = (n + 1) | 1
  if n < 212:
    m = binary_search(small_primes, n)
    return small_primes[m]

  # find our position in the sieve rotation via binary search
  x = int(n%210)
  m = binary_search(indices, x)
  i = r_int64(n + (indices[m] - x))

  # adjust offsets
  offs = offsets[m:] + offsets[:m]
  while True:
    for o in offs:
      if is_prime(i):
        return i
      i += o


# true if n is a prime power > 0
def is_prime_power(n):
  if n > 1:
    for p in small_primes:
      if n%p == 0:
        n /= p
        while n%p == 0: n /= p
        return n == 1

    r = isqrt(n)
    if r*r == n:
      return is_prime_power(r)

    s = icbrt(n)
    if s*s*s == n:
      return is_prime_power(s)

    p = r_int64(211)
    while p*p < r:
      for o in offsets:
        p += o
        if n%p == 0:
          n /= p
          while n%p == 0: n /= p
          return n == 1

    if n <= max_int:
      while p*p < n:
        for o in offsets:
          p += o
          if n%p == 0:
            return False
      return True

    return is_bpsw(n)
  return False


def next_good_prime(n):
  n = next_prime(n)
  d = (n-1 & 1-n)
  while not is_prime_power(n/d):
    n = next_prime(n)
    d = (n-1 & 1-n)
  return n


def main(argv):
  from time import time
  t0 = time()

  if len(argv) > 1:
    n = r_int64(int(argv[1]))
  else:
    n = r_int64(7)

  if len(argv) > 2:
    limit = int(argv[2])
  else:
    limit = 10

  m = 0
  e = 1
  q = n
  try:
    while True:
      e += 1
      p, q = q, next_good_prime(q)
      if q-p > m:
        m = q-p
        print m, "(%d - %d) %fs"%(q, p, time()-t0)
        n, q = p, n+p
        if log2(q) > 61:
          q >>= 2
        e = 1
        q = next_good_prime(q)
      elif e > limit:
        n, q = p, n+p
        if log2(q) > 61:
          q >>= 2
        e = 1
        q = next_good_prime(q)
  except KeyboardInterrupt:
    pass
  return 0

def target(*args):
  return main, None

if __name__ == '__main__':
  from sys import argv
  main(argv)

Вероятно, 4032, смешанное сито Аткина-Бернштейна и "детерминированный" Миллер-Рабин

Факторизация колес и хорошие простые числа

Совершенно очевидно, что за исключением 2, 3 и 5, каждое простое число взаимно простое с 2 * 3 * 5 = 60. Существует 16 классов эквивалентности по модулю 60, которые взаимно просты с 60, поэтому любой тест на простоту должен проверять только те 16 случаев.

Однако, когда мы ищем «хорошие» простые числа, мы можем еще больше развести стадо. Если gcd(x, 60) = 1мы обнаружим, что в большинстве случаев gcd(x-1, 60)это либо 6, либо 10. Существует 6 значений, xдля которых оно равно 2:

17, 23, 29, 47, 53, 59

Поэтому мы можем заранее вычислить «хорошие» простые числа формы 2^a 3^b + 1и 2^a 5^b + 1объединить их в результат простого генератора, который рассматривает только 10% чисел как даже потенциальные простые числа.

Замечания по реализации

Так как у меня уже была Java-реализация сита Аткина-Бернштейна, и он уже использует колесо в качестве ключевого компонента, казалось естественным убрать ненужные спицы и адаптировать код. Первоначально я пытался использовать архитектуру производителя-потребителя, чтобы использовать 8 ядер, но управление памятью было слишком запутанным.

Чтобы проверить, является ли простое число «хорошим» простым, я использую «детерминированный» тест Миллера-Рабина (который на самом деле означает тест Миллера-Рабина, который кто-то другой предварительно проверил по списку, сгенерированному детерминистически). Это, конечно, можно переписать, чтобы использовать также Atkin-Bernstein, с некоторым кэшированием для охвата диапазонов, соответствующих sqrt, cbrt и т. Д., Но я не уверен, будет ли это улучшением (потому что это будет тестирование множества чисел, которые Мне не нужно проверять), так что это эксперимент на другой день.

На моем довольно старом компьютере это работает

987417437 - 987413849 = 3588
1123404923 - 1123401023 = 3900
1196634239 - 1196630297 = 3942
1247118179 - 1247114147 = 4032

ровно через 2 минуты, а затем

1964330609 - 1964326433 = 4176
2055062753 - 2055058529 = 4224
2160258917 - 2160254627 = 4290

в 3:10, 3:20 и 3:30 соответственно.

import java.util.*;

public class PPCG65876 {
    public static void main(String[] args) {
        long[] specials = genSpecials();
        int nextSpecialIdx = 0;
        long nextSpecial = specials[nextSpecialIdx];
        long p = 59;
        long bestGap = 2;

        for (long L = 1; true; L += B) {

            long[][] buf = new long[6][B >> 6];
            int[] Lmodqq = new int[qqtab.length];
            for (int i = 0; i < Lmodqq.length; i++) Lmodqq[i] = (int)(L % qqtab[i]);

            for (long[] arr : buf) Arrays.fill(arr, -1); // TODO Can probably get a minor optimisation by inverting this
            for (int[] parms : elliptic) traceElliptic(buf[parms[0]], parms[1], parms[2], parms[3] - L, parms[4], parms[5], Lmodqq, totients[parms[0]]);
            for (int[] parms : hyperbolic) traceHyperbolic(buf[parms[0]], parms[1], parms[2], parms[3] - L, Lmodqq, totients[parms[0]]);

            // We need to filter down to squarefree numbers.
            long pg_base = L * M;
            squarefreeMid(buf, invTotients, pg_base, 247, 1, 38);
            squarefreeMid(buf, invTotients, pg_base, 253, 1, 38);
            squarefreeMid(buf, invTotients, pg_base, 257, 1, 38);
            squarefreeMid(buf, invTotients, pg_base, 263, 1, 38);
            squarefreeMid(buf, invTotients, pg_base, 241, 0, 2);
            squarefreeMid(buf, invTotients, pg_base, 251, 0, 2);
            squarefreeMid(buf, invTotients, pg_base, 259, 0, 2);
            squarefreeMid(buf, invTotients, pg_base, 269, 0, 2);

            // Turn bitmasks into primes
            long[] page = new long[150000]; // TODO This can almost certainly be smaller
            long[] transpose = new long[6];
            for (int j = 0, off = 0; j < (B >> 6); j++) {
                // Reduce cache locality problems by transposing.
                for (int k = 0; k < 6; k++) transpose[k] = buf[k][j];
                for (long mask = 1L; mask != 0; mask <<= 1) {
                    for (int k = 0; k < 6; k++) {
                        if ((transpose[k] & mask) == 0) page[off++] = pg_base + totients[k];
                    }

                    pg_base += M;
                }
            }

            // Insert specials and look for gaps.
            for (long q : page) {
                if (q == 0) break;

                // Do we need to insert one or more specials?
                while (nextSpecial < q) {
                    if (nextSpecial - p > bestGap) {
                        bestGap = nextSpecial - p;
                        System.out.format("%d - %d = %d\n", nextSpecial, p, bestGap);
                    }

                    p = nextSpecial;
                    nextSpecial = specials[++nextSpecialIdx];
                }

                if (isGood(q)) {
                    if (q - p > bestGap) {
                        bestGap = q - p;
                        System.out.format("%d - %d = %d\n", q, p, bestGap);
                    }

                    p = q;
                }
            }

        }
    }

    static long[] genSpecials() {
        // 2^a 3^b + 1 or 2^a 5^b + 1
        List<Long> tmp = new LinkedList<Long>();
        for (long threes = 3; threes <= 4052555153018976267L; threes *= 3) {
            for (long t = threes; t > 0; t <<= 1) tmp.add(t + 1);
        }
        for (long fives = 5; fives <= 7450580596923828125L; fives *= 5) {
            for (long f = fives; f > 0; f <<= 1) tmp.add(f + 1);
        }

        // Filter down to primes
        Iterator<Long> it = tmp.iterator();
        while (it.hasNext()) {
            long next = it.next();
            if (next < 60 || next > 341550071728321L || !isPrime(next)) it.remove();
        }

        Collections.sort(tmp);
        long[] specials = new long[tmp.size()];
        for (int i = 0; i < tmp.size(); i++) specials[i] = tmp.get(i);

        return specials;
    }

    private static boolean isGood(long p) {
        long d = p - 1;
        while ((d & 1) == 0) d >>= 1;

        if (d == 1) return false;

        // Is d a prime power?
        if (d % 3 == 0 || d % 5 == 0) {
            // Because of the way the filters before this one work, nothing should reach here.
            throw new RuntimeException("Should be unreachable");
        }

        // TODO Is it preferable to reuse the Atkin-Bernstein code, caching pages which correspond
        // to the possible power candidates?
        if (isPrime(d)) return true;
        for (int a = (d % 60 == 1 || d % 60 == 49) ? 2 : 3; (1L << a) < d; a++) {
            long r = (long)(0.5 + Math.pow(d, 1. / a));
            if (d == (long)(0.5 + Math.pow(r, a)) && isPrime(r)) return true;
        }

        return false;
    }

    /*---------------------------------------------------
               Deterministic Miller-Rabin
    ---------------------------------------------------*/
    public static boolean isPrime(int x) {
        // See isPrime(long). We pick bases which are known to work for the entire range of int.
        // Special case for the bases.
        if (x == 2 || x == 7 || x == 61) return true;

        int d = x - 1;
        int s = 0;
        while ((d & 1) == 0) { s++; d >>= 1; } // TODO Can be optimised

        if (!isSPRP(2, d, s, x)) return false;
        if (!isSPRP(7, d, s, x)) return false;
        if (!isSPRP(61, d, s, x)) return false;
        return true;
    }

    private static boolean isSPRP(int b, int d, int s, int x /* == d << s */) {
        int l = modPow(b, d, x);
        if (l == 1 || l == x - 1) return true;
        for (int r = 1; r < s; r++) {
            l = modPow(l, 2, x);
            if (l == x - 1) return true;
            if (l == 1) return false;
        }

        return false;
    }

    public static int modPow(int a, int b, int c) {
        int accum = 1;
        while (b > 0) {
            if ((b & 1) == 1) accum = (int)(accum * (long)a % c);
            a = (int)(a * (long)a % c);
            b >>= 1;
        }
        return accum;
    }

    public static boolean isPrime(long x) {
        if (x < Integer.MAX_VALUE) return isPrime((int)x);

        long d = x - 1;
        int s = 0;
        while ((d & 1) == 0) { s++; d >>= 1; } // TODO Can be optimised

        // If b^d == 1 (mod x) or (b^d)^(2^r) == -1 (mod x) for some r < s then we pass for base b.
        // We select bases according to Jaeschke, Gerhard (1993), "On strong pseudoprimes to several bases", Mathematics of Computation 61 (204): 915–926, doi:10.2307/2153262
        // TODO Would it be better to use a set of 5 bases from http://miller-rabin.appspot.com/ ?
        if (!isSPRP(2, d, s, x)) return false;
        if (!isSPRP(3, d, s, x)) return false;
        if (!isSPRP(5, d, s, x)) return false;
        if (!isSPRP(7, d, s, x)) return false;
        if (x < 3215031751L) return true;
        if (!isSPRP(11, d, s, x)) return false;
        if (x < 2152302898747L) return true;
        if (!isSPRP(13, d, s, x)) return false;
        if (x < 3474749660383L) return true;
        if (!isSPRP(17, d, s, x)) return false;
        if (x < 341550071728321L) return true;

        throw new IllegalArgumentException("Overflow");
    }

    private static boolean isSPRP(long b, long d, int s, long x /* == d << s */) {
        if (b * (double)x > Long.MAX_VALUE) throw new IllegalArgumentException("Overflow"); // TODO Work out more precise page bounds

        long l = modPow(b, d, x);
        if (l == 1 || l == x - 1) return true;
        for (int r = 1; r < s; r++) {
            l = modPow(l, 2, x);
            if (l == x - 1) return true;
            if (l == 1) return false;
        }

        return false;
    }

    /**
     * Computes a^b (mod c). We assume c &lt; 2^62.
     */
    public static long modPow(long a, long b, long c) {
        long accum = 1;
        while (b > 0) {
            if ((b & 1) == 1) accum = prodMod(accum, a, c);
            a = prodMod(a, a, c);
            b >>= 1;
        }
        return accum;
    }

    /**
     * Computes a*b (mod c). We assume c &lt; 2^62.
     */
    private static long prodMod(long a, long b, long c) {
        // The naive product would require 128-bit integers.

        // Consider a = (A << 32) + B, b = (C << 31) + D. Different shifts chosen deliberately.
        // Then ab = (AC << 63) + (AD << 32) + (BC << 31) + BD with intermediate values remaining in 63 bits.
        long AC = (a >> 32) * (b >> 31) % c;
        long AD = (a >> 32) * (b & ((1L << 31) - 1)) % c;
        long BC = (a & ((1L << 32) - 1)) * (b >> 31) % c;
        long BD = (a & ((1L << 32) - 1)) * (b & ((1L << 31) - 1)) % c;

        long t = AC;
        for (int i = 0; i < 31; i++) {
            t = (t + t) % c;
        }
        // t = (AC << 31)
        t = (t + AD) % c;
        t = (t + t) % c;
        t = (t + BC) % c;
        // t = (AC << 32) + (AD << 1) + BC
        for (int i = 0; i < 31; i++) {
            t = (t + t) % c;
        }
        // t = (AC << 63) + (AD << 32) + (BC << 31)
        return (t + BD) % c;
    }

    /*---------------------------------------------------
                      Atkin-Bernstein
    ---------------------------------------------------*/
    // Page size.
    private static final int B = 1001 << 6;
    // Wheel modulus for sharding between binary quadratic forms.
    private static final int M = 60;

    // Squares of primes 5 < q < 240
    private static final int[] qqtab = new int[] {
        49, 121, 169, 289, 361, 529, 841, 961, 1369, 1681, 1849, 2209, 2809,
        3481, 3721, 4489, 5041, 5329, 6241, 6889, 7921, 9409, 10201, 10609, 11449, 11881,
        12769, 16129, 17161, 18769, 19321, 22201, 22801, 24649, 26569, 27889, 29929, 32041, 32761,
        36481, 37249, 38809, 39601, 44521, 49729, 51529, 52441, 54289, 57121
    };
    // If a_i == q^{-2} (mod 60) is the reciprocal of qq[i], qq60tab[i] = qq[i] + (1 - a_i * qq[i]) / 60
    private static int[] qq60tab = new int[] {
        9, 119, 31, 53, 355, 97, 827, 945, 251, 1653, 339, 405, 515,
        3423, 3659, 823, 4957, 977, 6137, 1263, 7789, 1725, 10031, 1945, 2099, 11683,
        2341, 2957, 16875, 3441, 18999, 21831, 22421, 4519, 4871, 5113, 5487, 31507, 32215,
        35873, 6829, 7115, 38941, 43779, 9117, 9447, 51567, 9953, 56169
    };

    /**
     * Produces a set of parameters for traceElliptic to find solutions to ax^2 + cy^2 == d (mod M).
     * @param d The target residue.
     * @param a Binary quadratic form parameter.
     * @param c Binary quadratic form parameter.
     */
    private static List<int[]> initElliptic(final int[] invTotients, final int d, final int a, final int c) {
        List<int[]> rv = new ArrayList<int[]>();

        // The basic idea is that we maintain an invariant of the form
        //     M k = a x^2 + c y^2 - d
        // Therefore we increment x in steps F such that
        //     a((x + F)^2 - x^2) == 0 (mod M)
        // and similarly for y in steps G.
        int F = computeIncrement(a, M), G = computeIncrement(c, M);
        for (int f = 1; f <= F; f++) {
            for (int g = 1; g <= G; g++) {
                if ((a*f*f + c*g*g - d) % M == 0) {
                    rv.add(new int[] { invTotients[d], (2*f + F)*a*F/M, (2*g + G)*c*G/M, (a*f*f + c*g*g - d)/M, 2*a*F*F/M, 2*c*G*G/M });
                }
            }
        }

        return rv;
    }

    private static int computeIncrement(int a, int M) {
        // Find smallest F such that M | 2aF and M | aF^2
        int l = M / gcd(M, 2 * a);
        for (int F = l; true; F += l) {
            if (a*F*F % M == 0) return F;
        }
    }

    public static int gcd(int a, int b) {
        while (b != 0) {
            int t = b;
            b = a % b;
            a = t;
        }

        return a;
    }

    // NB This is generalised somewhat from primegen's implementation.
    private static void traceElliptic(final long[] buf, int x, int y, long start, final int cF2, final int cG2, final int[] Lmodqq, final int d) {
        // Bring the annular segment into the range of ints.
        start += 1000000000;
        while (start < 0) {
            start += x;
            x += cF2;
        }
        start -= 1000000000;
        int i = (int)start;

        while (i < B) {
            i += x;
            x += cF2;
        }

        while (true) {
            x -= cF2;
            if (x <= cF2 >> 1) {
                // It makes no sense that doing this in here should perform well, but empirically it does much better than
                // only eliminating the squares once.
                squarefreeTiny(buf, Lmodqq, d);
                return;
            }
            i -= x;

            while (i < 0) {
                i += y;
                y += cG2;
            }

            int i0 = i, y0 = y;
            while (i < B) {
                buf[i >> 6] ^= 1L << i;
                i += y;
                y += cG2;
            }
            i = i0;
            y = y0;
        }
    }

    // This only handles 3x^2 - y^2, and is closer to a direct port of primegen.
    private static void traceHyperbolic(final long[] a, int x, int y, long start, final int[] Lmodqq, final int d) {
        x += 5;
        y += 15;

        // Bring the segment into the range of ints.
        start += 1000000000;
        while (start < 0) {
            start += x;
            x += 10;
        }
        start -= 1000000000;
        int i = (int)start;

        while (i < 0) {
            i += x;
            x += 10;
        }

        while (true) {
            x += 10;
            while (i >= B) {
                if (x <= y) {
                    squarefreeTiny(a, Lmodqq, d);
                    return;
                }
                i -= y;
                y += 30;
            }

            int i0 = i, y0 = y;
            while (i >= 0 && y < x) {
                a[i >> 6] ^= 1L << i;
                i -= y;
                y += 30;
            }
            i = i0 + x - 10;
            y = y0;
        }
    }

    private static void squarefreeTiny(final long[] a, final int[] Lmodqq, final int d) {
        for (int j = 0; j < qqtab.length; ++j) {
            int qq = qqtab[j];
            int k = qq - 1 - ((Lmodqq[j] + qq60tab[j] * d - 1) % qq);
            while (k < B) {
                a[k >> 6] |= 1L << k;
                k += qq;
            }
        }
    }

    private static void squarefreeMid(long[][] buf, int[] invTotients, final long base, int q, int dqq, int di) {
        int qq = q * q;
        q = M * q + (M * M / 4);

        while (qq < M * B) {
            int i = qq - (int)(base % qq);
            if ((i & 1) == 0) i += qq;

            if (i < M * B) {
                int qqhigh = ((qq / M) << 1) + dqq;
                int ilow = i % M;
                int ihigh = i / M;
                while (ihigh < B) {
                    int n = invTotients[ilow];
                    if (n >= 0) buf[n][ihigh >> 6] |= 1L << ihigh;

                    ilow += di;
                    ihigh += qqhigh;
                    if (ilow >= M) {
                        ilow -= M;
                        ihigh += 1;
                    }
                }
            }

            qq += q;
            q += M * M / 2;
        }

        squarefreebig(buf, invTotients, base, q, qq);
    }

    private static void squarefreebig(long[][] buf, int[] invTotients, final long base, int q, long qq) {
        long bound = base + M * B;
        while (qq < bound) {
            long i = qq - (base % qq);
            if ((i & 1) == 0) i += qq;

            if (i < M * B) {
                int pos = (int)i;
                int n = invTotients[pos % M];
                if (n >= 0) {
                    int ihigh = pos / M;
                    buf[n][ihigh >> 6] |= 1L << ihigh;
                }
            }

            qq += q;
            q += M * M / 2;
        }
    }

    // The relevant totients of M - those which only have one forced prime factor.
    static final int[] totients = new int[] { 17, 23, 29, 47, 53, 59 };
    private static final int[] invTotients;
    // Parameters for tracing the hyperbolic BQF used for 59+60Z.
    private static final int[][] hyperbolic = new int[][] {
        {5, 1, 2, -1}, {5, 1, 8, -2}, {5, 1, 22, -9}, {5, 1, 28, -14}, {5, 4, 7, -1}, {5, 4, 13, -3}, {5, 4, 17, -5}, {5, 4, 23, -9},
        {5, 5, 4, 0}, {5, 5, 14, -3}, {5, 5, 16, -4}, {5, 5, 26, -11}, {5, 6, 7, 0}, {5, 6, 13, -2}, {5, 6, 17, -4}, {5, 6, 23, -8},
        {5, 9, 2, 3}, {5, 9, 8, 2}, {5, 9, 22, -5}, {5, 9, 28, -10}, {5, 10, 1, 4}, {5, 10, 11, 2}, {5, 10, 19, -2}, {5, 10, 29, -10}
    };

    // Parameters for tracing the elliptic BQFs used for all totients except 11 and 59.
    private static final int[][] elliptic;
    static {
        invTotients = new int[M];
        Arrays.fill(invTotients, -1);
        for (int i = 0; i < totients.length; i++) invTotients[totients[i]] = i;

        // Calculate the parameters for tracing the elliptic BQFs from a table of the BQF used for each totient.
        // E.g. for 17+60Z we use 5x^2 + 3y^2.
        int[][] bqfs = new int[][] {
            {17, 5, 3}, {23, 5, 3}, {29, 4, 1}, {47, 5, 3}, {53, 5, 3}
        };
        List<int[]> parmSets = new ArrayList<int[]>();
        for (int[] bqf : bqfs) parmSets.addAll(initElliptic(invTotients, bqf[0], bqf[1], bqf[2]));
        elliptic = parmSets.toArray(new int[0][]);
    }
}