Computing the first $n$ values of the Liouville function in linear time

Question

Is it possible to compute the first $n$ values of the Liouville function in linear time? Since we need to output $n$ values we clearly cannot do better than linear time, but the best I can figure out is something like $O(n \cdot \log{\log{n}})$: fill an array of size $n$ with ones, and for each prime power $p^a$, negate the value at the index of each of its multiples. I think it is possible to identify the prime powers as we count up using another table of $O(n)$ bits, essentially the sieve of Eratosthenes counting powers too, but there are still $\sum_{p^a \le n}{\frac{n}{p^a}} = n \cdot \log{\log{n}} + O(n)$ negation operations. Is it possible to do better than this?

Answer 1

Your problem is related to the problem of finding all the primes up to $n$, and it can be solved by similar means.

Although the traditional sieve of Eratosthenes has complexity $O(n \log \log n)$, there are improved versions, which work in $O(n)$ time. It is achieved by crossing out each composite number exactly once. For the purpose of finding primes, these algorithms can even be optimized to $O(n /\log \log n)$ time by the so-called wheel optimization. You can find details in e.g. this paper.

In order to solve your problem, you have to calculate the least prime factor function (denoted by $lpf$) during sieving, not only the primes. When the $lpf$ function is ready, you can compute completely multiplicative functions in $O(n)$ easily by dynamic programming, e.g.: $$ \lambda(x) = \begin{cases} -1 & x = lpf(x) \\ -\lambda \left( \frac{x}{lpf(x)} \right) & \text{otherwise} \end{cases} $$

One of the algorithms which calculates least prime factors in $O(n)$ time is described here. Below you can see its implementation in C++, with primes being a sorted list of all primes found so far, and lpf being an array representing the same-named function.

vector<int> primes;
vector<int> lpf(n, -1);
for (int x = 2; x < n; x++) {
  if (lpf[x] < 0) { //prime found
    lpf[x] = x;
    primes.push_back(x);
  }
  for (int i = 0; i < primes.size(); i++) {
    int p1 = primes[i], y = p1 * x;
    if (p1 > lpf[x] || y >= n)
        break;
    lpf[y] = p1;
  }
}

Consider a number $x$ having least prime factor $p = lpf(x)$. For each prime $p_1 \le p$, it is easy to see that $lpf(p_1 \cdot x) = p_1$. The algorithm simply applies this crossing-out rule for each number $x$ in increasing order. Of course, it needs the sorted list of primes in order to iterate over all the necessary primes $p_1$.

Every composite number $y$ is crossed out exactly once when considering number $x = \frac{y}{lpf(y)}$, so time complexity is $O(n)$. In practice however, it is slower than the traditional $O(n \log \log n)$ sieve, at least for the purpose of finding primes. I guess your approach would also be faster, especially with bitsets.

If you are interested in practical acceleration of your algorithm, you'd better think about memory/cache optimization, instead of improving asymptotic complexity by a double-log factor =)

Answer 2

Actually this cannot be true because you can't get the prime (power)s in $O(n)$:

The sieve of eratosthenes takes at least $O(nloglogn)$ Operations.

Assume we have a list of all the primes up to $n$ given. Then checking each $k$ for membership in the list takes at least (if we only check the first $\pi(\sqrt{k})$ elements) $\pi(\sqrt{k})$ operations each, which amounts to more than $O(1)$ so that's out. Then we could check each number for primality but the fastest known primality test runs in $O(log(n)^6)$, also more than $O(1)$. Since for each number we either have to run it against a list or check it by hand that amounts to all possibilites and we conclude it's impossible to calculate the value for all $k \leq n$ in $O(n)$

Answer 3

The number $N$ can not have more than one prime divisor $p$ that $p^2>N$. Such a divisor there, when $$N>\prod\limits_{\substack{p | N,\\ p^2\leq N}}p^a$$ This allows, similarly to the algorithm "Eratosthenes Sieve" (ES), to limit the enumeration prime divisors of N with ratio $p^2 \leq N$.

