A counter example for Polya's conjecture

Question

Polya's conjecture is a mathematical conjecture that suppose that the sum of the first (-1)^(Omega(n)) where Omega(n) is the number of prime divisors of n with multiplicity, is always negative or zero.

A counter example is 906316571, was found fifty years ago. I wonder how could they found it because it takes a massive amount of time, I tried to optimize my python algorithm but it still takes a massive time, Can you help me optimize it ?

Here's my code (I used memoization)

 >>> 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

As chepner told you, the original year 1958 proof was not done with brute force. Neither did it reveal the smallest number to break the rule, it was only found in 1980. I have not studied the case at all, but the 1980 proof may be have been done with a computer. It is more a question of the amount of RAM available, not a question of processing speed as such.

However, with modern computers it should not be overtly difficult to approach the problem with brute force. Python is not the best alternative here, but still the numbers can be found in a reasonable time.

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)

This is not the fastest or the most optimized function, the mathematics behind this should be quite obvious. In my machine (i7) running on one core it takes approximately 2800 seconds to calculate the full table of number of prime factors up to 1 x 10^9. (But beware, do not try this without 64-bit python and at least 8 GB RAM. The cumulative sum table consumes 4 GB.)

To prove that the above function works at least quite well, here is a plot of the interesting area:

Due to some problems with the first numbers, the results given by the code above are slightly off. To get the official summatory Liouville lambda, use cumulative[n-1] + 2. For the number mentioned in the question (906 316 571) the result is cumulative[906316570] + 2 equalling to 829, which is the maximum value in the region.