28
\$\begingroup\$

Following the fine tradition of questions such as Find the largest prime whose length, sum and product is prime , this is a variant on a largest prime challenge.

Input

Your code should not take any input.

Definition

We say a prime p is good if p-1 has exactly 2 distinct prime factors.

Output

Your code should output the absolute difference between consecutive good primes q and p so that |q-p| is as large as possible and q is the smallest good prime larger than p. You may output any number of good pairs and your last output will be taken as the score.

Example

The sequence of the first 55 good primes is https://oeis.org/A067466 .

Score

Your score is simply |q-p| for the pair of good primes you output.

Languages and libraries

You can use any language or library you like (that wasn't designed for this challenge) except for any library functions for primality testing or factoring integers. However, for the purposes of scoring I will run your code on my machine so please provide clear instructions for how to run it on Ubuntu.

My Machine The timings will be run on my machine. This is a standard Ubuntu install on an 8GB AMD FX-8350 Eight-Core Processor. This also means I need to be able to run your code.

Details

  • I will kill your code after 2 minutes unless it starts to run out of memory before that. It should therefore make sure to output something before the cut off.
  • You may not use any external source of primes.
  • You may use probabilistic prime testing methods although I am told by Mego that with good tables, Miller-Rabin can test up to 341,550,071,728,321 (or even higher) deterministically. See also http://miller-rabin.appspot.com/ .

Best entries that check all integers from 1

  • 756 by cat in Go
  • 756 by El'endia Starman in Python
  • 1932 by Adnan in C# (using mono 3.2.8)
  • 2640 by yeti in Python (using pypy 4.01)
  • 2754 by Reto Koradi in C++
  • 3486 by Peter Taylor in Java
  • 3900 by primo in RPython (using pypy 4.01)
  • 4176 by The Coder in Java

Best entries that may skip a large number of integers to find a large gap

  • 14226 by Reto Koradi in C++
  • 22596 by primo in RPython (using pypy 4.01). Record reached after 5 seconds!
\$\endgroup\$
  • \$\begingroup\$ This definition looks similar to the definition of Safe prime, and apart from 5 = 2*2 +1 every safe prime is a "good prime". (Though there are good primes which are not safe primes, like 13 = 2*2*3 + 1, so I guess this doesn't help with the challenge.) \$\endgroup\$ – Paŭlo Ebermann Dec 6 '15 at 21:37
  • 1
    \$\begingroup\$ Incredibly relevant paper by djb. \$\endgroup\$ – orlp Dec 7 '15 at 1:50
  • \$\begingroup\$ @PaŭloEbermann Am I right in thinking that is not even known for sure if there are an infinite number of safe primes? Would this mean we don't know for sure there are an infinite number of good primes? \$\endgroup\$ – user9206 Dec 7 '15 at 6:30
  • \$\begingroup\$ @Lembik I'm not really an expert about safe primes, I just noticed that the definitions are quite similar and looked the safe primes up. \$\endgroup\$ – Paŭlo Ebermann Dec 7 '15 at 19:58
  • \$\begingroup\$ i did it in Labview just now which i guess you wont be able to run. Im getting to 1686 right now, is there a way for me to get i the ranking? if yes id go and optimisze it a little. \$\endgroup\$ – Eumel Dec 10 '15 at 10:53

11 Answers 11

12
+50
\$\begingroup\$

RPython (PyPy 4.0.1), 4032

RPython is a restricted subset of Python, which can be translated to C and then compiled using the RPython Toolchain. Its expressed purpose is to aid in the creation of language interpreters, but it can also be used to compile simple programs.

To compile, download the current PyPy source (PyPy 4.0.1), and run the following:

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

The resulting executable will be named good-primes-c or similar in the current working directory.


Implementation Notes

The prime number generator primes is an unbounded Sieve of Eratosthenes, which uses a wheel to avoid any multiples of 2, 3, 5, or 7. It also calls itself recursively to generate the next value to use for marking. I'm quite satisfied with this generator. Line profiling reveals that the slowest two lines are:

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

so I don't think there's much room for improvement, other than perhaps using a larger wheel.

For the "goodness" check, first all factors of two are removed from n-1, using a bit-twiddling hack to find largest power of two which is a divisor (n-1 & 1-n). Because p-1 is necessarily even for any prime p > 2, it follows that 2 must be one of the distinct prime factors. What remains is sent to the is_prime_power function, which does what its name implies. Checking if a value is a prime power is "nearly free", since it is done simultaneously with the primality check, with at most O(logpn) operations, where p is the smallest prime factor of n. Trial division might seem a bit naïve, but by my testing it is the fastest method for values less than 232. I do save a bit by reusing the wheel from the sieve. In particular:

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

by iterating over a wheel of length 48, the p*p < n check will be skipped thousands of times, at the low, low price of no more than 48 additional modulo operations. It also skips over 77% of all candidates, rather than 50% by taking only odds.

The last few outputs are:

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

The code is also valid Python, and should reach 3588 ~ 3900 when run with a recent PyPy interpreter.


# 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

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.

To appease the math gods, I've decided to follow a Fibonacci-like sequence, having the next starting point as the sum of the previous two. If no new records are found after checking 10 pairs, the script moves on the the next.

The last few outputs are:

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

When compiled, 64-bit integers are used, although it is assumed in a few places that two integers may be added without overflow, so in practice only 63 are usable. Upon reaching 62 significant bits, the current value is halved twice, to avoid overflow in the calculation. The result is that the script shuffles through values on the 260 - 262 range. Not surpassing native integer precision also makes the script faster when interpreted.

The following PARI/GP script can be used to confirm this result:

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)
\$\endgroup\$
  • \$\begingroup\$ Your welome ;) Minor update, reaches 3330 about 15 seconds faster on my machine (and soon after runs out of memory...). \$\endgroup\$ – primo Dec 8 '15 at 17:15
  • 1
    \$\begingroup\$ It does indeed. \$\endgroup\$ – user9206 Dec 8 '15 at 19:13
  • 1
    \$\begingroup\$ @Lembik I think there may be some unexplored potential there. The best I've been able to locate by placing "random depth charges" (sequences that grow like n!) is 8274 (85773786705365303 - 85773786705357029). I may add it as a bonus submission. \$\endgroup\$ – primo Dec 11 '15 at 16:56
  • 1
    \$\begingroup\$ Using pypy (not compiled) I get 13386 (32770812521685383 - 32770812521671997) 21.64s . That's pretty fast! \$\endgroup\$ – user9206 Dec 12 '15 at 19:19
  • 1
    \$\begingroup\$ 22596 (918415168400717543 - 918415168400694947) 4.786576s :) \$\endgroup\$ – user9206 Dec 16 '15 at 6:43