Differences algorithm compared with ES as follows:
1. The processing of each factor changes the sign of each corresponding $\lambda$, and each corresponding product multiplied by a factor $p$ (in ES - only the number zero).
2. In addition to the prime factors, involved the same processing passage of their degree $p^a\geq N$. These changes complicate the algorithm proportionally, but do not increase the level of its complexity.

The proposed software implementation (PHP) pre-computed array primes dimension inter = 46340 (the SE algorithm), that allows to calculate the Liouville function for an arbitrary piece of data (val_from, val_to) in the range from 1 to maxval = 2147483647.

The test calculations for 1 000 000 points on the the most costly piece (maxval - 999999, maxval) on weak computer took 20 seconds.

In my opinion, such an algorithm is actually better..

The prorgam is:

ini_set('memory_limit','256M');
$maxint = 2147483647;
 $inter = 46340;      // =[sqrt($maxint)]
 $val_from = 1; //$maxint- 999;
 $val_to = 1000; //$maxint;

function print_p($arr, $width, $from){
  $curr = $from;
  printf("<br>(%05d : %05d)", $curr, ($curr+=$width)-1);
  $prn = $width;
  foreach($arr as $item){
   if ($prn-- == 0){
    $prn = $width-1;
    printf("<br>(%05d : %05d)",$curr, ($curr+=$width)-1);
   }
   printf("&emsp;%2d", $item);
    };
    print("<br>");
}

function base_sieve(){
    global $inter, $p_base;
    $p0 = array(    2, 3, 5, 7, 11,             13, 17, 19, 23, 29,         31, 37, 41, 43, 47,         53, 59, 61, 67, 71,
                    73, 79, 83, 89, 97,         101, 103, 107, 109, 113,    127, 131, 137, 139, 149,    151, 157, 163, 167, 173,
                    179, 181, 191, 193, 197,    199, 211, 223);

    $p_sieve = range(0,$inter-1);
    array_walk($p0, function($item) use(&$p_sieve, $inter) {
        for($j=$item*$item; $j<$inter; $j+=$item) $p_sieve[$j] = 0; 
  });
  $p_sieve[1]=0;
    $p_base = array();
  array_walk($p_sieve, function ($item) use(&$p_base){
        if($item) $p_base[]=$item;
    });
}

function liouville($from, $to){
    global $p_base;
  if(is_int($from) && is_int($to)){
   $lambda = array_fill($from, $to-$from+1, 1);
   $prod = $lambda;
    }else{
        exit(1);
    } 

    $count_base = count($p_base);
    for($key = 0; $key < $count_base-1; $key++){
        $factor = $p_base[$key];
   $factor2 = pow($factor,2);
   if($factor2 > $to) break;
   $item_old = 1;
        for($item = $factor; $item <= $to; $item *= $factor){
            if($item <= $item_old) break;     // overflow test
            $item_old = $item;
            for($j=$to-$to%$item; $j>= max($from, $factor2); $j-=$item) {
     $prod[$j] *= $factor;
     $lambda[$j] = -$lambda[$j];
            }
        }
    };

    printf("<br>");
    for($key = $from; $key <= $to; $key++){
   //printf("prod[%d]=%d ", $key, $prod[$key]);
   if($prod[$key] != $key) $lambda[$key] *= -1; 
  }
  return $lambda;
}

$p_base=0;
 base_sieve();
 print ("maxint=$maxint inter=$inter");
 printf("<br> base_sieve: p_first = %d p_last = %d", reset($p_base), end($p_base));
 printf("<br> liouville: val_from = %d&emsp;val_to = %d<br>", $val_from, $val_to);
 $lambda = liouville($val_from, $val_to);
print_p($lambda,25,$val_from);

Results:

