Speeding up the calculation of Liouville function

General Usage Performance
, ,
1

Hi!

In my journey to learn Julia, I was trying to find the smallest counter example of Polya’s conjecture which is a little less than 10^9, while using as little memory as possible. Any tips on how to optimize the algorithm below is highly appreciated.

Context:

Let n be a natural number and let \Omega(n) be the number of prime divisors of n with multiplicity. For instance, \Omega(3)=1, \Omega(6)=2 and \Omega(9)=2. The Liouville \lambda functions is defined as \lambda(n)=(-1)^{\Omega(n)}.

Polya conjectured that \displaystyle L(n)= \sum_{i=1}^n \lambda(k) < 0 for all n>1, which turned out to be false, but the smallest counterexample is huge (906,316,571).

Goal: Given n, find all values of \lambda(k) for 1\le k\le n. Equivalently, find the values of L. This allows us to find the smallest counterexample, as well as other values where L is positive.

My current approach with n=10^9:

Generate a bitvector PrimeQ where PrimeQ[i] is true if i is prime, false if not.

function prime_sieve(N::Integer)
    PrimeQ=trues(N); #assume all numbers are prime
    PrimeQ[1]=false; #Convention: 1 is not prime. 
    sN=Int(floor(sqrt(N)))+1 #go  up to the square root of N 
    for i=1:sN
        if PrimeQ[i]
            PrimeQ[2*i:i:N] .= false 
        end
    end

    return PrimeQ
end

The result of @btime is below

@btime prime_sieve(1000_000_000)
5.582 s (3 allocations: 119.21 MiB)

Next, to store the values of \lambda, we use a bitvector too (to save memory). First, we assign zero (false) to all values of \lambda (false corresponds to 1, true corresponds to -1). Then, we loop over all powers of primes p^j and flip the sign of \lambda(p^j). Here, I used this trick to speed up the bit-flipping part.

function fill_λ(PrimeQ::BitVector)

    N=length(PrimeQ);
    λ= falses(N);#Liouville lambda function

    p=2;
    while p<N
        q=p;
        while q<=N
            @views map!(!,λ[q:q:N],λ[q:q:N])
            q*=p;
        end 
        p=findnext(PrimeQ,p+1)
        if typeof(p)==Nothing #if p is the last one, then findnext returns nothing
            break
        end
    end
    return λ
end

The results of @btime are below

PrimeQ=prime_sieve(1000_000_000)
λ= fill_λ(PrimeQ)
@btime  fill_λ($PrimeQ)
  12.506 s (3 allocations: 119.21 MiB)

After this, finding the smallest counterexample is easy.

Question: Is there anything that I can do to speed up one or both of the functions above? I am not sure if I can optimize the memory usage, but if that’s possible, it would be the cherry on top.

2

For some reason this is faster on my machine…

function fill_λ_bis(PrimeQ::BitVector)

    N=length(PrimeQ);
    λ= falses(N);#Liouville lambda function

    p=2;
    while p<N
        q=p
        while q<=N
            @inbounds for i ∈ q:q:N
                λ[i]=!λ[i]
            end
            q*=p;
        end 
        p=findnext(PrimeQ,p+1)
        if typeof(p)==Nothing #if p is the last one, then findnext returns nothing
            break
        end
    end
    return λ
end
2 Likes
3

Apart from LaurentPlagne’s suggestion, try doing isnothing(p) instead of typeof(p)==Nothing. Apart from being nicer, it makes @code_warntype happier. That means it helps at least at an early compilation stage.

BTW: you don’t have to put a semicolon at the end of the line

2 Likes
4

@inbounds actually speeds up the code by 10%, thanks!

1 Like
5

I see, thanks! As for the semi-colon, it is an old habit from my Matlab days :smiley:.

2 Likes