Zgodnie z piękną tradycją pytań, takich jak Znajdź największą liczbę pierwszą, której długość, suma i iloczyn jest liczbą pierwszą , jest to wariant największego największego wyzwania.

Wkład

Twój kod nie powinien przyjmować żadnych danych wejściowych.

Definicja

Mówimy głównym pjest goodjeśli p-1ma dokładnie 2odmienne czynniki pierwszych.

Wydajność

Kod powinien wyjściowa bezwzględna różnica między kolejnymi liczbami pierwszymi dobrymi qi ptak, że |q-p|jest tak duży, jak to możliwe i qjest najmniejszym dobry prime większa niż p. Możesz wydać dowolną liczbę dobrych par, a twój ostatni wynik zostanie wykorzystany jako wynik.

Przykład

Sekwencja pierwszych 55 dobrych liczb pierwszych to https://oeis.org/A067466 .

Wynik

Twój wynik jest po prostu |q-p|za parę dobrych liczb pierwszych, które wyprowadzasz.

Języki i biblioteki

Możesz użyć dowolnego języka lub biblioteki (która nie została zaprojektowana do tego wyzwania), z wyjątkiem funkcji bibliotecznych do testowania pierwotności lub faktorowania liczb całkowitych. Jednak do celów punktacji uruchomię Twój kod na moim komputerze, więc podaj jasne instrukcje, jak uruchomić go na Ubuntu.

Moja maszyna Czasy zostaną uruchomione na moim komputerze. Jest to standardowa instalacja Ubuntu na 8-rdzeniowym procesorze AMD FX-8350 8 GB. Oznacza to również, że muszę być w stanie uruchomić Twój kod.

Detale

• Zabiję twój kod po 2 minutach, chyba że wcześniej zabraknie mu pamięci. Dlatego powinien upewnić się, że coś wyprowadził przed odcięciem.
• Nie możesz używać żadnego zewnętrznego źródła liczb pierwszych.
• Możesz użyć probabilistycznych metod testowania liczb pierwszych, chociaż Mego powiedział mi, że przy dobrych tabelach Miller-Rabin może testować do 341,550,071,728,321 (lub nawet wyżej) deterministycznie. Zobacz także http://miller-rabin.appspot.com/ .

Najlepsze wpisy, które sprawdzają wszystkie liczby całkowite od 1

• 756 przez cat in Go
• 756 autorstwa El'endii Starman w Pythonie
• 1932 Adnan w C # (używając mono 3.2.8)
• 2640 przez yeti w Pythonie (przy użyciu pypy 4.01)
• 2754 autorstwa Reto Koradi w C ++
• 3486 autorstwa Petera Taylora w Javie
• 3900 przez primo w RPython (przy użyciu pypy 4.01)
• 4176 autor: The Coder in Java

Najlepsze wpisy, które mogą pominąć dużą liczbę liczb całkowitych, aby znaleźć dużą lukę

• 14226 autorstwa Reto Koradi w C ++
• 22596 przez primo w RPython (używając pypy 4.01). Rekord osiągnięty po 5 sekundach!

RPython (PyPy 4.0.1), 4032

RPython jest ograniczonym podzbiorem Pythona, który można przetłumaczyć na C, a następnie skompilować za pomocą RPython Toolchain. Jego wyraźnym celem jest pomoc w tworzeniu tłumaczy językowych, ale można go również używać do kompilacji prostych programów.

Aby skompilować, pobierz bieżące źródło PyPy (PyPy 4.0.1) i uruchom następujące czynności:

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

Wynikowy plik wykonywalny zostanie nazwany good-primes-club podobny w bieżącym katalogu roboczym.

Uwagi dotyczące wdrażania

Generator liczb pierwszych primesjest nieograniczonym sito Eratostenesa, które używa koła, aby uniknąć wielokrotności 2 , 3 , 5 lub 7 . Wywołuje się także rekurencyjnie, aby wygenerować następną wartość do użycia przy oznaczaniu. Jestem bardzo zadowolony z tego generatora. Profilowanie linii ujawnia, że ​​dwie najwolniejsze linie to:

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

więc nie sądzę, że jest wiele miejsca na ulepszenia, inne niż być może przy użyciu większego koła.

W celu sprawdzenia „dobroci”, najpierw wszystkie czynniki dwóch są usuwane z n-1 , używając nieco kręcącego się hacka, aby znaleźć największą potęgę dwóch, która jest dzielnikiem (n-1 & 1-n). Ponieważ p-1 jest koniecznie nawet dla dowolnej liczby pierwszej p> 2 , wynika z tego, że 2 musi być jednym z odrębnych czynników pierwszych. To, co pozostaje, jest wysyłane do is_prime_powerfunkcji, która robi to, co sugeruje jej nazwa. Sprawdzanie, czy wartość jest siłą pierwszą, jest „prawie wolne”, ponieważ odbywa się to jednocześnie z kontrolą pierwotności, przy co najwyżej O (log _p n) operacjach, gdzie p jest najmniejszym pierwszym czynnikiem n. Podział próbny może wydawać się nieco naiwny, ale według moich testów jest to najszybsza metoda dla wartości mniejszych niż 2 ³² . Oszczędzam trochę, ponownie wykorzystując koło z sita. W szczególności:

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