maxint=2147483647 inter=46340
base_sieve: p_first = 2 p_last = 46337
liouville: val_from = 1 val_to = 1000


(00001 : 00025)  1 -1 -1  1 -1  1 -1 -1  1  1 -1 -1 -1  1  1  1 -1 -1 -1 -1  1  1 -1  1  1
(00026 : 00050)  1 -1 -1 -1 -1 -1 -1  1  1  1  1 -1  1  1  1 -1 -1 -1 -1 -1  1 -1 -1  1 -1
(00051 : 00075)  1 -1 -1  1  1  1  1  1 -1  1 -1  1 -1  1  1 -1 -1 -1  1 -1 -1 -1 -1  1 -1
(00076 : 00100) -1  1 -1 -1 -1  1  1 -1  1  1  1  1  1 -1  1  1 -1  1  1  1  1 -1 -1 -1  1
(00101 : 00125) -1 -1 -1  1 -1  1 -1 -1 -1 -1  1 -1 -1 -1  1 -1 -1  1  1 -1  1  1  1 -1 -1
(00126 : 00150)  1 -1 -1  1 -1 -1  1  1  1  1  1 -1 -1 -1  1  1  1  1  1  1  1 -1 -1 -1  1
(00151 : 00175) -1  1 -1 -1  1  1 -1  1  1  1  1 -1 -1 -1 -1  1 -1 -1  1 -1 -1 -1 -1 -1 -1
(00176 : 00200) -1  1  1 -1 -1 -1 -1  1  1  1 -1  1 -1  1 -1 -1 -1 -1  1 -1  1 -1  1 -1 -1
(00201 : 00225)  1  1  1  1  1  1 -1 -1  1  1 -1 -1  1  1  1  1  1  1  1  1  1 -1 -1  1  1
(00226 : 00250)  1 -1  1 -1 -1 -1  1 -1  1  1 -1  1 -1 -1  1 -1 -1 -1 -1 -1 -1  1  1  1  1
(00251 : 00275) -1 -1  1  1 -1  1 -1 -1  1  1 -1  1 -1 -1  1 -1  1 -1 -1 -1 -1 -1 -1  1 -1
(00276 : 00300)  1 -1  1 -1 -1 -1 -1 -1 -1 -1 -1  1 -1  1 -1  1 -1 -1  1  1  1  1  1  1 -1
(00301 : 00325)  1  1  1 -1  1  1 -1  1  1 -1 -1 -1 -1  1  1 -1 -1 -1  1 -1  1 -1  1  1 -1
(00326 : 00350)  1  1  1  1  1 -1 -1 -1  1  1  1 -1 -1  1  1  1  1 -1  1 -1  1 -1  1 -1  1
(00351 : 00375)  1  1 -1 -1  1 -1 -1  1 -1  1  1  1 -1  1  1 -1 -1 -1 -1 -1  1  1 -1 -1  1
(00376 : 00400)  1  1 -1 -1  1  1  1 -1  1 -1  1 -1 -1 -1  1  1 -1  1  1  1 -1 -1  1 -1  1
(00401 : 00425) -1 -1  1 -1 -1 -1  1 -1 -1 -1  1 -1  1  1  1  1  1 -1 -1 -1 -1  1 -1  1 -1
(00426 : 00450) -1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1 -1 -1 -1  1 -1 -1  1  1  1  1 -1 -1 -1
(00451 : 00475)  1 -1  1  1 -1 -1 -1  1  1  1 -1  1 -1 -1 -1  1 -1 -1  1 -1  1  1  1 -1 -1
(00476 : 00500)  1 -1  1 -1 -1  1  1 -1  1  1  1 -1  1  1  1 -1  1  1 -1  1 -1  1 -1 -1 -1
(00501 : 00525)  1  1 -1  1  1 -1 -1 -1 -1  1  1 -1  1  1  1  1  1 -1  1 -1 -1  1 -1 -1  1
(00526 : 00550)  1  1  1  1 -1 -1  1  1 -1  1  1  1  1 -1  1 -1  1  1  1  1  1 -1 -1 -1  1
(00551 : 00575)  1 -1  1  1 -1 -1 -1  1  1  1 -1  1 -1  1  1  1 -1  1 -1  1 -1  1  1 -1 -1
(00576 : 00600)  1 -1 -1  1  1  1 -1  1  1  1  1 -1 -1  1 -1  1 -1 -1 -1 -1 -1  1 -1 -1  1
(00601 : 00625) -1 -1 -1 -1 -1 -1 -1  1 -1 -1  1 -1 -1  1 -1 -1 -1 -1 -1  1  1  1  1  1  1
(00626 : 00650)  1 -1 -1  1 -1 -1  1  1  1  1  1 -1 -1 -1  1 -1 -1 -1  1 -1 -1 -1 -1  1  1
(00651 : 00675) -1 -1 -1 -1  1 -1 -1 -1 -1 -1 -1  1 -1  1 -1  1  1 -1  1 -1  1 -1 -1  1 -1
(00676 : 00700)  1 -1 -1  1 -1  1 -1 -1 -1  1  1  1 -1  1  1 -1 -1  1  1  1 -1  1  1  1 -1
(00701 : 00725) -1 -1  1 -1 -1  1  1  1 -1 -1 -1  1  1  1 -1 -1  1  1 -1 -1  1 -1  1 -1 -1
(00726 : 00750)  1 -1 -1  1 -1  1  1 -1  1  1  1  1  1 -1  1 -1 -1 -1 -1  1  1 -1  1  1 -1
(00751 : 00775) -1 -1  1 -1  1  1 -1  1 -1 -1 -1 -1  1 -1  1  1  1 -1 -1  1  1 -1 -1  1 -1
(00776 : 00800)  1 -1  1  1 -1  1 -1  1  1  1 -1 -1 -1  1 -1  1  1  1  1 -1 -1 -1  1  1 -1
(00801 : 00825) -1  1  1  1 -1 -1  1  1 -1  1 -1  1  1 -1  1  1  1  1  1  1 -1 -1 -1  1  1
(00826 : 00850) -1 -1 -1 -1 -1  1 -1 -1 -1  1  1  1  1 -1  1  1  1  1 -1 -1  1 -1 -1  1  1
(00851 : 00875)  1  1 -1 -1  1  1 -1  1 -1  1 -1  1 -1  1  1  1 -1  1  1  1  1  1 -1 -1  1
(00876 : 00900)  1 -1  1  1  1 -1 -1 -1  1 -1  1 -1 -1  1 -1 -1 -1  1 -1  1  1 -1  1  1  1
(00901 : 00925)  1 -1 -1  1  1 -1 -1 -1 -1  1 -1  1  1  1 -1 -1  1 -1 -1 -1  1  1  1 -1 -1
(00926 : 00950)  1 -1  1 -1  1 -1 -1  1  1 -1  1 -1 -1  1  1 -1 -1  1 -1 -1 -1 -1  1  1  1
(00951 : 00975)  1 -1 -1  1  1 -1 -1  1  1  1  1 -1 -1 -1  1  1 -1 -1 -1 -1 -1 -1  1  1  1
(00976 : 01000) -1 -1 -1  1 -1 -1  1 -1 -1  1 -1 -1  1  1 -1 -1  1  1 -1  1  1 -1  1  1  1

Answer 4

I know this question was posed 5 months ago but I thought I would add something. Not sure if this will increase computing time but,

$\lambda(n)=i^{\tau(n^{2})-1}$, where $\tau(n)$ is the divisor function and $i$ is the imaginary unit. From this of course the summatory Liouville function is $L(n)=\sum_{j=1}^{n}i^{\tau(j^{2})-1}$. Perhaps the divisor function can be calculated in less time.