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.2477822303771973
It still isn't that fast, 2.5 sec for 100 numbers. with my method for checking primality it took a couple of hours...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:
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 seq
Output:> 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