Mengikuti tradisi pertanyaan yang bagus seperti Temukan prime terbesar yang panjang, jumlah dan produknya prima , ini merupakan varian dari tantangan prime terbesar.

Memasukkan

Kode Anda tidak boleh mengambil input apa pun.

Definisi

Kami mengatakan prima padalah goodjika p-1memiliki 2faktor prima yang persis berbeda.

Keluaran

Kode Anda harus menampilkan perbedaan absolut antara bilangan prima baik berturut-turut qdan psehingga |q-p|sebesar mungkin dan qmerupakan bilangan prima terkecil terkecil lebih besar dari p. Anda dapat menampilkan sejumlah pasangan yang baik dan hasil terakhir Anda akan dianggap sebagai skor.

Contoh

Urutan dari 55 bilangan prima yang baik pertama adalah https://oeis.org/A067466 .

Skor

Skor Anda hanya |q-p|untuk pasangan bilangan prima bagus yang Anda hasilkan.

Bahasa dan perpustakaan

Anda dapat menggunakan bahasa atau pustaka yang Anda suka (yang tidak dirancang untuk tantangan ini) kecuali untuk fungsi pustaka apa pun untuk pengujian primality atau bilangan bulat anjak piutang. Namun, untuk keperluan penilaian saya akan menjalankan kode Anda di mesin saya jadi tolong berikan instruksi yang jelas untuk bagaimana menjalankannya di Ubuntu.

Mesin Saya Pengaturan waktu akan dijalankan pada mesin saya. Ini adalah instalasi standar Ubuntu pada Prosesor Delapan-Core AMD 8GB 8GB AMD. Ini juga berarti saya harus dapat menjalankan kode Anda.

Detail

• Saya akan membunuh kode Anda setelah 2 menit kecuali mulai kehabisan memori sebelum itu. Karena itu harus memastikan untuk mengeluarkan sesuatu sebelum dipotong.
• Anda tidak boleh menggunakan sumber bilangan prima eksternal.
• Anda dapat menggunakan metode pengujian prima probabilistik meskipun saya diberitahu oleh Mego bahwa dengan tabel yang baik, Miller-Rabin dapat menguji hingga 341.550.071.728.321 (atau bahkan lebih tinggi) secara deterministik. Lihat juga http://miller-rabin.appspot.com/ .

Entri terbaik yang memeriksa semua bilangan bulat dari 1

• 756 oleh kucing di Go
• 756 oleh El'endia Starman dengan Python
• 1932 oleh Adnan di C # (menggunakan mono 3.2.8)
• 2640 oleh yeti di Python (menggunakan pypy 4.01)
• 2754 oleh Reto Koradi di C ++
• 3486 oleh Peter Taylor di Jawa
• 3900 oleh primo dalam RPython (menggunakan pypy 4.01)
• 4176 oleh The Coder in Java

Entri terbaik yang dapat melompati sejumlah besar bilangan bulat untuk menemukan celah besar

• 14226 oleh Reto Koradi di C ++
• 22596 oleh primo dalam RPython (menggunakan pypy 4.01). Rekor tercapai setelah 5 detik!

RPython (PyPy 4.0.1), 4032

RPython adalah subset terbatas dari Python, yang dapat diterjemahkan ke C dan kemudian dikompilasi menggunakan RPython Toolchain. Tujuannya adalah untuk membantu dalam menciptakan penerjemah bahasa, tetapi juga dapat digunakan untuk menyusun program-program sederhana.

Untuk mengkompilasi, unduh sumber PyPy saat ini (PyPy 4.0.1), dan jalankan yang berikut:

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

Eksekusi yang dihasilkan akan diberi nama good-primes-catau serupa di direktori kerja saat ini.

Catatan Implementasi

Generator bilangan prima primesadalah Saringan Eratosthenes tanpa batas, yang menggunakan roda untuk menghindari kelipatan 2 , 3 , 5 , atau 7 . Itu juga menyebut dirinya secara rekursif untuk menghasilkan nilai berikutnya yang digunakan untuk menandai. Saya cukup puas dengan generator ini. Profil garis mengungkapkan bahwa dua baris paling lambat adalah:

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

jadi saya tidak berpikir ada banyak ruang untuk perbaikan, selain mungkin menggunakan roda yang lebih besar.

Untuk pemeriksaan "kebaikan", pertama semua faktor dua dihapus dari n-1 , menggunakan peretasan-twiddling untuk menemukan kekuatan terbesar dari dua yang merupakan pembagi (n-1 & 1-n). Karena p-1 harus bahkan untuk setiap prima p> 2 , maka 2 harus menjadi salah satu faktor prima yang berbeda. Apa yang tersisa dikirim ke is_prime_powerfungsi, yang sesuai dengan namanya. Memeriksa apakah nilai adalah kekuatan utama adalah "hampir bebas", karena hal itu dilakukan bersamaan dengan pemeriksaan primality, dengan paling O (log _p n) operasi, di mana p adalah faktor prima terkecil dari n. Divisi percobaan mungkin tampak naif sedikit, tapi dengan pengujian saya itu adalah metode tercepat untuk nilai-nilai kurang dari 2 ³² . Saya menghemat sedikit dengan menggunakan kembali roda dari ayakan. Khususnya:

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

