python - 波利亚猜想的反例

Question

Polya's conjecture是一个数学猜想，假设第一个 (-1)^(Omega(n)) 的和始终为负或零，其中 Omega(n) 是 n 的重数素因数的数量。

反例是 906316571，是五十年前发现的。我想知道他们是怎么找到它的，因为这需要大量的时间，我尝试优化我的 python 算法，但仍然需要大量的时间，你能帮我优化它吗？

这是我的代码(我使用了内存)

 >>> class Memoize:
def __init__(self, f):
    self.f = f
    self.memo = {}
def __call__(self, *args):
    if not args in self.memo:
        self.memo[args] = self.f(*args)
    return self.memo[args]

 >>> def sieve(m):
n=m+1;
s=[];
for i in range(2,n):
    s.append(i);
k=0;
while k<len(s):
    for i in range(2,int(n/s[k])+1):
        x=i*s[k];
        if s.count(x)==1:
            s.remove(x);
    k=k+1;
return s;
>>> s=sieve(100000);
>>> def omega(n):
k=0;
if n==1:
    return 0;
else :
    while k<len(s) and n%s[k]!=0 :
        k=k+1;
    if k<len(s):
        return omega(int(n/s[k]))+1;
    else :
        return 1;
>>> omega=Memoize(omega)
>>> def polya(n):
h=omega(n);
if n==1:
    return 0;
else :
    if omega(n)%2==0:
        return polya(n-1)+1;
    else :
        return polya(n-1)-1;
>>> polya=Memoize(polya);
>>> while polya(k)<=0 :
k=k+1;

Answer 1

正如 chepner 告诉您的那样，最初的 1958 年证明并不是用蛮力完成的。它也没有透露打破规则的最小数字，它是在1980年才发现的。我根本没有研究过这个案例，但1980年的证明可能是用计算机完成的。这更多的是可用 RAM 量的问题，而不是处理速度的问题。

然而，对于现代计算机来说，用蛮力解决这个问题应该不会太困难。 Python 在这里并不是最好的选择，但仍然可以在合理的时间内找到这些数字。

import numpy as np
import time

max_number = 1000000000

# array for results
arr = np.zeros(max_number, dtype='int8')

# starting time
t0 = time.time()

# some tracking for the time spent
time_spent = []

# go through all possible numbers up to the larges possible factor
for p in range(2, int(np.sqrt(max_number))):
    # if the number of factors for the number > 0, it is not a prime, jump to the next
    if arr[p] > 0:
        continue
    # if we have a prime, we will have to go through all its powers < max_number
    n = p
    while n < max_number:
         # increment the counter at n, 2n, 3n, ...
        arr[n::n] += 1
        # take the next power
        n = n * p
    # track the time spent

print "Time spent calculating the table of number of factors: {0} s".format(time.time()-t0)

# now we have the big primes left, but their powers are not needed any more
# they need to be taken into account up to max_number / 2
j = 0
for p in range(p + 1, (max_number + 1) / 2):
    if arr[p] > 0:
        continue
    arr[p::p] += 1
    if j % 10000 == 0:
        print "{0} at {1} s".format(p, time.time()-t0)
    j += 1

print "Primes up to {0} done at {1} s".format(p, time.time()-t0)
# now we have only big primes with no multiples left, they will be all 1
arr[arr == 0] = 1

print "Factor table done at {0} s".format(time.time() - t0)

# calculate the odd/even balance, note that 0 is not included and 1 has 0 factors
cumulative = np.cumsum(1 - 2 * (arr[1:] & 1), dtype='int32')
print "Total time {0} s".format(time.time()-t0)

这不是最快或最优化的函数，其背后的数学原理应该非常明显。在我的运行在一个核心上的机器 (i7) 中，计算最多 1 x 10^9 的质因数的完整表需要大约 2800 秒。 (但要注意，如果没有 64 位 python 和至少 8 GB RAM，请勿尝试此操作。累积和表消耗 4 GB。)

为了证明上述函数至少工作得很好，下面是有趣区域的图:

由于第一个数字存在一些问题，上面代码给出的结果略有偏差。要获取官方总结的 Liouville lambda，请使用 cumulative[n-1] + 2。对于问题中提到的数字 (906 316 571)，结果是 cumulative[906316570] + 2 等于 829，这是该区域的最大值。