19
\$\begingroup\$

Probably 4032, mixed Atkin-Bernstein sieve and "deterministic" Miller-Rabin

Wheel factorisation and good primes

It's highly obvious that with the exceptions of 2, 3, and 5, every prime is coprime to 2*3*5 = 60. There are 16 equivalence classes modulo 60 which are coprime to 60, so any primality test only needs to check those 16 cases.

However, when we're looking for "good" primes we can thin the herd even more. If gcd(x, 60) = 1, we find that in most cases gcd(x-1, 60) is either 6 or 10. There are 6 values of x for which it is 2:

17, 23, 29, 47, 53, 59

Therefore we can precompute the "good" primes of the form 2^a 3^b + 1 and 2^a 5^b + 1 and merge them into the result of a prime generator which only considers 10% of numbers as even potential primes.

Implementation notes

Since I already had a Java implementation of the Atkin-Bernstein sieve lying around, and that already uses a wheel as a key component, it seemed natural to strip out the unnecessary spokes and adapt the code. I originally tried using a producer-consumer architecture to exploit the 8 cores, but memory management was too messy.

To test whether a prime is a "good" prime, I'm using a "deterministic" Miller-Rabin test (which really means a Miller-Rabin test which someone else has pre-checked against a list generated deterministically). This can certainly be rewritten to also use Atkin-Bernstein, with some caching to cover the ranges corresponding to sqrt, cbrt, etc., but I'm not sure whether it would be an improvement (because it would be testing lots of numbers which I don't need to test), so that's an experiment for another day.

On my fairly old computer this runs to

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

in pretty much exactly 2 minutes, and then

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

at 3:10, 3:20, and 3:30 respectively.

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][]);
    }
}

Save as PPCG65876.java, compile as javac PPCG65876.java, and run as java -Xmx1G PPCG65876.

\$\endgroup\$
  • \$\begingroup\$ I thought you would probably do something that is way above my head. ;) Lembik's rules exclude library functions for prime testing, though, so I think you'll have to use your own. \$\endgroup\$ – Reto Koradi Dec 8 '15 at 0:22
  • \$\begingroup\$ @RetoKoradi, yes, on re-reading I agree that methods in "You may use probabilistic prime testing methods" means techniques rather than functions. Replacing it gives a notable speedup too, so extra thanks for pointing it out. \$\endgroup\$ – Peter Taylor Dec 8 '15 at 8:09
  • \$\begingroup\$ Thanks for this! Surprisingly it only gets to 3486 on my PC. On the command line I also don't seem to need -Xmx1G. \$\endgroup\$ – user9206 Dec 8 '15 at 10:08
  • \$\begingroup\$ Do you get much higher values if you let it run longer? I just stopped mine after about 40 hours. It found 6216 as the largest difference (with prime values around 16 billion) somewhere around 12-24 hours, and nothing more after that before I stopped it. The new "high scores" definitely get rarer and rarer after some time. \$\endgroup\$ – Reto Koradi Dec 8 '15 at 16:48
  • 1
    \$\begingroup\$ @RetoKoradi, I haven't let it run for much more than 15 minutes. I'm working on approaches to speed up the isGood check. \$\endgroup\$ – Peter Taylor Dec 8 '15 at 17:22
10
\$\begingroup\$

C++, 2754 (all values, brute force primality test)

This is brute force, but it's a start before our resident mathematicians might get to work with something more efficient.

I can add some more explanation if necessary, but it's probably very obvious from the code. Since if p is a prime other than 2, we know that p - 1 is even, and one of the two factors is always 2. So we enumerate the primes, reduce p - 1 by all factors 2, and check that the remaining value is either a prime, or that all its factors are the same prime.

Code:

#include <stdint.h>
#include <vector>
#include <iostream>

int main()
{
    std::vector<uint64_t> primes;
    uint64_t prevGoodVal = 0;
    uint64_t maxDiff = 0;

    for (uint64_t val = 3; ; val += 2)
    {
        bool isPrime = true;
        std::vector<uint64_t>::const_iterator itFact = primes.begin();
        while (itFact != primes.end())
        {
            uint64_t fact = *itFact;
            if (fact * fact > val)
            {
                break;
            }

            if (!(val % fact))
            {
                isPrime = false;
                break;
            }

            ++itFact;
        }

        if (!isPrime)
        {
            continue;
        }

        primes.push_back(val);

        uint64_t rem = val;
        --rem;
        while (!(rem & 1))
        {
            rem >>= 1;
        }

        if (rem == 1)
        {
            continue;
        }

        bool isGood = true;
        itFact = primes.begin();
        while (itFact != primes.end())
        {
            uint64_t fact = *itFact;
            if (fact * fact > rem)
            {
                break;
            }

            if (!(rem % fact))
            {
                while (rem > fact)
                {
                    rem /= fact;
                    if (rem % fact)
                    {
                        break;
                    }
                }

                isGood = (rem == fact);
                break;
            }

            ++itFact;
        }

        if (isGood)
        {
            if (prevGoodVal)
            {
                uint64_t diff = val - prevGoodVal;
                if (diff > maxDiff)
                {
                    maxDiff = diff;
                    std::cout << maxDiff
                              << " (" << val << " - " << prevGoodVal << ")"
                              << std::endl;
                }
            }

            prevGoodVal = val;
        }
    }

    return 0;
}

The program prints the difference as well as the corresponding two good primes any time a new maximum difference is found. Sample output from the test run on my machine, where the reported value of 2754 is found after about 1:20 minutes:

4 (11 - 7)
6 (19 - 13)
8 (37 - 29)
14 (73 - 59)
24 (137 - 113)
30 (227 - 197)
32 (433 - 401)
48 (557 - 509)
50 (769 - 719)
54 (1283 - 1229)
60 (1697 - 1637)
90 (1823 - 1733)
108 (2417 - 2309)
120 (3329 - 3209)
126 (4673 - 4547)
132 (5639 - 5507)
186 (7433 - 7247)
222 (8369 - 8147)
258 (16487 - 16229)
270 (32507 - 32237)
294 (34157 - 33863)
306 (35879 - 35573)
324 (59393 - 59069)
546 (60293 - 59747)
570 (145823 - 145253)
588 (181157 - 180569)
756 (222059 - 221303)
780 (282617 - 281837)
930 (509513 - 508583)
1044 (1046807 - 1045763)
1050 (1713599 - 1712549)
1080 (1949639 - 1948559)
1140 (2338823 - 2337683)
1596 (3800999 - 3799403)
1686 (6249743 - 6248057)
1932 (12464909 - 12462977)
2040 (30291749 - 30289709)
2160 (31641773 - 31639613)
2190 (34808447 - 34806257)
2610 (78199097 - 78196487)
2640 (105072497 - 105069857)
2754 (114949007 - 114946253)
^C

