# Longest Prime Sums

Sandbox

There are special sets S of primes such that $$\\sum\limits_{p\in S}\frac1{p-1}=1\$$. In this challenge, your goal is to find the largest possible set of primes that satisfies this condition.

Input: None

Output: A set of primes which satisfies the conditions above.

This challenge is a , where your score is the size of the set you output.

# Examples of Valid Outputs and their scores

{2} - 1
{3,5,7,13} - 4
{3,5,7,19,37} - 5
{3,5,7,29,31,71} - 6


Important: This challenge doesn't require any code, simply the set. However, it would be helpful to see how you generated your set, which probably involves some coding.

Edit: Found the math.stackexchange post that inspired this thanks to @Arnauld

• This challenge doesn't require any code - would this be better asked in puzzling.SE? – Digital Trauma Dec 11 '19 at 23:09
• @DonThousand Is it correct that you require the solution to be a set (e.g. all elements are distinct)? Or is a list of primes (with recurring values) also acceptable? I (possibly wrongly) assumed in my answer that repetitions were allowed. If not, I'll change my answer. – agtoever Dec 12 '19 at 11:30
• @agtoever It's a set. It must be distinct. I see you've already made that change, so that's good. – Don Thousand Dec 12 '19 at 14:35
• I agree with @DigitalTrauma that this should be on puzzling.SE. They specifically have a 'no-computers' tag which to me implies that puzzles requiring computers are acceptable there. – Sam Dean Dec 12 '19 at 15:55
• @SamDean No, that's not true. Open ended challenges like this are strictly off topic there. Problems must have a single definite answer to be posted there. Please read their FAQ – Don Thousand Dec 12 '19 at 15:58

# Score 100 8605

I used an algorithm that starts with one solution and repeatedly tries to split a prime $$\p\$$ in the solution into two other primes $$\q_ 1\$$ and $$\q_ 2\$$ that satisfy $$\\frac1{p-1} = \frac1{q_1-1}+\frac1{q_2-1}\$$.

It is known (and can be quickly checked) that the positive integer solutions to $$\\frac1n = \frac1x + \frac1y\$$ are in one-to-one correspondence with factorizations $$\n^2 = f_ 1 f_ 2\$$, the correspondence being given by $$\x = n + f_ 1\$$, $$\y = n + f_ 2\$$. We can search through the factorizations of $$\(p-1)^2\$$ to see if any of them yielded a solution where both new denominators $$\x,y\$$ were one less then primes; if so, then $$\p\$$ can be replaced by $$\q_1=x+1\$$, $$\q_2=y+1\$$. (If there were multiple such factorizations, I used the one where the minimum of $$\q_1\$$ and $$\q_2\$$ was smallest.)

I started with this seed solution of length 44:

seed = {3, 7, 11, 23, 31, 43, 47, 67, 71, 79, 103, 131, 139, 191, 211, 239, 331, 419, 443, 463, 547, 571, 599, 647, 691, 859, 911, 967, 1103, 1327, 1483, 1871, 2003, 2311, 2347, 2731, 3191, 3307, 3911, 4003, 4931, 6007, 6091, 8779}


