An Optimized Pre Computation

This post is an extension to make an optimization in precomputing the values of Totient and Mobius function. So make sure you first go through the previous post.

The extra initialization required for Algorithm 2 made me think for a while to eliminate the 'initializing' for loop. So I made (or atleast tried to make) a little optimization there. I realized that we don't need that loop. Instead, I initialized the values in the for loop only when the EulerPhi was zero. In other words, when we first reach a number, we do the initialization which should be sufficient enough. The following algorithm do that.

Algorithm 3

function PRECOMPUTATION(n)

            flag  $\gets$ [0] * (n)                                                                                            

            EulerPhi $\gets$ [0] * (n)

            Mu $\gets$ [0] * (n)

            Primes $\gets$ [0] * (n)

            pos $\gets$ 1

            for i = 2 to n do

if flag[i] = 0 then

                                    for j = i to n do

                                                flag[j] = 1

                                                if EulerPhi[j] = 0 then

                                                            EulerPhi[j] $\gets$ j

                                                            Mu[j] $\gets$ 1

                                                endif

                                                EulerPhi[j] $\gets$ ( EulerPhi[j] / i ) * ( i– 1 )

                                                if j/i%i = 0 then

                                                            Mu[j] $\gets$ 0

                                                else

                                                            Mu[j] $\gets$ Mu[j] * (-1)

                                                end if

                                                j $\gets$ j + i

                                    end for

                                    Primes[pos] $\gets$ i

                                    pos $\gets$ pos +1

                           end if

           end for

          return Primes

          return EulerPhi

          return Mu

end function

This pretty much completes our discussion of Precomputations. These algorithms can be effectively converted into codes in a programming language of your choice.

Yours Aye

Me

Pre Computation

The pre-computation of Totient function and Mobius function is one of the prerequisites of all the algorithms that we have seen so far. In fact, this pre-computation would form an excellent exercise in understand the standard Sieve of Erasthones.  In our sieving, we'll not only find the primes, but also the two aforementioned functions.

Let's first see the algorithm for Erasthones' sieve. It's not difficult understand or code. From there we'll build our Totient and Mobius functions.

Algorithm 1

function PRIMES(n)

            flag  $\gets$ [0] * (n)

            Primes $\gets$ [0] * (n)

            pos $\gets$ 1

            for i = 2 to n do

if flag[i] = 0 then

                                    for j = i to n do

                                                flag[j] $\gets$ 1
                                                j $\gets$ j + i
                                    end for

                                    Primes[pos] $\gets$ i

                                    pos $\gets$ pos +1

                           end if

            end for

            return Primes

end function

As you can see here, we start our algorithm by initializing a flag of 0's. Starting from 2, whenever we first encounter a 0, we label it as a prime. We then mark it's multiples as composites by changing their flags to 1. We then move on to find the next 0 and repeat the process.

Well, in order to bring in the totient function, we first note that the totient function of a particular number $n$ is set to its own value and changes only by its prime factorization. We encounter all the multiples of a prime in the inner for loop and use this loop to change the totient value as appropriate.

The logic applies to Mobius function as well. We initially set the value to be 1 and multiply it by -1 whenever we reach a number by a different prime. In addition, for Mobius we also have to check if its divisible by a square, and if it is, we change it to 0.

Algorithm 2

function PRECOMPUTATION(n)

            flag  $\gets$ [0] * (n)                                                                                            

            EulerPhi  $\gets$ [0, 1, 2, …, n]

            Mu  $\gets$ [1] * (n)

            Primes  $\gets$ [0] * (n)

            pos  $\gets$ 1

            for i = 2 to n do