real    1m20.233s
user    1m20.153s
sys 0m0.048s
\$\endgroup\$
7
\$\begingroup\$

C++, 14226 (high values only, Miller-Rabin test)

Posting this separately because it is entirely different from my initial solution, and I did not want to completely replace a post that had gotten a number of upvotes.

Thanks to @primo for pointing out a problem with the original version. There was an overflow for large numbers in the prime number test.

This takes advantage of some insights that have been gained during the evolution of other solutions. The main observations are:

  • Since the results clearly show that the gaps get larger as the primes themselves get larger, there is no point in bothering with small primes. Exploring large prime values is much more effective.
  • Probabilistic prime testing is required for primes of this size.

Based on this, the method employed here is fairly simple:

  • For primality testing, the Miller-Rabin test is used. The implementation is based on the pseudo code on the wikpedia page. With the bases used, it will deliver correct values up to 3825123056546413051 (see OEIS A014233), which is plenty for the range of values used here.
  • To determine if primes are good primes, the powers of 2 are stripped. Factorizing the remaining value would be very expensive, but is unnecessary. Instead, I calculate the much fewer possible roots using double math, and see if any of them produces an integer that is in fact the correct root.
  • Math is mostly using 64-bit unsigned values, with 128-bit unsigned values needed for some temporary values in the primality test.
  • Since I use double math for the roots, and a double can exactly represent integers of at most 53 bits, the maximum size that can safely be handled by this code is 54 bits (the number converted to double is at most half the size of the prime).
  • Since 54 bits was the maximum size for the number I was confident using, I start with a number that is somewhat smaller than the maximum 54-bit number. The code does report larger gaps for even larger start values, and they are probably correct, but I can't be as sure.

Results:

1266 (16888498602640739 - 16888498602639473)
1470 (16888498602645563 - 16888498602644093)
2772 (16888498602651629 - 16888498602648857)
2862 (16888498602655829 - 16888498602652967)
3120 (16888498602675053 - 16888498602671933)
3756 (16888498602685769 - 16888498602682013)
4374 (16888498602696257 - 16888498602691883)
5220 (16888498602745493 - 16888498602740273)
5382 (16888498603424039 - 16888498603418657)
5592 (16888498603511279 - 16888498603505687)
5940 (16888498603720697 - 16888498603714757)
6204 (16888498605020837 - 16888498605014633)
6594 (16888498605999017 - 16888498605992423)
14226 (16888498608108539 - 16888498608094313)
^C

real    0m26.335s
user    0m26.312s
sys 0m0.008s

Code:

#include <stdint.h>
#include <cmath>
#include <iostream>

uint64_t intRoot(uint64_t a, int p)
{
    double e = 1.0 / static_cast<double>(p);
    double dRoot = pow(a, e);

    return static_cast<uint64_t>(dRoot + 0.5);
}

uint64_t intPow(uint64_t a, int e)
{
    uint64_t r = 1;

    while (e)
    {
        if (e & 1)
        {
            r *= a;
        }

        e >>= 1;
        a *= a;
    }

    return r;
}

uint64_t modPow(uint64_t a, uint64_t e, uint64_t m)
{
    uint64_t r = 1;
    a %= m;

    while (e)
    {
        if (e & 1)
        {
            __uint128_t t = r;
            t *= a;
            t %= m;
            r = t;
        }

        e >>= 1;
        __uint128_t t = a;
        t *= a;
        t %= m;
        a = t;
    }

    return r;
}

bool isPrime(uint64_t n)
{
    const uint64_t a[] = {2, 3, 5, 7, 11, 13, 17, 19, 23};

    if (n < 2)
    {
        return false;
    }

    for (int k = 0; k < 9; ++k)
    {
        if (n == a[k])
        {
            return true;
        }

        if (n % a[k] == 0)
        {
            return false;
        }
    }

    int r = __builtin_ctzll(n - 1);
    uint64_t d = (n - 1) >> r;

    for (int k = 0; k < 9; ++k)
    {
        uint64_t x = modPow(a[k], d, n);

        if (x == 1 || x == n - 1)
        {
            continue;
        }

        bool comp = true;
        for (int i = 0; i < r - 1; ++i)
        {
            x = modPow(x, 2, n);
            if (x == 1)
            {
                return false;
            }
            if (x == n - 1)
            {
                comp = false;
                break;
            }
        }

        if (comp)
        {
            return false;
        }
    }

    return true;
}

int main()
{
    uint64_t prevGoodVal = 0;
    uint64_t maxDiff = 0;

    for (uint64_t val = (1ull << 54) - (1ull << 50) + 1; ; val += 2)
    {
        if (isPrime(val))
        {
            uint64_t d = static_cast<double>((val - 1) >> __builtin_ctzll(val - 1));
            bool isGood = false;

            if (isPrime(d))
            {
                isGood = true;
            }
            else
            {
                for (int e = 2; ; ++e)
                {
                    uint64_t r = intRoot(d, e);
                    if (r < 3)
                    {
                        break;
                    }

                    if (intPow(r, e) == d && isPrime(r))
                    {
                        isGood = true;
                        break;
                    }
                }
            }

            if (isGood)
            {
                if (prevGoodVal)
                {
                    uint64_t diff = val - prevGoodVal;
                    if (diff > maxDiff)
                    {
                        maxDiff = diff;
                        std::cout << maxDiff
                                  << " (" << val << " - " << prevGoodVal << ")"
                                  << std::endl;
                    }
                }

                prevGoodVal = val;
            }
        }
    }

    return 0;
}
\$\endgroup\$
  • \$\begingroup\$ @primo Should be correct now. There was an overflow where I multiplied two 64-bit numbers in the primality test, causing it to report "composite" for some large prime numbers. Thanks for pointing it out. Let me know if you still see a problem. \$\endgroup\$ – Reto Koradi Dec 15 '15 at 8:35
  • \$\begingroup\$ That's a good one. The race is on? ;) \$\endgroup\$ – primo Dec 15 '15 at 8:36
  • \$\begingroup\$ @primo I had some considerably larger values, but they would use primes that cannot be completely represented by a double. I think it would still give a precise enough approximation of the root to produce correct results. Or I could implement a root finding algorithm that does not use doubles. But I won't be able to spend more time on this before the bounty expires... \$\endgroup\$ – Reto Koradi Dec 15 '15 at 8:43
  • \$\begingroup\$ Your answer also reaches its maximum in 4 seconds! (Just like primo's.) \$\endgroup\$ – user9206 Dec 15 '15 at 16:35
6
\$\begingroup\$

PyPy-2.4.0

Python-2

The x files ...

Episode: "Look mom! Not a single division!"

;-)