This seed was found using an Egyptian fraction solver that a former research student of mine, Yue Shi, coded in ChezScheme. (The primes $$\p\$$ involved all have the property that $$\p-1\$$ is the product of distinct primes less than 30, which increased the likelihood of Shi's program finding a solution.)

The following Mathematica code continually updates a current solution by looking at its primes one by one, trying to split them in the manner described above. (The number 1000 in the third line is an arbitrary stopping point; in principle one could let the algorithm run forever.)

solution = seed;
j = 1; (* j is the index of the element of the solution that we'll try to split *)

While[j <= 1000 && j <= Length[solution],
currentP = solution[[j]];
allDivisors = Divisors[(currentP - 1)^2];
allFactorizations = {#, (currentP - 1)^2/#} & /@
Take[allDivisors, Floor[Length[allDivisors]/2]];
allSplits = currentP + allFactorizations;
goodFactorizations = Select[allSplits,
And @@ PrimeQ[#] && Intersection[#, solution] == {} &];
If[goodFactorizations == {},
j++,
solution = Union[Complement[solution, {currentP}], First@goodFactorizations]
]
]


The code above yields a solution of length 4126, whose largest element is about $$\8.7\times10^{20}\$$; by the end, it was factoring integers $$\(p-1)^2\$$ of size about $$\8.8\times10^{21}\$$.

In practice, I ran the code several times, using the previous output as the next seed in each case and increasing the cutoff for j each time; this allowed for the recovery of some small prime splits that had become non-redundant thanks to previous splitting, which somewhat mitigated the size of the integers the algorithm factored.

The final solution, which took about an hour to obtain, is too long to fit in this answer but has been posted online. It has length 8605 and largest element about $$\4.62\times10^{19}\$$.

Various runs of this code consistently found that the length of the solution was about 3–4 times as long as the set of primes that had been examined for splitting. In other words, the solution was growing much faster than the code scanned through the initial elements. It seems likely that this behavior would continue for a long time, yielding some gargantuan solutions.

• Reading this discouraged me from competing, because you know your stuff so well – Mark Jeronimus Dec 12 '19 at 16:35
• Oh I don't want anyone to be discouraged! Note that agtoever blew my original answer out of the water :) – Greg Martin Dec 12 '19 at 19:17
• @GregMartin In my opinion, you deserve the accepted answer. Not only have you discovered a much larger set, but also an even more effective and efficient algorithm. Well done. Kudos. – agtoever Dec 12 '19 at 19:38
• @agtoever Thanks for the kind words :) Your answer is very worthy and also was the first to have the code up; speed of posting is a consideration factor and I'm content for that to be the case. – Greg Martin Dec 12 '19 at 20:12
• @JollyJoker I think it's an extremely interesting question whether a given seed can generate an infinite sequence of splits. My heuristics tell me that a randomly chosen large prime is unlikely to have such splits; however, primes such that p–1 is very composite have a better chance to have such splits, and it might be the case that this algorithm tends to produce primes of that form. – Greg Martin Dec 13 '19 at 19:51

# Score 263385425 426 with only primes < 1.000.000 (was: non-competitive, now it is; score can be increased by running the program longer)

I followed the same path as Wheat Wizard: iteratively search for primes in the solution that can be replaced with a longer list of primes with the same result. I wrote that Python program that does exactly this. It starts with solution S = {2} and than iterates of all elements of that solution and tries to find a decomposition of that prime for which 1/(p-1) = sum(1/(q-1) for all q in the decomposition.

After I realized that S should be a set (and not a list), I altered the program to take this into account. I also added a ton of performance optimizations. The solution of 263 came up within 200 seconds or so (running under pypy3), but if you let it running, it steadily keeps coming with additional (longer) solutions.

Current best solution (425 elements, with all primes < 1M, calculated in ~ 15 min.):

S = [3, 5, 13, 19, 29, 37, 103, 151, 241, 281, 409, 541, 577, 593, 661, 701, 751, 1297, 1327, 2017, 2161, 2251, 2293, 2341, 2393, 2521, 2593, 2689, 2731, 3061, 3079, 3329, 3361, 3457, 6301, 6553, 7057, 7177, 7481, 7561, 8737, 9001, 9241, 9341, 10501, 11617, 12097, 12547, 14281, 14449, 14561, 15121, 17761, 17851, 18217, 18481, 20593, 21313, 22441, 23189, 23761, 24571, 26041, 26881, 28351, 28513, 29641, 30241, 36529, 37441, 46993, 49921, 51169, 57331, 58109, 58313, 58369, 58831, 59659, 60737, 60757, 61001, 61381, 61441, 61561, 61609, 63067, 63601, 64513, 64901, 65053, 65089, 65701, 65881, 66301, 66931, 67049, 69389, 69941, 70181, 72161, 72481, 72577, 72661, 73061, 73699, 74521, 77521, 78241, 79693, 81181, 86951, 88741, 90631, 98011, 100297, 102181, 107641, 108991, 109201, 109537, 114913, 117841, 118429, 121993, 122761, 123001, 124561, 127601, 128629, 130073, 130969, 131561, 133387, 133813, 138181, 138403, 139501, 146077, 149521, 159457, 160081, 162289, 163543, 166601, 174241, 175891, 176401, 177913, 180181, 182711, 189421, 199921, 201781, 206641, 218527, 223441, 227089, 229739, 234961, 238081, 238141, 238897, 239851, 246241, 250057, 261577, 266401, 267961, 280321, 280837, 280897, 281233, 283501, 283861, 284161, 287233, 288049, 291721, 297601, 299053, 302221, 306853, 309629, 313153, 316681, 322057, 325921, 332489, 342211, 342241, 349981, 352273, 354961, 355321, 360977, 365473, 379177, 390097, 390961, 394717, 395627, 401057, 404251, 404489, 412127, 412651, 416881, 417649, 418027, 424117, 427681, 428221, 428401, 429409, 430921, 434521, 435481, 441937, 443873, 444641, 451441, 453601, 454609, 455149, 459649, 466201, 468001, 473617, 474241, 480737, 481693, 483883, 496471, 498301, 498961, 499141, 499591, 499969, 501601, 501841, 502633, 513067, 514513, 517609, 523261, 524521, 525313, 529381, 538721, 540541, 545161, 550117, 552553, 560561, 562633, 563501, 563851, 568177, 570781, 575723, 587497, 590669, 591193, 599281, 601801, 601903, 604001, 605551, 607993, 609589, 611389, 617401, 621007, 627301, 628561, 628993, 629281, 635449, 637201, 639211, 642529, 645751, 651361, 651857, 653761, 654853, 655453, 657091, 662941, 664633, 667801, 669121, 669901, 670177, 673201, 673921, 675109, 688561, 689921, 691363, 692641, 694033, 695641, 697681, 698293, 700591, 703081, 703561, 705169, 705181, 707071, 709921, 713627, 732829, 735373, 737413, 739861, 742369, 745543, 750121, 750721, 754771, 756961, 757063, 758753, 759001, 760321, 761671, 762721, 766361, 773501, 774181, 776557, 779101, 782461, 784081, 784981, 786241, 788317, 794641, 795601, 797273, 800089, 801469, 808081, 808177, 810151, 813121, 815671, 819017, 823481, 823621, 825553, 831811, 833281, 833449, 836161, 839161, 840911, 846217, 859657, 859861, 860609, 863017, 865801, 869251, 870241, 875521, 876929, 878011, 880993, 884269, 891893, 895681, 898921, 899263, 902401, 904861, 905761, 907369, 908129, 914861, 917281, 917317, 921601, 922321, 923833, 926377, 939061, 941641, 942401, 943009, 943273, 944161, 944821, 944833, 949621, 949961, 950041, 950401, 953437, 953443, 954001, 957349, 957529, 960121, 960961, 963901, 964783, 967261, 967627, 967751, 968137, 971281, 973561, 973591, 984127, 984341, 984913, 986437, 991381, 992941, 994561, 995347, 996001]


Some of the decompositions used:

1/(  2-1) = sum(1/(p-1)) for p in: {(5, 7, 13, 3)}
1/(  3-1) = sum(1/(p-1)) for p in: {(7, 13, 5)}
1/(  5-1) = sum(1/(p-1)) for p in: {(13, 7)}
1/(  7-1) = sum(1/(p-1)) for p in: {(13, 19, 37)}
1/( 13-1) = sum(1/(p-1)) for p in: {(19, 37)}
1/( 19-1) = sum(1/(p-1)) for p in: {(29, 71, 181)}
1/( 29-1) = sum(1/(p-1)) for p in: {(37, 127)}
1/( 37-1) = sum(1/(p-1)) for p in: {(43, 281, 2521)}
1/( 43-1) = sum(1/(p-1)) for p in: {(53, 223, 13469)}
...
1/(8779-1) = sum(1/(p-1)) for p in: {(8969, 739861, 941641)}
1/(8807-1) = sum(1/(p-1)) for p in: {(9001, 773501, 865801)}
1/(8821-1) = sum(1/(p-1)) for p in: {(8941, 657091)}
1/(8941-1) = sum(1/(p-1)) for p in: {(9041, 808177)}
1/(8969-1) = sum(1/(p-1)) for p in: {(9463, 227089, 705169)}
1/(9001-1) = sum(1/(p-1)) for p in: {(9109, 759001)}
1/(9041-1) = sum(1/(p-1)) for p in: {(9241, 417649)}
1/(9109-1) = sum(1/(p-1)) for p in: {(10891, 55661)}
1/(9181-1) = sum(1/(p-1)) for p in: {(9343, 529381)}
1/(9241-1) = sum(1/(p-1)) for p in: {(9341, 863017)}
1/(9341-1) = sum(1/(p-1)) for p in: {(15121, 25219, 784561)}
1/(9343-1) = sum(1/(p-1)) for p in: {(9689, 261577)}
1/(9463-1) = sum(1/(p-1)) for p in: {(10957, 69389)}
1/(9689-1) = sum(1/(p-1)) for p in: {(12457, 43597)}
1/(10333-1) = sum(1/(p-1)) for p in: {(10501, 645751)}
...
1/(131561-1) = sum(1/(p-1)) for p in: {(237361, 295153)}
1/(166601-1) = sum(1/(p-1)) for p in: {(243433, 527851)}
1/(266401-1) = sum(1/(p-1)) for p in: {(386401, 857809)}
1/(355321-1) = sum(1/(p-1)) for p in: {(639577, 799471)}


Python3 code:

import itertools
import sympy
import cProfile
import functools
import bisect
import operator

def sundaram3(max_n):
# Returns a list of all primes under max_n
numbers = list(range(3, max_n + 1, 2))
half = (max_n) // 2
initial = 4
for step in range(3, max_n + 1, 2):
for i in range(initial, half, step):
numbers[i - 1] = 0
initial += 2 * (step + 1)
if initial > half:
return  + list([_f for _f in numbers if _f])

# Precalculate all primes up to a million to speed things up
PRIMES_TO_1M = list(sundaram3(1000000))

def nextprime(number):
# partly precalculated fast version for calculating the
# first (e.g. smallest) prime that is largest than numer
global PRIMES_TO_1M

if number <= PRIMES_TO_1M[-2]:
return PRIMES_TO_1M[bisect.bisect(PRIMES_TO_1M, number)]
return sympy.nextprime(number)

def isprime(number):
# partly precalculated fast version to determine of number is prime
global PRIMES_TO_1M

if number < 1000000:
return number in PRIMES_TO_1M
return sympy.isprime(number)

def upper_limit(prime, length=2):
# Returns the largest prime q in the decomposition of prime with the given
# length such that 1 / (prime - 1) = sum( 1 / (q - 1)) for q in
# set V, with V has the given length.
# ASSUMPTION: all q are unique; this assumption is not validated,
# but for this codegolf, the solution must be a set, so this is safe.
if length == 1:
return prime
nextp = nextprime(prime)
largestprime = (prime * nextp - 2 * prime + 1 ) // (nextp - prime)
if not isprime(largestprime):
largestprime = nextprime(largestprime)
return upper_limit(largestprime, length - 1)

def find_decomposition(prime, length=2):
# Returns a list of primes V = {q1, q2, q3, ...} for which holds that:
# 1 / (prime - 1) = sum(1 / (q - 1)) for q in V.
# Note that there may be more than one V of a given prime, but this
# function returns the first found V (implementation note: the sortest one)
print(f"Searching decomposition for prime {prime} with length {length} in range ({prime}, {upper_limit(prime, length)}]")
prime_range = PRIMES_TO_1M[bisect.bisect(PRIMES_TO_1M, prime) + 1:
bisect.bisect(PRIMES_TO_1M, upper_limit(prime, length) + 1)]
# we only search for combinations of length -1; the last factor is calculated
for combi in itertools.combinations(prime_range, length - 1):
# we find the common factor of prime and all primes in combi
# and use that to calculate the remaining prime. This is faster
# than trying all prime combinations.
factoritems = [-prime + 1] + [c - 1 for c in combi]
factor = -functools.reduce(operator.mul, factoritems)
remainder = - factor / sum(factor // p for p in factoritems) + 1
if remainder == int(remainder) and isprime(remainder) and remainder not in combi:
combi = combi + (int(remainder),)
print(f"Found decomposition: {combi}")
return combi
return None

def find_solutions():
# Finds incrementally long solutions for the set of primes S, for which:
# sum(1/(p-1)) == 1 for p in S.
# We do this by incrementally searching for primes in S that can be
# replaced by longer subsets that have the same value of sum(1/p-1).
# These replacements are stored in a dictionary "decompositions".
# Starting with the base solution S =  and all decompositions,
# you can construct S.
decompositions = {}  # prime: ([decomposition], max tried length)
S = 
old_solution_len = 0

# Keep looping until there are no decompositions that make S longer
while len(S) > old_solution_len:
# Loop over all primes in S to search for decompositions
for p in sorted(set(S)):
# If prime p is not in the decompositions dict, add it
if p not in decompositions:
decompositions[p] = (find_decomposition(p, 2), 2)
# If prime p is in the decompositions dict, but without solution,
# try to find a solution 1 number longer than the previous try
elif not decompositions[p]:
length = decompositions[p] + 1
decompositions[p] = (find_decomposition(p, length), length)
# If prime p is in decompositions and it has a combi
# and the combi is not already in S, replace p with the combi
elif all(p not in S for p in decompositions[p]):
old_solution_len = len(S)
print(f"Removing occurence of {p} with {decompositions[p]}")
S.remove(p)
S.extend(decompositions[p])
S = sorted(S)
break  # break out of the for loop
print(f"Found S with length {len(S)}: S = {S}")
print(f"Decompositions: ")
for prime in sorted(decompositions.keys()):
print(f"    1/({prime:3}-1) = sum(1/(p-1)) for p in: \u007b{decompositions[prime]}\u007d")

def main():
cProfile.run("find_solutions()")

if __name__ == "__main__":
main()

• Using your function, you could decompose the twenty 43s into twenty 47, 491, 33811, increasing the score by 40. Also, the sixteen 73s can go to sixteen times 79, 1093, 6553, adding another 32 to your score – Mathias711 Dec 12 '19 at 10:47
• @LuisMendo altered the code, added solution. To all who downvoted: please let me know you find anything non-competing about this answer. I believe this is now a valid (and sound) answer. I'll edit some code comments later to explain the algorithms used. – agtoever Dec 12 '19 at 12:54
• @Arnauld you are right, I think. Fixed that. – agtoever Dec 12 '19 at 12:55
• Impressive score! – Luis Mendo Dec 12 '19 at 13:18
• Woops, sorry we edited at the same time. Feel free to changes the removing/replacing again, my bad. I added Python highlighting to your answer so it's easier to read the code and comments of your program. – Kevin Cruijssen Dec 12 '19 at 13:38

# Score 3234 36

{5, 7, 11, 13, 17, 23, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 101, 113, 131, 137, 151, 211, 229, 241, 281, 313, 379, 401, 433, 457, 491, 521, 571, 601, 25117, 293362609}


This is an improvement of Arnauld's answer. I just noticed that

$$\ \dfrac{1}{19-1}=\dfrac{1}{73-1}+\dfrac{1}{61-1}+\dfrac{1}{41-1} \$$

But 41 and 61 were already used in Arnauld's answer so I had to then figure out that

$$\ \dfrac{1}{61-1} = \dfrac{1}{151-1} + \dfrac{1}{101-1} \$$ and $$\ \dfrac{1}{41-1} = \dfrac{1}{281-1}+\dfrac{1}{211-1}+\dfrac{1}{151-1}+\dfrac{1}{101-1} \$$

But now I am using 151 and 101 twice. So I spent some time and discovered that

$$\ \dfrac{1}{151-1} = \dfrac{1}{401-1} + \dfrac{1}{241-1} \$$ and $$\ \dfrac{1}{101-1} = \dfrac{1}{601-1} + \dfrac{1}{571-1} + \dfrac{1}{457-1} + \dfrac{1}{229-1} \$$

So now I can just replace the 19 with 71, 151, 101, 151, 211, 229, 241, 281, 401, 457, 471, 601 and the sequence will maintain it's properties.

I also discovered that I can replace 19 with the sequence 71, 151, 101, 151, 211, 241, 241, 281, 401, 433, 541, 601, but that has 241 twice.

After that improvement I also noticed that 79 could be replaced with 521, 313, 131, to increase the size by 2 more. And 73 can be replaced with 113, 379, 433 for another 2.

# Score 22

Just to get the ball rolling.

{5,7,11,13,17,19,23,31,37,41,43,47,53,59,61,67,71,79,137,491,25117,293362609}


Compute the fraction

I suspect that the sequence can be made arbitrary large, but my code is currently too messy and inefficient for anything significantly better than that.

• Most mathematicians suspect that as well, however, that claim is equivalent to a weak form of a long unsolved conjecture. – Don Thousand Dec 12 '19 at 2:36
• @DonThousand Can you link to what is known/suspected mathematically? What is the conjecture you refer to? – Anush Dec 12 '19 at 7:59
• @Anush I'll try to find the paper. It was quite a while ago. – Don Thousand Dec 12 '19 at 8:03

This non-answer expands on the modifications to @GregMartin's answer referenced in this comment there.

The allDivisors to goodFactorizations block can be sped up noticeably by not asking Divisors to factor the square of a number (and also a few other changes). The @GregMartin's original code:

divMethod[p_] := Module[{},
allDivisors = Divisors[(p - 1)^2];
allFactorizations = {#, (p - 1)^2/#} & /@
Take[allDivisors, Floor[Length[allDivisors]/2]];
allSplits = p + allFactorizations;
goodFactorizations = Select[allSplits,
And @@ PrimeQ[#] &]
]


The intersection check at the end is removed, since it is common to both the above and my proposed replacement:

g[p_] := Module[
{pm1s, factorization, divisors},
pm1s = (p - 1)^2;
(* Since the original seed contains only odd
primes, we know one factor of p-1. *)
factorization = FactorInteger[(p - 1)/2];
(* For divisors d1, d2, such that d1 d2 = pm1s,
we require d1+p and d2+p are prime.  p>2, so p
is odd and both d1+p and d2+p are odd, so d1
and d2 are necessarily both even.  This means
we only consider splitting the powers of 2 so
that at least one falls in each divisor; we
want to know the power of 2 in pm1s, minus 1.
*)
If[factorization[[1, 1]] != 2,
pow2m1 = 1,
pow2m1 = 2 factorization[[1, 2]] + 1
];
divisors = Outer[Times,
2^Range[1, pow2m1],
Sequence @@ (
Select[
#[]^Range[0, 2 #[]],
(# < p - 1 &)] & /@
Rest[factorization])];
(* The Join@@ is a hack to deal with Reap's
denormalized output on no-Sow runs.  See
https://mathematica.stackexchange.com/questions/67625/what-is-shorthand-way-of-reap-list-that-may-be-empty-because-of-zero-sow *)
Join @@ Reap[
Map[
(If[#1 < p - 1 && PrimeQ[#2],
(If[PrimeQ[#2],
Sow[{#1, #2}]
] &)[#2, pm1s/#1 + p]
] &)[#, # + p] &,
divisors,
{-1}];
][]
]


Note: the replacement makes no attempt to sort the collection of pairs. The following tests do not demonstrate the potential difference in order of results from divMethod and g.

divMethod
g
(* {} *)
(* {} *)

divMethod
g
(* {{31, 421}, {37, 127}} *)
(* {{31, 421}, {37, 127}} *)

p = NextPrime[10^3]
RepeatedTiming[divMethod[p]]
RepeatedTiming[g[p]]
(* {0.00041, {{1201, 6301}}} *)
(* {0.000487, {{1201, 6301}}} *)

p = NextPrime[10^6]
RepeatedTiming[divMethod[p]]
RepeatedTiming[g[p]]
(* {0.000271, {}} *)
(* {0.0000619, {}} *)

p = NextPrime[10^9]
RepeatedTiming[divMethod[p]]
RepeatedTiming[g[p]]
(* {0.000271, {}} *)
(* {0.000103, {}} *)


Let's collect timing data for sets of 1000 consecutive primes starting at powers of 10 and see what trends we see.

divTiming = Table[{10^k, RepeatedTiming[
p = NextPrime[10^k];
For[count = 1, count <= 1000, count++,
divMethod[p];
p = NextPrime[p];
]
][]}, {k, 3, 22}];

fiTiming = Table[{10^k, RepeatedTiming[
p = NextPrime[10^k];
For[count = 1, count <= 1000, count++,
g[p];
p = NextPrime[p];
]
][]}, {k, 3, 22}];

ListLogLinearPlot[{divTiming, fiTiming},
PlotLegends -> {"Divisors", "FactorInteger"}] ListLogLinearPlot[
Transpose[{
divTiming[[All, 1]],
(divTiming/fiTiming)[[All, 2]]}],
PlotLegends -> {"Divisors/FactorInteger"}] So around 10^7 or 10^8, the Divisors method is noticeably slower than the FactorInteger method. The ratio of times rises to 3-ish for primes around 10^15. There is a dip, for slightly larger primes, but the trend suggests the speed-up will improve as we go to larger primes than tested.

• Wow, this is some really cool analysis. Is this all done in R? – Don Thousand Dec 17 '19 at 1:07
• @DonThousand : This is Mathematica code, as is \@GregMartin's. – Eric Towers Dec 17 '19 at 4:09