if flag[i] = 0 then

                                    for j = i to n do

                                                flag[j]  $\gets$ 1

                                                EulerPhi[j]  $\gets$ ( EulerPhi[j] / i ) * ( i– 1 )

                                                if  j/i%i = 0 then

                                                            Mu[j]  $\gets$ 0

                                                else

                                                            Mu[j]  $\gets$ Mu[j] * (-1)

                                                end if

                                                j $\gets$ j + i

                                    end for

                                    Primes[pos]  $\gets$ i

                                    pos  $\gets$ pos +1

                           end if

          end for

          return Primes

          return EulerPhi

          return Mu

end function

This should be the end of this discussion. But I noted that it is as easy to initiate an array full of 1's (or any other value for that matter) as it is to initiate with 0. For example, when I tried to program the above algorithm in C, I first had to run the loop for the entire value of $n$, just to initialize the EulerPhi and Mobius function values. I therefore made a small optimization to the above algorithm which you can see here.

UPDATE (17/Nov/2019): Yet another interesting, and possibly faster, sieve is described in the Codeforces blog: Math note - linear sieve.

Yours Aye
Me

Liouville Summatory function

The Liouville Summatory function starts with the following well known identity

$\lambda*1=\epsilon_2$

where $\lambda$ is the Liouville's function and $\epsilon_2$ is the characteristic function of squares.

Like in the case of Merten's function, this property simplifies $\hat{F}(n)$ to a simple value.

Let $f(n)=\lambda(n)$. Then,

$F(n)=\displaystyle\sum\limits_{i=1}^n\lambda(i)=L(n)$ and $\hat{F}(n)=\displaystyle\sum\limits_{i=1}^n\sum_{d|i}\lambda(d)=\lfloor\sqrt n\rfloor$

Using these values in A special case of Dirichlet's Hyperbola method, we have

$\lfloor\sqrt n\rfloor=\displaystyle\sum\limits_{i=1}^{n/(u+1)} \left\lfloor\frac{n}{i}\right\rfloor \lambda(i) + \sum_{d=1}^u L\left(\left\lfloor\frac{n}{d}\right\rfloor\right) -u^{\text{ }}L\left(\left\lfloor\frac{n}{u+1}\right\rfloor\right)$, $u=\lfloor \sqrt{n}\rfloor$

Solving for the first term of the right summation,

$L(n)=\lfloor\sqrt n\rfloor-\displaystyle\sum\limits_{i=1}^{n/(u+1)} \left\lfloor\frac{n}{i}\right\rfloor \lambda(i) - \sum_{d=2}^u L\left(\left\lfloor\frac{n}{d}\right\rfloor\right) +u^{\text{ }}L\left(\left\lfloor\frac{n}{u+1}\right\rfloor\right)$, $u=\lfloor \sqrt{n}\rfloor$

We can use the intermediate result we obtained in Dirichlet's hyperbola method to write

$\lfloor\sqrt n\rfloor=\displaystyle\sum\limits_{k=1}^nL\left(\left\lfloor\frac{n}{k}\right\rfloor\right)$

$L(n)=\lfloor\sqrt n\rfloor-\displaystyle\sum\limits_{k=2}^nL\left(\left\lfloor\frac{n}{k}\right\rfloor\right)$


Yours Aye
Me

Moment of a function

Define the $m$th moment of a function $f(n)$ as

$F_{(m)}(n)=\displaystyle\sum\limits_{k=1}^n k^mf(k)$

Now the same procedure also relates the moments of the three functions, since

$f(n)*g(n)=h(n)\implies n^mf(n)*n^mg(n)=n^mh(n)$

For example, since we know that $\varphi(n)*1=n$, we can use this to calculate the Summatory Totient moments.

Choosing $f(n)=\varphi(n)$ and $g(n)=1$, we have $n^k\varphi(n)*n^k=n^{k+1}$. The corresponding summation functions are

$F_k(n)=\displaystyle\sum\limits_{m=1}^n m^k\varphi(m)$, $G(n)=S_k(n)$ and $H(n)=S_{k+1}(n)$. The corresponding results are

$F_k(n)=S_{k+1}(n)-\displaystyle\sum\limits_{m=2}^n m^k F_k\left(\left\lfloor\frac{n}{m}\right\rfloor\right)$