M = g = 0
B = L = {}
n = 2
while 1:
        if n in L:
                B = P = L[n]
                del L[n]
        else:
                if len(B) == 2:
                        if g:
                                m = n - g
                                if M < m:
                                        M = m
                                        print n, g, m
                        g = n
                P = [n]
        for p in P:
                npp = n + p
                if npp in L:
                        if p not in L[npp]:
                                L[npp] += [p]
                else:
                        L[npp] = [p]
        n += 1

I tested it on Debian8 with PyPy-2.4.0 and Python2 started like:

timeout 2m pypy -O x
timeout 2m python2 -O x

If there is really plenty of RAM, the del L[n] line may be deleted.


The basic prime number generator is this:

L = {}
n = 2

while 1:

        if n in L:
                P = L[n]
                del L[n]
        else:
                print n
                P = [n]

        for p in P:
                npp = n + p
                if npp in L:
                        if p not in L[npp]:
                                L[npp] += [p]
                else:
                        L[npp] = [p]

        n += 1

It basically does exactly what the sieve of Eratosthenes does but in a different order.

L is dictionary but can be seen as list (tape) of lists of numbers. Nonexistent cells L[n] are interpreted as n has no known prime divisiors up to now.

The while loop does a prime or not prime decission on each turn for L[n].

  • If L[n] exists (same as n in L), P = L[n] is a list of distinct prime divisors of n. So n is not a prime.

  • If L[n] does not exist, no prime divisor was found. So n must be prime then with P = [n] being the known divisor.

Now P is the list of known prime divisors for both cases.

The for p in P loop moves every entry of P forward by the distance of it's value on the tape of numbers.

This is how the divisors jump on the tape and this is the reason why these jumping numbers have to be prime. New numbers only get on the tape by the else decission above and those are numbers with no known divisors other than them self. Nonprimes never get into these lists L[n].

The primes jumping on the list are all distinct because every number n is looked at only once and only is added as divisor (if not prime:) 0 or (if prime:) 1 times. Known prime divisors only will move forward but never get duplicated. So L[n] always will hold distinct prime divisors or will be empty.


Back to the upper program for the good primes gaps:

    if n in L:
            B = P = L[n]

...keeps the prime divisors of n in B when n is known not to be prime.

If n is recognised to be prime, B holds the list of prime divisors of the previous loop pass looking at n-1:

    else:
            if len(B) == 2:

...so len(B) == 2 means n - 1 has two distinct prime divisors.

                        if g:
                                m = n - g
                                if M < m:
                                        M = m
                                        print n, g, m
                        g = n

g just remembers the last seen good prime before the new one, M is the length of the previous maximum good prime gap and m the length of the newly found gap.


Happy end.

\$\endgroup\$
  • \$\begingroup\$ Nice solution. For me, this hits 2640 in about 117s. \$\endgroup\$ – primo Dec 10 '15 at 6:51
  • 1
    \$\begingroup\$ Could you add a little explanation please. \$\endgroup\$ – user9206 Dec 10 '15 at 8:04
  • 1
    \$\begingroup\$ @Lembik: Done... \$\endgroup\$ – user19214 Dec 10 '15 at 14:09
4
\$\begingroup\$

C#, probably 1932

I did find out that, the faster your algorithm is for finding primes, the better your score. I'm also pretty sure that my algorithm is not the most optimal method for prime searching.

using System;
using System.Collections.Generic;

namespace GoodPrimes
{
    class Program
    {
        static void Main(string[] args)
        {
            int[] list_of_primes = new int[168]{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, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997};
            bool is_last_prime = false;
            int last_prime = 0;
            int max_value = 0;
            int old_max_value = 1000000;
            int old_min_value = 3;
            HashSet<int> primeSet = new HashSet<int>();
            primeSet.Add(2);
            int X = 0;
            Console.WriteLine("Initialize primes until " + old_max_value);
            for (int i = old_min_value; i < old_max_value; i++)
            {
                if (IsPrime(i, primeSet))
                    primeSet.Add(i);
            }
            old_min_value = old_max_value;
            for (int i = 3; ; i += 2)
            {
                if (i > old_max_value)
                {
                    old_max_value += 500000;
                    Console.WriteLine("Initialize primes until " + old_max_value);
                    for (int j = old_min_value; j < old_max_value; j++)
                    {
                        for(int k = 0; k < list_of_primes.Length; k++)
                            if(j % list_of_primes[k] == 0 && j > list_of_primes[k])
                                continue;
                        if (IsPrime(j, primeSet))
                            primeSet.Add(j);
                    }
                    old_min_value = old_max_value;
                }
                if (primeSet.Contains(i))
                {
                    is_last_prime = false;
                    X = (i - 1) / 2;
                    while (X % 2 == 0)
                        X = X / 2;
                    if (IsPrime(X, primeSet))
                        is_last_prime = true;
                    for (int j = 3; j < i; j++)
                    {
                        if (j % 2 == 0 && j > 2)
                            continue;
                        if (j % 3 == 0 && j > 3)
                            continue;
                        if (j % 5 == 0 && j > 5)
                            continue;
                        if (j % 7 == 0 && j > 7)
                            continue;
                        if (j % 11 == 0 && j > 11)
                            continue;
                        if (j % 13 == 0 && j > 13)
                            continue;
                        if (j % 17 == 0 && j > 17)
                            continue;

                        if (X % j == 0 || is_last_prime)
                        {
                            while (X % j == 0)
                                X = X / j;
                            if ((primeSet.Contains(j) && X == 1) || is_last_prime)
                            {
                                while (X % j == 0)
                                    X = X / j;
                                if (X == 1 || is_last_prime)
                                {
                                    if (i - last_prime > max_value)
                                    {
                                        max_value = i - last_prime;
                                        Console.WriteLine("New max value: " + max_value.ToString() + " (" + i.ToString() + "-" + last_prime.ToString() + ")");
                                    }
                                    last_prime = i;
                                }
                            }
                            break;
                        }
                    }
                }
            }
            Console.ReadLine();
        }