poprzez iterację na kole o długości 48, p*p < nczek zostanie pominięty tysiące razy, przy niskiej, niskiej cenie nie większej niż 48 dodatkowych operacji modulo. Pomija również ponad 77% wszystkich kandydatów, a nie 50%, biorąc tylko szanse.

Ostatnie kilka wyników to:

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

Kod jest również prawidłowym Pythonem i powinien osiągnąć 3588 ~ 3900, gdy jest uruchamiany z najnowszym interpreterem 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

To przesłanie różni się nieco od innych opublikowanych do tej pory, ponieważ nie sprawdza wszystkich dobrych liczb pierwszych, ale wykonuje stosunkowo duże skoki. Wadą tego jest to, że nie można użyć sit ^{[Stoję skorygowany?]} , Więc trzeba polegać całkowicie na testach pierwotności, które w praktyce są nieco wolniejsze. Istnieje również szczęśliwy środek między stopą wzrostu a liczbą wartości sprawdzanych za każdym razem. Mniejsze wartości są znacznie szybsze do sprawdzenia, ale większe wartości mają większe luki.

Aby uspokoić bogów matematyki, postanowiłem zastosować sekwencję podobną do Fibonacciego, mając następny punkt początkowy jako sumę dwóch poprzednich. Jeśli po sprawdzeniu 10 par nie zostaną znalezione żadne nowe rekordy, skrypt przejdzie do następnego.

Ostatnie kilka wyników to:

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

Podczas kompilacji używane są 64-bitowe liczby całkowite, chociaż w kilku miejscach zakłada się, że można dodać dwie liczby całkowite bez przepełnienia, więc w praktyce można użyć tylko 63 liczb całkowitych. Po osiągnięciu 62 bitów znaczących wartość prądu jest dwukrotnie zmniejszana o połowę, aby uniknąć przepełnienia obliczeń. Rezultatem jest to, że tasuje skryptów poprzez wartości na 2 ⁶⁰ - 2 ⁶² przedziale. Nieprzekraczanie natywnej precyzji liczb całkowitych powoduje również, że skrypt jest interpretowany szybciej.

Do potwierdzenia tego wyniku można użyć następującego skryptu 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)

Prawdopodobnie 4032, mieszane sito Atkin-Bernstein i „deterministyczny” Miller-Rabin

Rozkład na koła i dobre liczby pierwsze

Jest bardzo oczywiste, że z wyjątkiem 2, 3 i 5, każda liczba pierwsza jest chroniona prawem autorskim do 2 * 3 * 5 = 60. Istnieje 16 klas równoważności modulo 60, które są coprime do 60, więc każdy test pierwotności musi tylko sprawdzić te 16 przypadków.

Kiedy jednak szukamy „dobrych” liczb pierwszych, możemy jeszcze bardziej rozrzedzić stado. Jeśli gcd(x, 60) = 1okaże się, że w większości przypadków gcd(x-1, 60)jest to 6 lub 10. Istnieje 6 wartości, xdla których jest to 2:

17, 23, 29, 47, 53, 59

Dlatego możemy wstępnie obliczyć „dobre” liczby pierwsze formy 2^a 3^b + 1i 2^a 5^b + 1scalić je w wynik generatora liczb pierwszych, który uznaje tylko 10% liczb za nawet potencjalne liczby pierwsze.

Uwagi dotyczące wdrożenia

Ponieważ miałem już implementację Java sita Atkin-Bernstein, która używa już koła jako kluczowego elementu, naturalnie wydawało się, że usuwam niepotrzebne szprychy i dostosowuję kod. Początkowo próbowałem wykorzystać architekturę producent-konsument w celu wykorzystania 8 rdzeni, ale zarządzanie pamięcią było zbyt nieuporządkowane.

Aby sprawdzić, czy liczba pierwsza jest „dobrą” liczbą pierwszą, używam „deterministycznego” testu Millera-Rabina (co tak naprawdę oznacza test Millera-Rabina, który ktoś wcześniej sprawdził na podstawie listy generowanej deterministycznie). Można to z pewnością przepisać, aby użyć również Atkina-Bernsteina, z pewnym buforowaniem, aby pokryć zakresy odpowiadające sqrt, cbrt itp., Ale nie jestem pewien, czy byłoby to ulepszenie (ponieważ testowałoby wiele liczb, które Nie muszę testować), więc to eksperyment na kolejny dzień.

Na moim dość starym komputerze to działa

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

a właściwie dokładnie 2 minuty

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

odpowiednio o 3:10, 3:20 i 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][]);
    }
}