这个新答案将
big_n
的总时间缩短到
4分6秒。
对我的
初始答案进行分析,发现以下事实:
- 总时间:1小时42分33秒
- 分解数字所花费的时间:几乎占据了全部时间
相比之下,从
3
到
sqrt(2*N - 1)
生成所有素数仅需
38.5秒(使用Atkin筛法)。
因此,我决定尝试一种版本,其中我们将所有数字
m
生成为已知质数乘积。也就是说,生成器产生数字本身以及涉及的不同质因数。无需因式分解。
结果仍然是500_000_000_002_841
,与@Koder发现的相差4。我还不知道问题来自哪里。编辑:在修正
xmax
边界后(
isqrt(2 * N-m ** 2)
而不是
isqrt(2 * N-m ** 2-1)
,因为我们确实希望包括直角三角形的斜边等于
N
),我们现在得到了正确的结果。
质数生成器的代码包含在最后。基本上,我使用了Atkin筛法,并将其适应(没有花费太多时间)Python。我相当确定它可以加速(例如使用numpy
,甚至可能使用numba
)。
为了从质数生成整数(我们知道这是可能的,多亏了
算术基本定理),我们只需遍历所有可能的积
prod(p_i**k_i)
,其中
p_i
是第
i
个质数,
k_i
是任何非负整数。
最简单的公式是递归公式:
def gen_ints_from_primes(p_list, upto):
if p_list and upto >= p_list[0]:
p, *p_list = p_list
pk = 1
p_tup = tuple()
while pk <= upto:
for q, p_distinct in gen_ints_from_primes(p_list, upto=upto // pk):
yield pk * q, p_tup + p_distinct
pk *= p
p_tup = (p, )
else:
yield 1, tuple()
不幸的是,我们很快就会遇到内存限制(和递归限制)。因此,这里有一个非递归版本,除了素数列表本身之外不使用额外的内存。基本上,当前值q
(正在生成的整数)和列表中的索引是我们生成下一个整数所需的所有信息。当然,这些值未排序,但这并不重要,只要它们都被覆盖即可。
def rem_p(q, p, p_distinct):
q0 = q
while q % p == 0:
q //= p
if q != q0:
if p_distinct[-1] != p:
raise ValueError(f'rem({q}, {p}, ...{p_distinct[-4:]}): p expected at end of p_distinct if q % p == 0')
p_distinct = p_distinct[:-1]
return q, p_distinct
def add_p(q, p, p_distinct):
if len(p_distinct) == 0 or p_distinct[-1] != p:
p_distinct += (p, )
q *= p
return q, p_distinct
def gen_prod_primes(p, upto=None):
if upto is None:
upto = p[-1]
if upto >= p[-1]:
p = p + [upto + 1]
q = 1
i = 0
p_distinct = tuple()
while True:
while q * p[i] <= upto:
i += 1
while q * p[i] > upto:
yield q, p_distinct
if i <= 0:
return
q, p_distinct = rem_p(q, p[i], p_distinct)
i -= 1
q, p_distinct = add_p(q, p[i], p_distinct)
例子-
>>> p_list = list(primes(20))
>>> p_list
[2, 3, 5, 7, 11, 13, 17, 19]
>>> sorted(gen_prod_primes(p_list, 20))
[(1, ()),
(2, (2,)),
(3, (3,)),
(4, (2,)),
(5, (5,)),
(6, (2, 3)),
(7, (7,)),
(8, (2,)),
(9, (3,)),
(10, (2, 5)),
(11, (11,)),
(12, (2, 3)),
(13, (13,)),
(14, (2, 7)),
(15, (3, 5)),
(16, (2,)),
(17, (17,)),
(18, (2, 3)),
(19, (19,)),
(20, (2, 5))]
如您所见,我们不需要分解任何数字,因为它们很方便地与涉及的不同质数一起出现。
要仅获取奇数,请从素数列表中删除2
:
>>> sorted(gen_prod_primes(p_list[1:]), 20)
[(1, ()),
(3, (3,)),
(5, (5,)),
(7, (7,)),
(9, (3,)),
(11, (11,)),
(13, (13,)),
(15, (3, 5)),
(17, (17,)),
(19, (19,))]
为了利用这个数字和因子的展示方式,我们需要对原始答案中给出的函数进行一些修改:
def phi(n, upto=None, p_list=None):
if upto is None or upto > n:
upto = n
if p_list is None:
p_list = list(distinct_factors(n))
if upto < n:
cnt = upto
for q in products_of(p_list, upto):
cnt += upto // q if q > 0 else -(upto // -q)
return cnt
cnt = n
for p in p_list:
cnt = cnt * (p - 1) // p
return cnt
通过以上内容,我们现在可以重新编写计数函数:
def pt_count_m(N):
mmax = isqrt(2*N - 1)
p_list = list(primes(mmax))[1:]
for m, p_distinct in gen_prod_primes(p_list, upto=mmax):
if m < 3:
continue
xmax = isqrt(2*N - m*m)
cnt_m = phi(m+1, upto=xmax, p_list=(2,) + tuple(p_distinct))
if cnt_m > 0:
yield m, cnt_m
def pt_count(N, progress=False):
mmax = isqrt(2*N - 1)
it = pt_count_m(N)
if progress:
it = tqdm(it, total=(mmax - 3 + 1) // 2)
return sum(cnt_m for m, cnt_m in it)
现在:
%timeit pt_count(100_000_000)
31.1 ms ± 38.9 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit pt_count(1_000_000_000)
104 ms ± 299 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%%time
big_n = 3_141_592_653_589_793
N = big_n
res = pt_count(N)
CPU times: user 4min 5s, sys: 662 ms, total: 4min 6s
Wall time: 4min 6s
>>> res
500000000002845
附录:阿特金筛法
正如承诺的那样,这是我版本的阿特金筛法。它肯定可以加速。
def primes(limit):
sqrtLimit = isqrt(limit) + 1
is_prime = [False, False, True, True, False] + [False for _ in range(5, limit + 1)]
for x in range(1, sqrtLimit):
x2 = x * x
for y in range(1, sqrtLimit):
y2 = y*y
n = 4 * x2 + y2
if n <= limit and (n % 12 == 1 or n % 12 == 5): is_prime[n] ^= True
n = 3 * x2 + y2
if n <= limit and (n % 12 == 7): is_prime[n] ^= True
n = 3*x2-y2
if n <= limit and x > y and n % 12 == 11: is_prime[n] ^= True
for n in range(5, sqrtLimit):
if is_prime[n]:
sqN = n**2
for i in range(1, int(limit/sqN) + 1):
k = i * sqN
is_prime[k] = False
for i, truth in enumerate(is_prime):
if truth: yield i
numba
和我第二个答案的变体。我会在有时间的时候添加一些说明。不确定是将其添加到第二个答案中还是写一个第三个答案。你有什么想法? - Pierre D