        private static bool IsPrime(int i, HashSet<int> j)
        {
            if (i == 2)
                return true;
            for (int m = 2; m < Math.Sqrt(System.Convert.ToDouble(i)) + 1; m++)
            {
                if (j.Contains(m))
                {
                    if (m % 2 == 0 && m > 2)
                        continue;
                    if (m % 3 == 0 && m > 3)
                        continue;
                    if (m % 5 == 0 && m > 5)
                        continue;
                    if (m % 7 == 0 && m > 7)
                        continue;
                    if (m % 11 == 0 && m > 11)
                        continue;
                    if (m % 13 == 0 && m > 13)
                        continue;
                    if (m % 17 == 0 && m > 17)
                        continue;
                    if (i % m == 0)
                        return false;
                }
            }
            return true;
        }
    }
}
\$\endgroup\$
4
\$\begingroup\$

Python 3, 546

...in two minutes on my machine, which I think is significantly less powerful than yours.

def getPrimes_parallelized(): #uses sieve of Sundaram
        yield 2
        yield 3
        P = [[4,1]]
        i = 2
        while 1:
            if P[0][0] <= i:
                while P[0][0] <= i:
                    P[0][0] += 2*P[0][1]+1
                    P.sort()
            elif P[0][0] > i:
                yield 2*i+1
                P.append([2*(i+i*i), i])
                P.sort()
            i += 1

def goodPrimes(x):
    P = getPrimes_parallelized()
    primes = []

    for p in P:
        primes.append(p)
        n = p-1
        factors = []

        for p2 in primes:
            if n%p2 == 0:
                factors.append(p2)
                while n%p2 == 0: n //= p2

            if len(factors) > x: break

        if len(factors) <= x: yield p

maxdiff = 0
GP = goodPrimes(2)
p1 = next(GP)
gp = next(GP)
gps = [(p1,gp)]

while 1:
    if gp-p1 > maxdiff:
        maxdiff = gp-p1
        print("p: %d, q: %d, |q-p|: %d" % (p1,gp,gp-p1))
    p1,gp = gp,next(GP)

I probably could make this more efficient by optimizing for the x=2 case, but eh. Good enough. :P

\$\endgroup\$
  • \$\begingroup\$ Your code just outputs p: 2, q: 3, |q-p|: 1 for me. \$\endgroup\$ – user9206 Dec 7 '15 at 17:18
  • 1
    \$\begingroup\$ @Lembik: Ah, whoops. I pared this down from the version that had plotting stuff, and left out a crucial line. Fixed. \$\endgroup\$ – El'endia Starman Dec 7 '15 at 20:46
4
\$\begingroup\$

Go, probably 756

For shame! I'm such a novice that I only naively reused some old code and expected it to work and be fast! If I reimplemented this and actually built it around good primes, it would be so much faster, but alas, I am learning. (I'll probably answer again tomorrow with a fully rebuilt solution that is purpose-built.)

package main

import "fmt"

func mkPrime(ch chan<- int) {
    for i := 2; ; i++ {
        ch <- i // Send 'i' to channel 'ch'.
    }
}

// Copy the values from channel 'in' to channel 'out',
// removing those divisible by 'prime'.
func filterPrm(in <-chan int, out chan<- int, prime int) {
    for {
        i := <-in // Receive value from 'in'.
        if i%prime != 0 {
            out <- i // Send 'i' to 'out'.
        }
    }
}

func mkPFac(max int, ch chan<- int) {
    ch <- 2
    for i := 3; i <= max; i += 2 {
        ch <- i
    }
    ch <- -1 // signal that the limit is reached
}

// Copy the values from channel 'in' to channel 'out',
// removing those divisible by 'prime'.
func filterPFac(in <-chan int, out chan<- int, prime int) {
    for i := <-in; i != -1; i = <-in {
        if i%prime != 0 {
            out <- i
        }
    }
    out <- -1
}

func calcPFactors(numToFac int) []int {
    rv := []int{}
    ch := make(chan int)
    go mkPFac(numToFac, ch)
    for prime := <-ch; (prime != -1) && (numToFac > 1); prime = <-ch {
        for numToFac%prime == 0 {
            numToFac = numToFac / prime
            rv = append(rv, prime)
        }
        ch1 := make(chan int)
        go filterPFac(ch, ch1, prime)
        ch = ch1
    }
    return rv
}

func rmDup(list []int) []int {
    var nlist []int
    for _, e := range list {
        if !isIn(e, nlist) {
            nlist = append(nlist, e)
        }
    }
    return nlist
}

func isIn(a int, list []int) bool {
    for _, b := range list {
        if b == a {
            return true
        }
    }
    return false
}

// The prime sieve: Daisy-chain Filter processes.
func main() {
    var diff, prev, high int
    ch := make(chan int) // Create a new channel.
    go mkPrime(ch)       // Launch Generate goroutine.
    for i := 0; i < 10000000000; i++ {
        prime := <-ch
        list := rmDup(calcPFactors(prime - 1))
        if len(list) == 2 {
            //fmt.Println(list, prime)
            diff = prime - prev
            //fmt.Println(diff)
            prev = prime
            if diff > high {
                high = diff
                fmt.Println(high)
            }
        }
        ch1 := make(chan int)
        go filterPrm(ch, ch1, prime)
        ch = ch1
    }
}

Uses concurrency, obviously.

\$\endgroup\$
  • 1
    \$\begingroup\$ Go is always welcome :) \$\endgroup\$ – user9206 Dec 16 '15 at 6:42
4
\$\begingroup\$

Java, 4224 (99.29 s)

Heavily Customized Sieve of Eratosthenes with advantage of BitSet

import java.util.BitSet;

public class LargeGoodPrimeGap {

    // Use this to find upto Large Gap of 4032 - Max 4032 found in 55.17 s
    // static int    limit         = 125_00_00_000;

    // Use this to find upto Large Gap of 4224 - Max 4224 found in 99.29 s
    static int    limit         = Integer.MAX_VALUE - 1;

