Mar-14-2021, 08:57 AM
I implemented a method to generate fortunate numbers but it was too slow due to my naive implementation of checking primes and to find the smallest prime larger than a certain number, which is required for fortunate number generation. So I looked around and found an astonishingly quick method on codegolf.stackexchange where "fortunately" the winner, "primo", used Python. It is a probability based method that speeds up primality tests several magnitudes. It is still a little slow if you generate many numbers. A hundred is OK but with your ideas one probably uses many more. In any case, here goes:
So I tried to tweak it to make it faster but with this final version of the function "fortunate" it still takes 2.16 sec:
from time import time
def fortunate(n):
# https://en.wikipedia.org/wiki/Fortunate_number
pprod = 2
lastp = 2
seq = []
while len(seq) < n:
m = 2
while not is_prime(pprod + m):
m +=1
seq.append(m)
lastp = next_prime(lastp)
pprod *= lastp
return seq
#### ---- "primos" method for finding the next prime ---- ####
# see: https://codegolf.stackexchange.com/questions/10701/fastest-code-to-find-the-next-prime
# legendre symbol (a|m)
# note: returns m-1 if a is a non-residue, instead of -1
def legendre(a, m):
return pow(a, (m-1) >> 1, m)
# strong probable prime
def is_sprp(n, b=2):
d = n-1
s = 0
while d&1 == 0:
s += 1
d >>= 1
x = pow(b, d, n)
if x == 1 or x == n-1:
return True
for r in range(1, s):
x = (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):
P = 1
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 = 0
while s > 0:
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 > 0:
if t&1 == 1:
# U, V of n+1
U, V = ((U + V) * inv_2)%n, ((D*U + V) * inv_2)%n
q = (q * Q)%n
t -= 1
else:
# U, V of n*2
U, V = (U * V)%n, (V * V - 2 * q)%n
q = (q * q)%n
t >>= 1
# double s until we have the 2**r*sth Lucas number
while r > 0:
U, V = (U * V)%n, (V * V - 2 * q)%n
q = (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 = set([
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]
max_int = 2147483647
# an 'almost certain' primality check
def is_prime(n):
if n < 212:
return n in small_primes
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:
i = 211
while i*i < n:
for o in offsets:
i += o
if n%i == 0:
return False
return True
# Baillie-PSW
# this is technically a probabalistic test, but there are no known pseudoprimes
if not is_sprp(n): return False
a = 5
s = 2
while legendre(a, n) != n-1:
s = -s
a = s-a
return is_lucas_prp(n, a)
# 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:
while True:
if n in small_primes:
return n
n += 2
# find our position in the sieve rotation via binary search
x = int(n%210)
s = 0
e = 47
m = 24
while m != e:
if indices[m] < x:
s = m
m = (s + e + 1) >> 1
else:
e = m
m = (s + e) >> 1
i = int(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
#### ---- End of "primos" method for finding the next prime ---- ####
start = time()
print(fortunate(100))
print(time()-start)And the output:Output:[serafim:~/pyprogs]
> python fortunate.py
[3, 5, 7, 13, 23, 17, 19, 23, 37, 61, 67, 61, 71, 47, 107, 59, 61, 109, 89, 103, 79, 151, 197, 101, 103, 233, 223, 127, 223, 191, 163, 229, 643, 239, 157, 167, 439, 239, 199, 191, 199, 383, 233, 751, 313, 773, 607, 313, 383, 293, 443, 331, 283, 277, 271, 401, 307, 331, 379, 491, 331, 311, 397, 331, 353, 419, 421, 883, 547, 1381, 457, 457, 373, 421, 409, 1061, 523, 499, 619, 727, 457, 509, 439, 911, 461, 823, 613, 617, 1021, 523, 941, 653, 601, 877, 607, 631, 733, 757, 877, 641]
2.2477822303771973So I tried to tweak it to make it faster but with this final version of the function "fortunate" it still takes 2.16 sec:
def fortunate(n):
# https://en.wikipedia.org/wiki/Fortunate_number
pprod = 2
lastp = 2
seq = [0]*n
k = 0
while seq[-1] == 0:
nextp = next_prime(pprod + 2)
seq[k] = nextp-pprod
lastp = next_prime(lastp)
pprod *= lastp
k += 1
return seqOutput:> python fortunate.py
[3, 5, 7, 13, 23, 17, 19, 23, 37, 61, 67, 61, 71, 47, 107, 59, 61, 109, 89, 103, 79, 151, 197, 101, 103, 233, 223, 127, 223, 191, 163, 229, 643, 239, 157, 167, 439, 239, 199, 191, 199, 383, 233, 751, 313, 773, 607, 313, 383, 293, 443, 331, 283, 277, 271, 401, 307, 331, 379, 491, 331, 311, 397, 331, 353, 419, 421, 883, 547, 1381, 457, 457, 373, 421, 409, 1061, 523, 499, 619, 727, 457, 509, 439, 911, 461, 823, 613, 617, 1021, 523, 941, 653, 601, 877, 607, 631, 733, 757, 877, 641]
2.1603410243988037