$F_k(n)=S_{k+1}(n)-\displaystyle\sum\limits_{m=1}^{n/(u+1)}m^k\varphi(m)S_k\left(\left\lfloor\frac{n}{m}\right\rfloor\right)-\sum_{m=2}^u m^kF_k\left(\left\lfloor\frac{n}{m}\right\rfloor\right)+S_k(u)F_k\left(\left\lfloor\frac{n}{u+1}\right\rfloor\right)$


Yours Aye
Me

Jordan Summatory function

Now that we have seen Dirichlet's hyperbola method, we try to use it find sub-linear algorithms for different functions.

Jordan's totient function, denoted as $J_k(n)$, generalizes Euler's totient function. Let's use $\mathbf{J}_k(n)$ to denote Jordan summatory function. That is

$\mathbf{J}_k(n)=\displaystyle\sum\limits_{m=1}^nJ_k(m)$

We know that,$1*J_k(n)=n^k$, where '$*$' denotes Dirichlet convolution.

Choosing $f(n)=J_k(n)$ and $g(n)=1$, we have $h(n)=f*g=n^k$. The corresponding summatory functions are

$F(n)=\mathbf{J}_k(n)$, $G(n)=n$ and $H(n)=S_k(n)$.

Using the results obtained in Dirichlet's hyperbola method, we get

$\mathbf{J}_k(n)=S_k(n)-\displaystyle\sum\limits_{m=2}^n\mathbf{J}_k\left(\left\lfloor\frac{n}{m}\right\rfloor\right)$

$\mathbf{J}_k(n)=S_k(n)-\displaystyle\sum\limits_{m=1}^{n/(u+1)}  \left\lfloor\frac{n}{m}\right\rfloor J_k(n) - \sum_{m=2}^u \mathbf{J}_k\left(\left\lfloor\frac{n}{m}\right\rfloor\right) +u^\text{ }\mathbf{J}_k\left(\left\lfloor\frac{n}{u+1}\right\rfloor\right)$


Yours Aye
Me

Dirichlet's Hyperbola method

Let $f(n)$ and $g(n)$ be two arithmetic functions. Define $h(n)=(f*g)(n)$ as the Dirichlet convolution of $f$ and $g$. That is,

$h(n)=\displaystyle\sum\limits_{d|n}f(d)g\left(\frac{n}{d}\right)$

We now define the summatory functions $F(n)$, $G(n)$ and $H(n)$.

$F(n)=\displaystyle\sum\limits_{k=1}^nf(k)$ and likewise for the other two functions.

We'll now see how Dirichlet's hyperbola method relates the three summatory functions. If we expand $H(n)$ and collect the $f$ terms (or the $g$ terms) and use the definition of $G(n)$ (or that of $F(n)$), we'll get

$H(n)=\displaystyle\sum\limits_{k=1}^nf(k)G\left(\left\lfloor\frac{n}{k}\right\rfloor\right)=\displaystyle\sum\limits_{k=1}^ng(k)F\left(\left\lfloor\frac{n}{k}\right\rfloor\right)$

This by itself is a powerful result. We'll see later how this serves to give reccurence relations for calculating summatory functions. Now we know that the floor function remains constant over a long range and it is something to take advantage of. Using the same techniques that were used in A special case of Dirichlet's hyperbola method, we can write

$H(n)=\displaystyle\sum\limits_{k=1}^{n/(u+1)}  f(k)G\left(\left\lfloor\frac{n}{k}\right\rfloor\right) + \sum_{k=1}^u g(k)F\left(\left\lfloor\frac{n}{k}\right\rfloor\right) -G(u)F\left(\left\lfloor\frac{n}{u+1}\right\rfloor\right)$

or solving for $F(n)$,