    // BitSet is highly efficient against boolean[] when Billion numbers were involved
    // BitSet uses only 1 bit for each number
    // boolean[] uses 8 bits aka 1 byte for each number which will produce memory issues for large numbers
    static BitSet primes        = new BitSet(limit + 1);
    static int    limitSqrt     = (int) Math.ceil(Math.sqrt(limit));

    static int    maxAllowLimit = Integer.MAX_VALUE - 1;

    static long   start         = System.nanoTime();

    public static void main(String[] args) {

        genPrimes();

        findGoodPrimesLargeGap();

    }

    // Generate Primes by Sieve of Eratosthenes
    // Sieve of Eratosthenes is much efficient than Sieve of Atkins as
    // Sieve of Atkins involes Division, Modulus, Multiplication, Subtraction, Addition but
    // Sieve of Eratosthenes involves only addition and multiplication
    static void genPrimes() {

        // Check if the Given limit exceeds the Permitted Limit 2147483646 (Integer.MAX_VALUE - 1)
        // If the limit exceeded, Out the Error Message and Exit the Program
        if ( limit > maxAllowLimit ) {
            System.err.printf(String.format("Limit %d should not be Greater than Max Limit %d", limit, maxAllowLimit));
            System.exit(0);
        }

        // Mark numbers from 2 to limit + 1 as Prime
        primes.set(2, limit + 1);

        // Now all Values in primes will be true except 0 and 1,
        // True  represents     prime number 
        // False represents not prime number

        // Set the First Prime number
        int prime = 2;
        // Set the First multiple of prime
        int multiple = prime;
        // Reduce the limit by 1 if limit == Interger.MAX_VALUE - 1 to prevent
        // Integer overflow on multiple variable
        int evenLimit = limit == Integer.MAX_VALUE - 1 ? limit - 1 : limit;

        // Mark all Even Numbers as Not Prime except 2
        while ( (multiple += prime) <= evenLimit ) {
            primes.clear(multiple);
        }

        // If evenLimit != limit, set last even number as Not Prime
        if ( evenLimit != limit ) {
            primes.clear(limit);
        }

        int primeAdd;

        // Set odd multiples of each Prime as not Prime;
        // prime <= limitSqrt -> Check Current Prime <= SQRT(limit)
        // prime = primes.nextSetBit(prime + 1) -> Assign the next True (aka Prime) value as Current Prime
        //  ^ - Above initialisation is highly efficient as Next True check is only based on bits
        // prime > 0 -> To handle -ve values returned by above True check if no more True is to be found
        for ( prime = 3; prime > 0 && prime <= limitSqrt; prime = primes.nextSetBit(prime + 1) ) {
            // All Prime Numbers except 2 were odd numbers
            // Adding a Prime number with itself will result in an Even number,
            // but all the Even numbers were already marked as not Prime.
            // So every odd multiple (3rd, 5th, 7th, ...) of Current Prime will only be marked as not Prime
            // and skipping all the even multiples (2nd, 4th, 6th, ...)
            // This reduces the time for prime calculation by ~50% when comparing with all multiples marking
            primeAdd = prime + prime;
            // multiple = prime * prime -> Unmarked Prime will appear only from this number as previous values
            // are already marked as Non Prime by previous prime multiples
            // multiple += primeAdd -> Increases the multiple by multiple + (CurrentPrime x 2) which will
            // always be a odd multiple (5th, 7th, 9th, ...)
            for ( multiple = prime * prime; multiple <= limit && multiple > 0; multiple += primeAdd ) {
                // Clear or False the multiple if it True
                primes.clear(multiple);
            }
        }

        double end = (System.nanoTime() - start) / 1000000000.0;
        System.out.printf("Total Primes upto %d = %d in %.2f s", limit, primes.cardinality(), end);

    }

    static void findGoodPrimesLargeGap() {

        int prevGP = 7;
        int prevDiff = 0;

        for ( int i = 11; i <= limit && i > 0; i = primes.nextSetBit(i + 1) ) {
            int gp = i - 1;
            int distPrimes = 0;
            for ( int j = 2; j <= limitSqrt && distPrimes < 3 && j > 0; j = primes.nextSetBit(j + 1) ) {
                if ( gp % j == 0 ) {
                    ++distPrimes;
                    while ( gp % j == 0 ) {
                        gp = gp / j;
                    }
                    if ( gp <= 1 ) {
                        break;
                    }
                }
                if ( primes.get(gp) ) {
                    ++distPrimes;
                    break;
                }
            }
            if ( distPrimes == 2 ) {
                int currDiff = i - prevGP;
                if ( currDiff > prevDiff ) {
                    System.out.println(
                            String.format("(%d - %d) %d (%.2f s)", i, prevGP, prevDiff = currDiff, (System.nanoTime() - start) / 1000000000.0));
                }
                prevGP = i;
            }
        }

    }

}

Time taken is dependent on the Max Limit of the Prime numbers which are going be calculated.

For

static int    limit         = Integer.MAX_VALUE - 1;

