Return to Answer

This submission is slightly different than the others posted so far, in that it doesn't check all good primes, but instead makes relatively large jumps. One disadvantage of doing this is that sieves cannot be used ^{[I stand corrected?]}, so one has to rely entirely on primality testing which in practice is quite a bit slower. There's also a happy medium to be found between the rate of growth, and the number of values checked each time. Smaller values are much faster to check, but larger values are more likely to have larger gaps.

This submission is slightly different than the others posted so far, in that it doesn't check all good primes, but instead makes relatively large jumps. One disadvantage of doing this is that sieves cannot be used ^{[I stand corrected?]}, so one has to rely entirely on primality testing which in practice is quite a bit slower. There's also a happy medium to be found between the rate of growth, and the number of values checked each time. Smaller values are much faster to check, but larger values are more likely to have larger gaps.

from time import time

# 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
      if n not in sieve:
        if n < q:
          yield n
        else:
          sieve[q + dindices[p%105]*p] = p
          p = pg.next()
          q = p*p
      else:
        stp = sieve.pop(n)
        nxt = n/stp
        nxt += dindices[nxt%105]
        while nxt*stp in sieve: nxt += dindices[nxt%105]
        sieve[nxt*stp] = stp

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):
  m = 0
  p = q = 7
  t0 = time()
  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)

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

6420 (12519586667324027 - 12519586667317607) 2.656000s
6720 (707871808582625903 - 707871808582619183) 5.015000s
8880 (626872872579606869 - 626872872579597989) 6.750000s
10146 (1206929709956703809 - 1206929709956693663) 31.812000s
22596 (918415168400717543 - 918415168400694947) 57.218000s

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:

    # modular multiplication a*b (mod m)
    # Schrage's method
    def mul_mod(a, b, m):
      b %= m
      if b:
        a %= m
        q = m/b; s = a/q
        if s:
          r = m - q*b; t = a - s*q
          if r >= q:
            # unrolled once for speed
            a1, b1 = r, s
            q1 = m/b1; s1 = a1/q1
            if s1:
              r1 = m - q1*b1; t1 = a1 - s1*q1
              if r1 >= q1:
                return (b*t - (b1*t1 - mul_mod(r1, s1, m)))%m
              return (b*t - (b1*t1 - r1*s1))%m
          return (b*t - r*s)%m
        return a*b
      return 0

  # 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:
  r_int64 = int
  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)

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

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)