Find the largest gap between good primes
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 < 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 < 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.
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)
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