Total Primes upto 2147483646 = 105097564 in 17.65 s
(11 - 7) 4 (17.71 s)
(19 - 13) 6 (17.71 s)
(37 - 29) 8 (17.71 s)
(73 - 59) 14 (17.71 s)
(137 - 113) 24 (17.71 s)
(227 - 197) 30 (17.71 s)
(433 - 401) 32 (17.71 s)
(557 - 509) 48 (17.71 s)
(769 - 719) 50 (17.71 s)
(1283 - 1229) 54 (17.71 s)
(1697 - 1637) 60 (17.71 s)
(1823 - 1733) 90 (17.71 s)
(2417 - 2309) 108 (17.71 s)
(3329 - 3209) 120 (17.71 s)
(4673 - 4547) 126 (17.71 s)
(5639 - 5507) 132 (17.71 s)
(7433 - 7247) 186 (17.71 s)
(8369 - 8147) 222 (17.71 s)
(16487 - 16229) 258 (17.71 s)
(32507 - 32237) 270 (17.72 s)
(34157 - 33863) 294 (17.72 s)
(35879 - 35573) 306 (17.72 s)
(59393 - 59069) 324 (17.72 s)
(60293 - 59747) 546 (17.72 s)
(145823 - 145253) 570 (17.73 s)
(181157 - 180569) 588 (17.73 s)
(222059 - 221303) 756 (17.73 s)
(282617 - 281837) 780 (17.73 s)
(509513 - 508583) 930 (17.74 s)
(1046807 - 1045763) 1044 (17.75 s)
(1713599 - 1712549) 1050 (17.77 s)
(1949639 - 1948559) 1080 (17.77 s)
(2338823 - 2337683) 1140 (17.78 s)
(3800999 - 3799403) 1596 (17.80 s)
(6249743 - 6248057) 1686 (17.85 s)
(12464909 - 12462977) 1932 (17.96 s)
(30291749 - 30289709) 2040 (18.31 s)
(31641773 - 31639613) 2160 (18.34 s)
(34808447 - 34806257) 2190 (18.41 s)
(78199097 - 78196487) 2610 (19.40 s)
(105072497 - 105069857) 2640 (20.07 s)
(114949007 - 114946253) 2754 (20.32 s)
(246225989 - 246223127) 2862 (24.01 s)
(255910223 - 255907313) 2910 (24.31 s)
(371348513 - 371345567) 2946 (27.97 s)
(447523757 - 447520673) 3084 (30.50 s)
(466558553 - 466555373) 3180 (31.15 s)
(575713847 - 575710649) 3198 (35.00 s)
(606802529 - 606799289) 3240 (36.13 s)
(784554983 - 784551653) 3330 (42.89 s)
(873632213 - 873628727) 3486 (46.39 s)
(987417437 - 987413849) 3588 (50.97 s)
(1123404923 - 1123401023) 3900 (56.60 s)
(1196634239 - 1196630297) 3942 (59.70 s)
(1247118179 - 1247114147) 4032 (61.88 s)
(1964330609 - 1964326433) 4176 (94.89 s)
(2055062753 - 2055058529) 4224 (99.29 s)
\$\endgroup\$
  • \$\begingroup\$ This is surprisingly faster than the other Java submission! \$\endgroup\$ – user9206 Dec 17 '15 at 20:07
  • \$\begingroup\$ @Lembik, I'll add more detailed explanation later today.. \$\endgroup\$ – The Coder Dec 18 '15 at 2:50
  • \$\begingroup\$ @Lembik , Highly Customized the Sieve Logic. Now the time taken to generate all primes is reduced by ~50%. So within 100s, max large diff within Integer.MAX_VALUE can be found \$\endgroup\$ – The Coder Dec 19 '15 at 18:00
3
\$\begingroup\$

Python 3, 1464

With help from Lembik, whose idea was to just check for the first two good primes after a power of two and when found immediately move on to the next power of two. If someone can use this as a jumping point, feel free. A portion of my results are below after running this in IDLE. The code follows.

Credit to primo as I grabbed their list of small primes for this code.

Edit: I have edited the code to fit the actual specifications of the problem (two distinct prime divisors not exactly two distinct prime divisors), and I implemented not moving on to the next power of two until the good primes the program has found have a gap larger than that of the last two good primes it found. I should also give credit to Peter Taylor, as I used his idea that good primes could only be a few values mod 60.

Again, I have run this on a slow computer in IDLE, so the results may be faster in something like PyPy, but I have not been able to check.

A sample of my results (p, q, q-p, time):

8392997 8393999 1002 2.6750288009643555
16814663 16815713 1050 7.312098026275635
33560573 33561653 1080 8.546097755432129
67118027 67119323 1296 10.886202335357666
134245373 134246753 1380 20.37420392036438
268522349 268523813 1464 59.23987054824829
536929187 536931047 1860 95.36681914329529

My code:

from time import time

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,223,227,229,233,239]

def good(n=0):
    end = n or 100
    time0 = time()
    x,y = 0,0
    recent_max = 0
    for i in range(2,end):
        two = 2**i
        for j in range(two+3,2*two,2):
            m=j%60
            if not(m==17or m==23or m==29or m==47or m==53or m==59): continue
            comp = 0
            for p in small_primes:
                if j % p == 0 and j != p:
                    comp = 1
                    break
            for p in range(241,int(pow(j,.5))+1,2):
                if j % p == 0 and j != p:
                    comp = 1
                    break
            if comp: continue
            d = j-1 & 1-j
            if is_prime_power((j-1)/d):
                x,y = y,j
                if x and y and y-x > recent_max:
                    print(x,y,y-x,time()-time0)
                    recent_max = y-x
                    x,y=0,0
                    break

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
    for p in range(241,int(pow(n,.5))+1,2):
        if n%p == 0:
            n //= p
            while n % p == 0: n //= p
            return n == 1
    return n > 1

good()
\$\endgroup\$
  • \$\begingroup\$ I don't think your code is correct. Do you have any justification for incrementing j by 4 rather than 2? And you seem to reject unconditionally if j-1 isn't a prime times a power of two, where you should be testing whether it's a prime power times a power of two. \$\endgroup\$ – Peter Taylor Dec 11 '15 at 17:50
  • \$\begingroup\$ @PeterTaylor Oh sweet Jesus thank you. I knew I was missing something. Two distinct prime factors, not exactly two distinct prime factors. I'll correct this in the morning. \$\endgroup\$ – Sherlock9 Dec 11 '15 at 17:58
  • \$\begingroup\$ Indeed. The next good prime after 549755815199 is 549755816417 (2^5 x 17179869263), a gap of only 1218. \$\endgroup\$ – primo Dec 11 '15 at 18:00
2
\$\begingroup\$

Go: All Integers: 5112

max |q-p| 5112 q 4278566879 p 4278561767

good.go:

// Find the largest gap between good primes
// https://codegolf.stackexchange.com/questions/65876/
//
// We say a prime p is good if p-1 has exactly 2 distinct prime factors.
//
// Your code should output the absolute difference between consecutive
// good primes q and p so that |q-p| is as large as possible and
// q is the smallest good prime larger than p. You may output any number of
// good pairs and your last output will be taken as the score.
//
// The timings will be run on a standard Ubuntu install on
// an 8GB AMD FX-8350 eight-core processor.
// http://products.amd.com/en-us/search/CPU/AMD-FX-Series/AMD-FX-8-Core-Black-Edition/FX-8350/92
//
// I will kill your code after 2 minutes unless it starts to
// run out of memory before that. It should therefore make sure to
// output something before the cut off.
//
// A067466 Primes p such there are 2 distinct prime factors in p-1.
// https://oeis.org/A067466
//
// 7, 11, 13, 19, 23, 29, 37, 41, 47, 53, 59, 73, 83, 89, 97, 101, 107, ...
//
// peterSO: max |q-p| 5112 q 4278566879 p 4278561767
// https://codegolf.stackexchange.com/a/73770/51537
//
// p is a good prime number, if
//
//   p-1 = x**a * y**b
//
// Where p is a prime number, x and y are are distinct prime numbers,
// and a and b are positive integers.
//
// For p > 2, p is odd and (p-1) is even. Therefore, either x or y = 2.