$g(1)F(n)=H(n)- \displaystyle\sum\limits_{k=2}^u g(k)F\left(\left\lfloor\frac{n}{k}\right\rfloor\right) -\displaystyle\sum\limits_{k=1}^{n/(u+1)}  f(k)G\left(\left\lfloor\frac{n}{k}\right\rfloor\right) + G(u)F\left(\left\lfloor\frac{n}{u+1}\right\rfloor\right)$

where $u=\lfloor \sqrt{n}\rfloor$. This method can be used to obtain sub-linear algorithms for Totient summatory function, Mertens function, Liouville Summatory function and other functions.

With a similar procedure, we can create many more just by knowing the Dirichlet convolution between the two functions. Though am not familiar with the analysis of algorithms, I think the first formula in each case is an $O(n^{\frac{3}{4}})$ algorithm and the second one is a $O(n^{\frac{2}{3}})$ algorithm.

UPDATE (17/Nov/2019): An interesting and similar algorithm is described in this Codeforces blog: Looking for Extended Eratosthenes sieve tutorial.

Yours Aye'
Me

Mertens Function

Mertens function can be informally called as the Moebius Summatory function. Well because,

$M(n)=\displaystyle\sum\limits_{k=1}^n\mu(n)$

where $\mu(n)$ is the Mobius function.

Very similar to computation of Totient summatory function is the idea of Merten's function. Like in the case of Totient summatory function, applying A special case of Dirichlet's Hyerbola method  exploits a property of Mobius function, $\mu(n)$. For integer $n$, we have

$\displaystyle\sum\limits_{d|n}\mu(n)=
\begin{cases}
1,&\text{if }n=1\\
0,&\text{if }n>1
\end{cases}$

In other words, $\mu*1=\epsilon$, where $\epsilon$ is the multiplicative identity. (i.e. $\epsilon(1)=1$, all other values $0$).

What this property does is simplify $\hat{F}(n)$ to a incredibly simple value. Let $f(n)=\mu(n)$. Then,

$F(n)=\displaystyle\sum\limits_{i=1}^n\mu(i)=M(n)$ and $\hat{F}(n)=\displaystyle\sum\limits_{i=1}^n\sum_{d|i}\mu(d)=1$

Using these values in A special case of Dirichlet's Hyperbola method, we have

$1=\displaystyle\sum\limits_{i=1}^{n/(u+1)} \left\lfloor\frac{n}{i}\right\rfloor \mu(i) + \sum_{d=1}^u M\left(\left\lfloor\frac{n}{d}\right\rfloor\right) -u^{\text{ }}M\left(\left\lfloor\frac{n}{u+1}\right\rfloor\right)$, $u=\lfloor \sqrt{n}\rfloor$

Solving for the first term of the right summation,

$M(n)=1-\displaystyle\sum\limits_{i=1}^{n/(u+1)} \left\lfloor\frac{n}{i}\right\rfloor \mu(i) - \sum_{d=2}^u M\left(\left\lfloor\frac{n}{d}\right\rfloor\right) +u^{\text{ }}M\left(\left\lfloor\frac{n}{u+1}\right\rfloor\right)$, $u=\lfloor \sqrt{n}\rfloor$

Precomputing the values of $\mu(k)$ for $k \leq \sqrt{n}$ and memoizing, gives a clean method to calculate the Mertens function. Using a similar code, I can computer $M(10^9)=-222$ in $70$ seconds, $M(10^8)=1928$ in $15$ seconds and $M(10^7)=1037$ in $2$ seconds.

Again as in Totient summatory function, we can use the intermediate result we obtained in Dirichlet's hyperbola method to write

$1=\displaystyle\sum\limits_{k=1}^nM\left(\left\lfloor\frac{n}{k}\right\rfloor\right)$

$M(n)=1-\displaystyle\sum\limits_{k=n}^nM\left(\left\lfloor\frac{n}{k}\right\rfloor\right)$

Again this method avoids precomputation at cost of being time-expensive.


Yours Aye
Me