dengan mengulangi roda panjang 48, p*p < ncek akan dilewati ribuan kali, pada harga rendah, rendah tidak lebih dari 48 operasi modulo tambahan. Itu juga melompati lebih dari 77% dari semua kandidat, daripada 50% dengan hanya mengambil peluang.

Beberapa keluaran terakhir adalah:

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

Kode juga Python yang valid, dan harus mencapai 3588 ~ 3900 ketika dijalankan dengan juru bahasa PyPy baru-baru ini.

# 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

Kiriman ini sedikit berbeda dari yang lain yang diposting sejauh ini, dalam hal ini tidak memeriksa semua bilangan prima yang baik, tetapi malah membuat lompatan yang relatif besar. Salah satu kelemahan dari melakukan ini adalah bahwa saringan tidak dapat digunakan ^{[saya berdiri dikoreksi?]} , Jadi kita harus bergantung sepenuhnya pada pengujian primality yang dalam praktiknya sedikit lebih lambat. Ada juga media bahagia yang bisa ditemukan antara tingkat pertumbuhan, dan jumlah nilai yang diperiksa setiap kali. Nilai yang lebih kecil jauh lebih cepat untuk diperiksa, tetapi nilai yang lebih besar lebih cenderung memiliki kesenjangan yang lebih besar.

Untuk menenangkan dewa-dewa matematika, saya telah memutuskan untuk mengikuti urutan Fibonacci, memiliki titik awal berikutnya sebagai jumlah dari dua sebelumnya. Jika tidak ada catatan baru ditemukan setelah memeriksa 10 pasangan, skrip bergerak pada yang berikutnya.

Beberapa keluaran terakhir adalah:

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

Ketika dikompilasi, integer 64-bit digunakan, meskipun diasumsikan di beberapa tempat bahwa dua integer dapat ditambahkan tanpa overflow, jadi dalam praktiknya hanya 63 yang dapat digunakan. Setelah mencapai 62 bit signifikan, nilai saat ini dibelah dua, untuk menghindari kelebihan dalam perhitungan. Hasilnya adalah bahwa skrip mengocok nilai-nilai pada rentang 2 ⁶⁰ - 2 ⁶² . Tidak melebihi presisi integer asli juga membuat skrip lebih cepat ketika ditafsirkan.

Skrip PARI / GP berikut dapat digunakan untuk mengonfirmasi hasil ini:

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)

Mungkin 4032, campuran Atkin-Bernstein dan "deterministik" Miller-Rabin

Faktorisasi roda dan bilangan prima yang baik

Sangat jelas bahwa dengan pengecualian 2, 3, dan 5, setiap prime adalah koprime hingga 2 * 3 * 5 = 60. Ada 16 kelas ekivalen modulo 60 yang coprime hingga 60, jadi setiap tes primality hanya perlu memeriksa 16 kasus.

Namun, ketika kita mencari bilangan prima yang "baik", kita bisa semakin menipiskan kawanan. Jika gcd(x, 60) = 1, kami menemukan bahwa dalam kebanyakan kasus gcd(x-1, 60)adalah 6 atau 10. Ada 6 nilai xyang 2:

17, 23, 29, 47, 53, 59

Oleh karena itu kita dapat melakukan precompute bilangan prima "baik" dari formulir 2^a 3^b + 1dan 2^a 5^b + 1menggabungkannya menjadi hasil dari generator prima yang hanya menganggap 10% dari bilangan sebagai bilangan prima potensial .

Catatan implementasi

Karena saya sudah memiliki implementasi Java dari saringan Atkin-Bernstein yang tergeletak di sekitar, dan yang sudah menggunakan roda sebagai komponen utama, rasanya wajar untuk mengeluarkan jari-jari yang tidak perlu dan mengadaptasi kode. Saya awalnya mencoba menggunakan arsitektur produsen-konsumen untuk mengeksploitasi 8 core, tetapi manajemen memori terlalu berantakan.

Untuk menguji apakah suatu prima adalah prima "baik", saya menggunakan uji Miller-Rabin "deterministik" (yang benar-benar berarti uji Miller-Rabin yang telah diperiksa orang lain terhadap daftar yang dihasilkan secara deterministik). Ini tentu saja dapat ditulis ulang untuk juga menggunakan Atkin-Bernstein, dengan beberapa caching untuk mencakup rentang yang sesuai dengan sqrt, cbrt, dll., Tetapi saya tidak yakin apakah itu akan menjadi peningkatan (karena akan menguji banyak angka yang Saya tidak perlu menguji), jadi itu percobaan untuk hari lain.

Di komputer saya yang cukup lama ini berjalan ke

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

dalam waktu 2 menit, dan kemudian

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

masing-masing pada 3:10, 3:20, dan 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][]);
    }
}