# Speeding up the calculation of Liouville function

**URL:** <https://discourse.julialang.org/t/speeding-up-the-calculation-of-liouville-function/105447>\
**Category:** Performance\
**Tags:** question, performance, boolean\
**Created:** [October 26, 2023, 4:05pm UTC](https://discourse.julialang.org/t/speeding-up-the-calculation-of-liouville-function/105447 "2023-10-26T16:05:09Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![HMegh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hmegh/32/216684_2.png) [@HMegh](https://discourse.julialang.org/u/HMegh)\
**Post date:** [October 26, 2023, 4:05pm UTC](https://discourse.julialang.org/t/speeding-up-the-calculation-of-liouville-function/105447/1 "2023-10-26T16:05:09Z")

</div>

Hi!

In my journey to learn Julia, I was trying to find the smallest counter example of [Polya’s conjecture](https://en.wikipedia.org/wiki/P%C3%B3lya_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.

```julia-auto
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

```julia-auto
@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](https://discourse.julialang.org/t/fastest-way-to-negate-some-booleans-in-an-array/89261/4) to speed up the bit-flipping part.

```julia-auto
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

```julia-auto
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.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [October 26, 2023, 5:11pm UTC](https://discourse.julialang.org/t/speeding-up-the-calculation-of-liouville-function/105447/2 "2023-10-26T17:11:02Z")

</div>

For some reason this is faster on my machine…

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

```

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [October 26, 2023, 8:27pm UTC](https://discourse.julialang.org/t/speeding-up-the-calculation-of-liouville-function/105447/3 "2023-10-26T20:27:58Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![HMegh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hmegh/32/216684_2.png) [@HMegh](https://discourse.julialang.org/u/HMegh)\
**Post date:** [October 27, 2023, 5:13pm UTC](https://discourse.julialang.org/t/speeding-up-the-calculation-of-liouville-function/105447/4 "2023-10-27T17:13:09Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![HMegh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hmegh/32/216684_2.png) [@HMegh](https://discourse.julialang.org/u/HMegh)\
**Post date:** [October 27, 2023, 5:15pm UTC](https://discourse.julialang.org/t/speeding-up-the-calculation-of-liouville-function/105447/5 "2023-10-27T17:15:20Z")

</div>

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