RPython (PyPy 4.0.1), 4032
RPython jest ograniczonym podzbiorem Pythona, który można przetłumaczyć na C, a następnie skompilować za pomocą RPython Toolchain. Jego wyraźnym celem jest pomoc w tworzeniu tłumaczy językowych, ale można go również używać do kompilacji prostych programów.
Aby skompilować, pobierz bieżące źródło PyPy (PyPy 4.0.1) i uruchom następujące czynności:
$ pypy /pypy-4.0.1-src/rpython/bin/rpython --opt=3 good-primes.py
Wynikowy plik wykonywalny zostanie nazwany good-primes-club podobny w bieżącym katalogu roboczym.
Uwagi dotyczące wdrażania
Generator liczb pierwszych primesjest nieograniczonym sito Eratostenesa, które używa koła, aby uniknąć wielokrotności 2 , 3 , 5 lub 7 . Wywołuje się także rekurencyjnie, aby wygenerować następną wartość do użycia przy oznaczaniu. Jestem bardzo zadowolony z tego generatora. Profilowanie linii ujawnia, że dwie najwolniejsze linie to:
37> n += o
38> if n not in sieve:
więc nie sądzę, że jest wiele miejsca na ulepszenia, inne niż być może przy użyciu większego koła.
W celu sprawdzenia „dobroci”, najpierw wszystkie czynniki dwóch są usuwane z n-1 , używając nieco kręcącego się hacka, aby znaleźć największą potęgę dwóch, która jest dzielnikiem (n-1 & 1-n). Ponieważ p-1 jest koniecznie nawet dla dowolnej liczby pierwszej p> 2 , wynika z tego, że 2 musi być jednym z odrębnych czynników pierwszych. To, co pozostaje, jest wysyłane do is_prime_powerfunkcji, która robi to, co sugeruje jej nazwa. Sprawdzanie, czy wartość jest siłą pierwszą, jest „prawie wolne”, ponieważ odbywa się to jednocześnie z kontrolą pierwotności, przy co najwyżej O (log p n) operacjach, gdzie p jest najmniejszym pierwszym czynnikiem n. Podział próbny może wydawać się nieco naiwny, ale według moich testów jest to najszybsza metoda dla wartości mniejszych niż 2 32 . Oszczędzam trochę, ponownie wykorzystując koło z sita. W szczególności:
59> while p*p < n:
60> for o in offsets:
poprzez iterację na kole o długości 48, p*p < nczek zostanie pominięty tysiące razy, przy niskiej, niskiej cenie nie większej niż 48 dodatkowych operacji modulo. Pomija również ponad 77% wszystkich kandydatów, a nie 50%, biorąc tylko szanse.
Ostatnie kilka wyników to:
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
Kod jest również prawidłowym Pythonem i powinien osiągnąć 3588 ~ 3900, gdy jest uruchamiany z najnowszym interpreterem PyPy.
# 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
To przesłanie różni się nieco od innych opublikowanych do tej pory, ponieważ nie sprawdza wszystkich dobrych liczb pierwszych, ale wykonuje stosunkowo duże skoki. Wadą tego jest to, że nie można użyć sit [Stoję skorygowany?] , Więc trzeba polegać całkowicie na testach pierwotności, które w praktyce są nieco wolniejsze. Istnieje również szczęśliwy środek między stopą wzrostu a liczbą wartości sprawdzanych za każdym razem. Mniejsze wartości są znacznie szybsze do sprawdzenia, ale większe wartości mają większe luki.
Aby uspokoić bogów matematyki, postanowiłem zastosować sekwencję podobną do Fibonacciego, mając następny punkt początkowy jako sumę dwóch poprzednich. Jeśli po sprawdzeniu 10 par nie zostaną znalezione żadne nowe rekordy, skrypt przejdzie do następnego.
Ostatnie kilka wyników to:
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
Podczas kompilacji używane są 64-bitowe liczby całkowite, chociaż w kilku miejscach zakłada się, że można dodać dwie liczby całkowite bez przepełnienia, więc w praktyce można użyć tylko 63 liczb całkowitych. Po osiągnięciu 62 bitów znaczących wartość prądu jest dwukrotnie zmniejszana o połowę, aby uniknąć przepełnienia obliczeń. Rezultatem jest to, że tasuje skryptów poprzez wartości na 2 60 - 2 62 przedziale. Nieprzekraczanie natywnej precyzji liczb całkowitych powoduje również, że skrypt jest interpretowany szybciej.
Do potwierdzenia tego wyniku można użyć następującego skryptu PARI / GP:
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)