package main

import (
    "fmt"
    "math"
    "runtime"
    "time"
)

var start = time.Now()

const (
    primality = 0x80
    prime     = 0x00
    notprime  = 0x80
    distinct  = 0x7F
)

func oddPrimes(n uint64) (sieve []uint8) {
    // odd prime numbers
    sieve = make([]uint8, (n+1)/2)
    sieve[0] = notprime
    p := uint64(3)
    for i := p * p; i <= n; i = p * p {
        for j := i; j <= n; j += 2 * p {
            sieve[j/2] = notprime
        }
        for p += 2; sieve[p/2] == notprime; p += 2 {
        }
    }
    return sieve
}

func maxGoodGap(n uint64) {
    // odd prime numbers
    sieve := oddPrimes(n)
    // good prime numbers
    fmt.Println("|q-p|", " = ", "q", "-", "p", ":", "t")
    m := ((n + 1) + 1) / 2
    var max, px, qx uint64
    for i, s := range sieve {
        if s == prime {
            p := 2*uint64(i) + 1
            if p < m {
                // distinct odd prime number factors
                for j := p + 2*p; j <= m; j += 2 * p {
                    sieve[j/2]++
                }
            }
            // Remove factors of 2 from p-1.
            p1 := p - 1
            for ; p1&1 == 0; p1 >>= 1 {
            }
            // Does p-1 have exactly 2 distinct prime factors?
            // That is, one distinct prime factor other than 2.
            if sieve[p1/2]&distinct <= 1 {
                // maximum consecutive good prime gap
                px, qx = qx, p
                if max < qx-px {
                    max = qx - px
                    if px != 0 {
                        fmt.Println(max, " = ", qx, "-", px, " : ", time.Since(start))
                    }
                }
            }
        }
    }
}

func init() {
    runtime.GOMAXPROCS(1)
}

func main() {
    // Two minutes: max |q-p| 5112 q 4278566879 p 4278561767
    var n uint64 = math.MaxUint32 // 4294967295
    fmt.Println("n =", n)
    maxGoodGap(n)
    fmt.Println("n =", n, "real =", time.Since(start))
}

Output:

$ go build good.go && ./good
n = 4294967295
|q-p|  =  q - p : t
4  =  11 - 7  :  18.997478838s
6  =  29 - 23  :  19.425839298s
8  =  37 - 29  :  19.5924487s
14  =  73 - 59  :  20.351329953s
24  =  137 - 113  :  21.339752269s
30  =  227 - 197  :  22.310449147s
32  =  433 - 401  :  23.511560468s
48  =  557 - 509  :  23.904677275s
50  =  769 - 719  :  24.518310365s
54  =  1283 - 1229  :  25.350700584s
60  =  1697 - 1637  :  25.782520338s
90  =  1823 - 1733  :  25.883049102s
108  =  2417 - 2309  :  26.300049556s
120  =  3329 - 3209  :  26.735575056s
126  =  4673 - 4547  :  27.190597227s
132  =  5639 - 5507  :  27.420936586s
186  =  7433 - 7247  :  27.761805597s
222  =  8369 - 8147  :  27.909656781s
258  =  16487 - 16229  :  28.710626512s
270  =  32507 - 32237  :  29.469193619s
294  =  34157 - 33863  :  29.525197303s
306  =  35879 - 35573  :  29.578355515s
324  =  59393 - 59069  :  30.11620771s
546  =  60293 - 59747  :  30.131928104s
570  =  145823 - 145253  :  31.014864294s
588  =  181157 - 180569  :  31.223246627s
756  =  222059 - 221303  :  31.415507367s
780  =  282617 - 281837  :  31.640006297s
930  =  509513 - 508583  :  32.169485481s
1044  =  1046807 - 1045763  :  32.783669616s
1050  =  1713599 - 1712549  :  33.186784964s
1080  =  1949639 - 1948559  :  33.290533456s
1140  =  2338823 - 2337683  :  33.434568615s
1596  =  3800999 - 3799403  :  33.810580195s
1686  =  6249743 - 6248057  :  34.183678793s
1932  =  12464909 - 12462977  :  34.683651976s
2040  =  30291749 - 30289709  :  35.296022077s
2160  =  31641773 - 31639613  :  35.325773748s
2190  =  34808447 - 34806257  :  35.390646164s
2610  =  78199097 - 78196487  :  35.878632519s
2640  =  105072497 - 105069857  :  36.018381898s
2754  =  114949007 - 114946253  :  36.058571726s
2862  =  246225989 - 246223127  :  36.337844257s
2910  =  255910223 - 255907313  :  36.351442541s
2946  =  371348513 - 371345567  :  36.504506082s
3084  =  447523757 - 447520673  :  36.60250012s
3180  =  466558553 - 466555373  :  36.626346413s
3198  =  575713847 - 575710649  :  36.761306175s
3240  =  606802529 - 606799289  :  36.799984807s
3330  =  784554983 - 784551653  :  37.014430956s
3486  =  873632213 - 873628727  :  37.121270926s
3588  =  987417437 - 987413849  :  37.25618423s
3900  =  1123404923 - 1123401023  :  37.417362803s
3942  =  1196634239 - 1196630297  :  37.504784859s
4032  =  1247118179 - 1247114147  :  37.565187304s
4176  =  1964330609 - 1964326433  :  38.39652816s
4224  =  2055062753 - 2055058529  :  38.502515034s
4290  =  2160258917 - 2160254627  :  38.625633674s
4626  =  2773400633 - 2773396007  :  39.324109323s
5112  =  4278566879 - 4278561767  :  41.022658954s
n = 4294967295 real = 41.041491885s
$

For comparison: peterSO max 5112 in 41.04s versus The Coder max 4176 in 51.97s.

Coder: max |q-p| 4176 q 1964330609 p 1964326433

Output:

$ javac coder.java && java -Xmx1G coder
Total Primes upto 2147483646 = 105097564 in 11.61 s
(11 - 7) 4 (11.64 s)
<< SNIP >>
(1247118179 - 1247114147) 4032 (34.86 s)
(1964330609 - 1964326433) 4176 (51.97 s)
$
\$\endgroup\$
  • \$\begingroup\$ This looks very impressive. \$\endgroup\$ – user9206 Mar 8 '16 at 11:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy