Sunday, November 29, 2015

C-code for the PrimePi function

This code uses the algorithm described here. The same algorithm has been used for SumOfPrimes function here.


#include<stdio.h>
#include<math.h>

#define M 1000000
#define Mod 1000000000

int Lower[M+1], Higher[M+1];
char Used[M+1];

long long SumNat(long long n)
{
    long long p = n, q = n+1;
    if(p%2)
        q /= 2;
    else
        p /= 2;
    return (p % Mod) * (q % Mod) % Mod;
}

long long PrimePi(long long n)
{
    long long v = sqrt(n), p, temp, q, j, end, i, d, t;
    Higher[1] = n % Mod - 1;
    for(p = 2; p <= v; p++) {
        Lower[p] = p - 1;
        Higher[p] = (n/p) % Mod - 1;
    }
    for(p = 2; p <= v; p++) {
        if(Lower[p] == Lower[p-1])
            continue;
        temp = Lower[p-1];
        q = p * p;
        Higher[1] -= Higher[p] - temp;
        if(Higher[1] < 0)
            Higher[1] += Mod;
        j = 1 + (p & 1);
        end = (v <= n/q) ? v : n/q;
        for(i = p + j; i <= 1 + end; i += j) {
            if(Used[i])
                continue;
            d = i * p;
            if(d <= v)
                Higher[i] -= Higher[d] - temp;
            else {
                t = n/d;
                Higher[i] -= Lower[t] - temp;
            }
            if(Higher[i] < 0)
                Higher[i] += Mod;
        }
        if(q <= v)
            for(i = q; i <= end; i += p*j)
                Used[i] = 1;
        for(i = v; i >= q; i--) {
            t = i/p;
            Lower[i] -= Lower[t] - temp;
        }
    }
    return Higher[1] % Mod;
}

int main() {
    printf("%lld", PrimePi(1000000000000LL));
    return 0;
}

Well, I conclude the discussion on the PrimePi functions with this post, I'll think of something to continue blogging which may take a while.



Until then
Yours Aye
Me

C-code for the SumOfPrimes function


We saw the algorithm for the generalized Prime counting function in the previous post. Here is the C code for the same, with some slight optimizations.


#include<stdio.h>
#include<math.h>

#define M 1000000
#define Mod 1000000000

long long Lower[M+1], Higher[M+1];
char Used[M+1];

long long SumNat(long long n)
{
    long long p = n, q = n+1;
    if(p%2)
        q /= 2;
    else
        p /= 2;
    return (p % Mod) * (q % Mod) % Mod;
}

long long SumPrimePi(long long n)
{
    long long v = sqrt(n), p, temp, q, j, end, i, d, t;
    Higher[1] = SumNat(n) - 1;
    for(p = 2; p <= v; p++) {
        Lower[p] = SumNat(p) - 1;
        Higher[p] = SumNat(n/p) - 1;
    }
    for(p = 2; p <= v; p++) {
        if(Lower[p] == Lower[p-1])
            continue;
        temp = Lower[p-1];
        q = p * p;
        Higher[1] -= (Higher[p] - temp) % Mod * p % Mod;
        if(Higher[1] < 0)
            Higher[1] += Mod;
        j = 1 + (p & 1);
        end = (v <= n/q) ? v : n/q;
        for(i = p + j; i <= 1 + end; i += j) {
            if(Used[i])
                continue;
            d = i * p;
            if(d <= v)
                Higher[i] -= (Higher[d] - temp) % Mod * p % Mod;
            else {
                t = n/d;
                Higher[i] -= (Lower[t] - temp) % Mod * p % Mod;
            }
            if(Higher[i] < 0)
                Higher[i] += Mod;
        }
        if(q <= v)
            for(i = q; i <= end; i += p*j)
                Used[i] = 1;
        for(i = v; i >= q; i--) {
            t = i/p;
            Lower[i] -= (Lower[t] - temp) % Mod * p % Mod;
            if(Lower[i] < 0)
                Lower[i] += Mod;
        }
    }
    return Higher[1] % Mod;
}

int main() {
    printf("%lld", SumPrimePi(1000000000000LL));
    return 0;
}


The code computes the value of $\pi_1(10^{12})$ in about 6 seconds which is way, way better than Mathematica's computation time. If you are looking for the C-code for $\pi(x)$, you can find it in this post.

Saturday, November 28, 2015

Algorithm for the Prime counting function


In Counting Primes - Prime Counting function - Part III, we saw the function used to calculate the sum of all primes less than a given value and it's Mathematica code. Well, the Mathematica code is good enough but we can do better with C. But if we try to convert the code given into C, we run into a problem. In the given code, we have relied heavily upon Mathematica's in-built functions. But we can sort the problem easily by using the some of the properties $\varphi_k(x,\pi(\sqrt{x}))$. First note that,

Result 1 : $\pi_k(x)=\varphi_k(x,\pi(\sqrt{x}))=\varphi_k(x,\sqrt{x})=\varphi_k(x,x-2)=\varphi_k(x,x-1)$, since $\pi(x) \leq x$.

Secondly, $\pi_k(x)$ increases only at prime values of $x$. This can be used to test the primality of $x$ in our C-code. Or casting it in different terms, we can say

Result 2 : $x$ is prime, if and only if $\pi_k(x) > \pi_k(x-1)$ (or) $\varphi_k(x,x-2) > \varphi_k(x-1,x-2)$.



The Algorithm to find $\pi_k(n)$ will now proceed as follows.

Input : Get the value of $n$.

Step 1: Create a table with 2 rows. Mark the first row as $m$ and fill it as $1,2,\cdots,\lfloor\sqrt{n}\rfloor, \lfloor n/\lfloor\sqrt{n}\rfloor\rfloor,\cdots,\lfloor n/2\rfloor,\lfloor n/1 \rfloor$. Mark the second row as $\varphi_k(m)$ and the fill the corresponding $\varphi_k(m,0)$ values.
Step 2 : Initiate $r=2$.
Step 3 : Check whether $r$ prime or not. To do this, compare the corresponding $\varphi$ entries of $r$ and $r-1$ (making use of Result 2).
Step 4 : If they are same, skip to Step 6. Else, conclude $r$ to be a prime.
Step 5 : For each $m\geq r^2$ , starting from the last entry and moving left, update the values of $\varphi_1$ using the recurrence, $\varphi_k(m)=\varphi_k(m)-r^k*(\varphi_k(\lfloor m/r \rfloor)-\varphi_k(r-1))$
Step 6 : Increment $r$.
Step 7: If $r$ equals $\lfloor\sqrt{n}\rfloor$, output the rightmost entry of the second row and Stop. Else, goto Step 3.



For example, the tables while computing $\pi_1(20)$ are given below.

TABLE 1 :: $\varphi_1(m,0)$
$
\begin{matrix}
m&1&2&3&4&5(=\lfloor20/4\rfloor)&6(=\lfloor20/3\rfloor)&10(=\lfloor20/2\rfloor)&20(=\lfloor20/1\rfloor)\\
\varphi_1(m)&0&2&5&9&14&20&54&209\\
\end{matrix}
$

TABLE 2 :: $\varphi_1(m,1)$
$
\begin{matrix}
m&1&2&3&4&5&6&10&20\\
\varphi_1(m)&0&2&5&5&10&10&26&101\\
\end{matrix}
$

TABLE 3 :: $\varphi_1(m,2)$
$
\begin{matrix}
m&1&2&3&4&5&6&10&20\\
\varphi_1(m)&0&2&5&5&10&10&17&77\\
\end{matrix}
$

Note a couple of things here. The table is modified $\pi(\lfloor\sqrt{n}\rfloor)$ times after the initialization. As a byproduct of the algorithm, we get the SumOfPrimes for all the entries in the first row. Also, we generate the list of primes $\leq \lfloor\sqrt{n}\rfloor$.

After I wrote this, I found that no text can really equal a cute demonstration in terms of clarity. So I found this, which will visualizes the algorithm's execution. Hope that clarifies things. I'll post the C-code for this algorithm in the next post.


Until then
Yours Aye
Me

Friday, November 27, 2015

Pólya's Enumeration Theorem with Applications - Part II

This post continues what we discussed here.

Pólya's Enumeration theorem allows us to do much more. In all of the examples above, we did not care about the number of balls we had for each color. In fact, that information was not needed to answer the questions we asked ourselves. Also, in both the necklace and the urn problems, we allowed only one ball at a given position. What if there can be more than one ball at one particular position? Consider the following question.

Example 5: 8 indistinguishable Urns arranged in a circle. What is the number of ways to distribute $n$ blue balls into these urns discounting rotational symmetry?

To make the question clear, number the urns from 1 to 8 in clockwise direction. Let's say we have 16 blue balls. Then according to the question, the distributions {1,2,3,4,5,1,0,0} and {5,1,0,0,1,2,3,4} are identical because one can be transformed to other by rotations. To solve this problem with Pólya's Enumeration Theorem we first find the cycle index of Cyclic group of degree 8.

$Z_{C_8}(x_1,x_2,...,x_8)=\frac{1}{8}(x_1^8+x_2^4+2x_4^2+4x_8)$

Pólya's Enumeration theorem says the answer we are looking for is given by

$\frac{1}{8}((X_1^8)_n+(X_2^4)_n+2(X_4^2)_n+4(X_8)_n)$

where $X_a^b=\frac{1}{(1-z^a)^b}$ and $(F(z))_n$ is the coefficient of $z^n$ in the power series of $F(z)$. With Mathematica's help we find the answer to be $30667$ when $n=16$.

The same procedure can be applied for all the permutations groups and each will give the required answer according to their nature. One thing of interest here is that when there are $n$ balls of the same color, $Z_{S_r}$ counts the number of partitions of $n$ into $r$ parts ($0$ included) and $Z_{-S_r}$ counts the number of partitions of $n$ into $r$ distinct parts ($0$ included). This interpretation has an interesting application in the multiplicative partition of $n$ which will be final piece of our discussion.

Now, the previous problem had balls of the same color. What if there are $n_i$ balls of color $c_i$. This is the second main point that I wished to convey in this post. Let's look at the next problem.

Example 6: 4 indistinguishable Urns are arranged in a circle. What is the number of ways to distribute $m$ different colored balls, with $n_i$ balls of color $c_i$, into these urns discounting rotational symmetry?

The answer to this question proceeds in pretty much the same way as the previous answer. We find the cycle index of Cyclic group of degree 4. The important distinction comes in the definition of $(F(z))_n$. Write $n=(n_1,n_2,..,n_m)$ and define,

$(F(z))_n=\Pi_{j=1}^m(F(z))_{n_j}$ where $F(z)$ will be of the form $\Pi_{i=1}^rX_i^\alpha$.

In case of 4 Urns with 1 white balls and 1 black balls discounting rotational symmetry, we use

$Z_{C_4}(x_1,x_2,x_3,x_4)=\frac{1}{4}(x_1^4+x_2^2+2x_4)$

and the final answer will be $\frac{1}{4}((X_1^4)_{(1,1)}+(X_2^2)_{(1,1)}+2(X_4)_{(1,1)}) = 4$

because {bg,$0$,$0$,$0$}, {b,$0$,g,$0$}, {b,g,$0$,$0$} and {g,b,$0$,$0$} are the possible arrangements.

Example 7: What is the number of ways to distribute 2 white balls and 1 black balls into 3 indistinguishable Urns without taking order into account?

This question can be interpreted in at least two (equivalent) ways. First, it is the number of ways of partitioning the 2-partite number (2,1) into three 2-partite numbers. Second, it is the number of multiplicative partitions of $12(=2^2\times3^1)$. That's easy. $12 = 1\times 2\times 6 = 1\times 3\times 4 = 2\times 2\times 3 = 1\times 1\times 12$. Let's check whether Pólya's Enumeration Theorem gives us the same answer.

$Z_{S_3}=\frac{1}{6}(x_1^3+3x_1x_2+2x_3)$

Following the same procedure as above we see the answer is $\frac{1}{6}((X_1^3)_{(2,1)}+3(X_1X_2)_{(2,1)}+2(X_3)_{(2,1)})$.  Let's just evaluate one of these terms, preferably $(X_1X_2)_{(2,1)}$.

$(X_1X_2)_{(2,1)}=[z^2]\frac{1}{(1-z)(1-z^2)}\times[z^1]\frac{1}{(1-z)(1-z^2)}=2\times1=2$.

The final answer is $\frac{1}{6}(6\times3+3\times2\times1+2\times0\times0)=4$ which we saw already. In terms of balls and Urns the distribution would look like {$0$,w,wb}, {$0$,b,ww}, {w,w,b} and {$0$,$0$,wwb}.

Example 8: What is the answer to the same question if we want no two urns to have the same pattern of balls?

In this case it asks for the multiplicative partition of 12 with distinct numbers in which case the last two multiplicative partitions listed above will be eliminated. In terms of Urns and balls, {w,w,b} will be eliminated because two urns have one white ball each and {$0$,$0$,wwb} will be eliminated because two urns are empty. To get the answer via Pólya's Enumeration theorem, we'll be using $Z_{-S_3}$ instead of $Z_{S_3}$.

$Z_{-S_3}=\frac{1}{6}(x_1^3-3x_1x_2+2x_3)$

Hence the answer will be  $\frac{1}{6}((X_1^3)_{(2,1)}-3(X_1X_2)_{(2,1)}+2(X_3)_{(2,1)})=2$.

Well, that's it from my end. I'm quite impressed with the way I progressed but it's still upto a first time reader to really judge how good I was in explaining what I intended to explain. Any suggestions and corrections are welcome. It feels extremely good to have come up with this post and I enjoyed it. Hope you'll do too. C ya later.

Until then
Yours Aye,
Me

Sunday, September 13, 2015

Counting primes - Prime Counting function - Part III


In the lines of Divisor Summatory function, let us define a Prime Summatory function $\pi_k(x)$. That is,

$\pi_k(x)=\displaystyle\sum_{p\leq x}p^k$

where $p$ is prime.

Then we can use the same function we used in Counting primes - Prime Counting function - Part II to calculate $\pi_k(x)$. We can define a function $\varphi_k(x,a)$ which gives the sum of the $k$th powers of integers between $2$ and $x$ (both inclusive), that remains unmarked after sieving with the first $a$ primes. Then as before,

$\pi_k(x) = \varphi_k(x,\pi(\sqrt{x}))$

and

$\varphi_k(x,a)=
\begin{cases}
\varphi_k(x,a-1)-p_a^k(\varphi_k(x/p_a,a-1)-\varphi_k(p_a-1,a-1)),&\text{if }x\geq p_a^2\\
\varphi_k(x,\pi(\sqrt{x})),&\text{otherwise}
\end{cases}$

such that $\varphi_k(0,a)=\varphi_k(1,a)=0$ and $\varphi_k(x,0)=2^k + 3^k + \cdots + x^k$.

Even the same Mathematica code used in the previous post with a little modification does the job. For example, the following Mathematica code gives the sum of primes less than a given value.

Clear[S, PrimeP, Su]
S[n_, m_] :=
 S[n, m] =
  Which[n <= 1, 0, m == 0, n - 1, n < Power[Prime[m], 2], PrimeP[n],
   True, S[n,
     m - 1] - (S[Floor[n/Prime[m]], m - 1] - S[Prime[m] - 1, m - 1])]
PrimeP[n_] :=
 PrimeP[n] = If[n <= 3, n - 1, S[n, PrimeP[Floor[Sqrt[n]]]]]
Su[n_, m_] :=
 Su[n, m] =
  Which[n <= 1, 0, m == 0, n (n + 1)/2 - 1, n < Power[Prime[m], 2],
   Su[n, PrimeP[Floor[Sqrt[n]]]], True,
   Su[n, m - 1] -
    Prime[m] (Su[Floor[n/Prime[m]], m - 1] - Su[Prime[m] - 1, m - 1])]
n = Power[10, 9];
Timing[Block[{RecursionLimit = Power[10, 8]}, Su[n, PrimeP[n]]]]

This code computes the value of $\pi_1(10^8)$ in $7$ seconds and the value of $\pi_1(10^9)$ in about $30$ seconds. Though the code does not have any significant changes from that of $\pi(x)$, the code $\pi_1(x)$ still relies on $\pi(x)$ which increases the calculation time.

Won't it be nice if we can somehow eliminate that and thereby increase the performance? We'll see about that in the next post.


Until then
Yours Aye
Me

Saturday, September 12, 2015

Counting primes - Prime Counting function - Part II


We discussed the Legendre's formula associated with the prime counting function in the previous post. This post discusses an alternative formulation whose importance will be apparent when we discuss the SumOfPrimes function.

Define a new function $\varphi(x,a)$ such that it counts the number of integers between $2$ and $x$ that remain unmarked after sieving with the first $a$ primes. Such a function will be related to the prime-counting function as

$\pi(x) = \varphi(x,\pi(\sqrt{x}))$

This function satisfies the recurrence relation

$\varphi(x,a)=
\begin{cases}
\varphi(x,a-1)-(\varphi(x/p_a,a-1)-\varphi(p_a-1,a-1)),&\text{if }x\geq p_a^2\\
\pi(x),&\text{otherwise}
\end{cases}$

such that $\varphi(0,a)=\varphi(1,a)=0$ and $\varphi(x,0)=x-1$.

Note that whereas $\phi(10,2)=3$ (the three numbers being $1, 5$ and $7$), $\varphi(10,2)=4$ (the four numbers are $2, 3, 5$ and $7$). Because when sieve with the $2$, every even number except $2$ goes marked, with $3$ every mulitple of $3$ except $3$ itself goes marked and so on.

The following Mathematica code uses these properties of $\varphi(x,a)$ to calculate $\pi(x)$.

Clear[S, PrimeP, Su]
S[n_, m_] :=
 S[n, m] =
  Which[n <= 1, 0, m == 0, n - 1, n < Power[Prime[m], 2], PrimeP[n],
   True, S[n,
     m - 1] - (S[Floor[n/Prime[m]], m - 1] - S[Prime[m] - 1, m - 1])]
PrimeP[n_] :=
 PrimeP[n] = If[n <= 3, n - 1, S[n, PrimeP[Floor[Sqrt[n]]]]]
Su[n_, m_] :=
 Su[n, m] =
  Which[n <= 1, 0, m == 0, n (n + 1)/2 - 1, n < Power[Prime[m], 2],
   Su[n, PrimeP[Floor[Sqrt[n]]]], True,
   Su[n, m - 1] -
    Prime[m] (Su[Floor[n/Prime[m]], m - 1] - Su[Prime[m] - 1, m - 1])]
n = Power[10, 9];
Timing[Block[{RecursionLimit = Power[10, 8]}, PrimeP[n]]]

This code computes $\pi(10^7)$ in a sec, takes approx. 4 secs for $\pi(10^8)$ and 18 secs for $\pi(10^9)$. As we can see, there is no significant increase in performance between the two codes. But as I said at the beginning of this post, this function can be easily tweaked to find the SumOfPrimes function because of its lack of dependence on $a$ in the second case of the recurrence.

We'll see about that in detail in the next post.


Until then
Yours Aye
Me

Counting primes - The Prime-counting function - Legendre's Phi function


The prime-counting function, denoted as $\pi(x)$, counts the number of primes less than or equal to some value x. Though there are many formulas to calculate the value of this function approximately, finding the exact value is a challenging problem, especially for beginner level programmer.

Apart from naive counting, the next best formula used for counting primes was given by Legendre. You can read about it in the wiki page. This formula introduces a new function called the Legendre's phi function, which is then used to compute $\pi(x)$.

However, in this post we'll use a closely related function that is slightly more efficient than the Lengendre's formula given in Wiki or Math-world page.

According to Math-world, Legendre's Phi function, denoted as $\phi(x,a)$, counts the number of integers less or equal to $x$ which are not divisible by any of the first $a$ primes. This function is then related to the prime-counting function as

$\phi(x,\pi(\sqrt{x}))=\pi(x)-\pi(\sqrt{x})+1$

Also, it follows the following recurrence relation

$\phi(x,a)=
\begin{cases}
\phi(x, a-1) - \phi(x/p_a, a-1), &\text{if }x\geq p_a^2 \\
\pi(x)-a+1,&\text{otherwise}
\end{cases}$

where $p_a$ denotes the $a$th prime. Also, $\phi(0, a) = 0$, $\phi(1, a) = 1$ and $\phi(x, 0)=x$.

The following Mathematica code uses these relations to compute $\pi(x)$.

Clear[S, PrimeP]
S[n_, m_] :=
 S[n, m] =
  Which[n <= 1, n, m == 0, n, n < Power[Prime[m], 2],
   PrimeP[n] - m + 1, True, S[n, m - 1] - S[Floor[n/Prime[m]], m - 1]]
PrimeP[n_] :=
 PrimeP[n] =
  If[n <= 3, n - 1,
   S[n, PrimeP[Floor[Sqrt[n]]]] + PrimeP[Floor[Sqrt[n]]] - 1]
n = Power[10, 9];
Timing[Block[{RecursionLimit = Power[10, 8]}, PrimeP[n]]]

The code takes approx a second to calculate $\pi(10^7)$, 4 secs for $\pi(10^8)$ and 20 secs for $\pi(10^9)$. Though these timings are reasonable for Mathematica, and will speed up considerably if we do in C, we'll see an alternative formulation in the next post.

Until then
Yours Aye
Me

Friday, May 29, 2015

C-Code for Merten's Function


The C-Code for Merten's function proceeds much in the same way as that of Totient Summatory function and uses the formula derived in Merten's function. This code evaluates the value of $M(10^{11})$ in approx. $8$ seconds and $M(10^{12})$ in approx. $80$ seconds. Here is the code.

#include<stdio.h>
#include<math.h>

#define M 60000000
#define Mod 1000000007

int MertCache[M+1], Mert[M+1];
char Mu[M+1], flag[M+1];

void precalc()
{
    int i, j;
    Mu[1] = 1;
    Mert[1] = 1;
    MertCache[1] = 56789;
    for(i = 2; i <= M; i++) {
        if(!flag[i])
            for(j = i; j <= M; j += i) {
                if( !flag[j] )
                    Mu[j] = 1;
                if(j/i%i == 0)
                    Mu[j] = 0;
                else
                    Mu[j] = -Mu[j];
                flag[j] = 1;
        }
        Mert[i] = Mert[i-1] + Mu[i];
        MertCache[i] = 56789;
    }
}

long long Merten(long long n, long long k)
{
    long long m = n/k;
    if(m <= M)
        return Mert[m];
    else if(MertCache[k] != 56789)
        return MertCache[k];
    long long p = sqrt(m), q = m/(p+1), i, res = 0;
    for(i = 2; i <= p; i++)
        res += Merten(n, k*i);
    for(i = 1; i <= q; i++)
        res = res + (m/i) * Mu[i];
    res = 1 - res;
    res += p * Merten(n, k*(p+1));
    return MertCache[k] = res;
}

int main()
{
    precalc();
    printf("%lld", Merten(100000000000LL,1));
    return 0;
}

Note: $\text{Merten}(n,k) = M\left(\left\lfloor\dfrac{n}{k}\right\rfloor\right)$ in the code.

This code doesn't include the Phi function at all. But if you plan to use both the Totient summatory function and the Merten's function in the same code, I think code may slow down. I'vnt tried it yet but I'll also post the code that uses both the functions.

UPDATE (13 May 2017): As in Totient Summatory function, I made an iterative version of the Merten's function. You can check the performance of the code in this Online IDE.

from timeit import default_timer

start = default_timer()

n = 10 ** 6
s = int(n ** 0.5)


def sieve():
    Lows, Flag = [0] + [1] * s, [0] * (1 + s)
    Lows[1] = 1
    for i in range(2, 1 + s):
        if Flag[i] == 0:
            for j in range(i, 1 + s, i):
                Flag[j] = 1
                if (j // i % i) == 0:
                    Lows[j] = 0
                else:
                    Lows[j] *= -1
        Lows[i] += Lows[i - 1]
    return Lows


Lows = sieve()
Highs = [0] * (1 + s)

for k in range(s, 0, -1):
    tempn = n // k
    u = int(tempn ** 0.5)
    v = tempn // (1 + u)
    res = 1
    for i in range(1, 1 + v):
        res -= (Lows[i] - Lows[i - 1]) * (tempn // i)
    for i in range(2, 1 + u):
        res -= Highs[k * i] if k * i <= s else Lows[tempn // i]
    i = 1 + u
    res += Lows[tempn // i] * u
    Highs[k] = res

print(Highs[1])
print(default_timer() - start)

Note: $\text{Highs}[k] = M\left(\left\lfloor\dfrac{n}{k}\right\rfloor\right)$ in the code. Also, the value $-962982$ is purely arbitrary just to make sure that value has not been touched yet.


Until then
Yours Aye
Me

C Code for Totient Summatory function


I thought of creating a post on Dirichlet Generating functions which would be right continuation of the direction my blog is heading. But this post is about the C code for Totient Summatory function.

The code uses the algorithm described in Totient Summatory function. If you haven't visited that post, I suggest you to do so. The code computes the answer for $10^{11}$ in approx. $7$ seconds and for $10^{12}$ in approx. $25$ seconds. Not bad compared to Mathematica. Here's the code for you.

#include<stdio.h>
#include<math.h>

#define M 60000000
#define Mod 1000000007

int Phi[M+1], SumPhi[M+1], TotCache[M+1];

void precalc()
{
    int i, j;
    Phi[1] = 1;
    SumPhi[1] = 1;
    for(i = 2; i <= M; i++) {
        if(!Phi[i])
            for(j = i; j <= M; j += i) {
                if(!Phi[j])
                    Phi[j] = j;
                Phi[j] = (Phi[j] / i) * (i - 1); }
        SumPhi[i] = (SumPhi[i-1] + Phi[i]) % Mod; }
}

long long SumNat(long long n)
{
    long long p = n, q = n+1;
    if(p%2)
        q /= 2;
    else
        p /= 2;
    return (p % Mod) * (q % Mod) % Mod;
}

long long TotientSum(long long n, long long k)
{
    long long m = n / k;
    if(m <= M)
        return SumPhi[m];
    else if(TotCache[k] > 0)
        return TotCache[k];
    long long p = sqrt(m), q = m / (p+1), i, res = 0;
    for(i = 2; i <= p; i++) {
        res += TotientSum(n, k * i);
        if(res >= Mod)
            res %= Mod;
    }
    for(i = 1; i <= q; i++) {
        res = res + ((m / i) % Mod * Phi[i]) % Mod;
        if(res >= Mod)
            res %= Mod;
    }
    res = SumNat(m) - res;
    res += (p % Mod) * TotientSum(n, k * (p + 1));
    res %= Mod;
    TotCache[k] = res;
    return res;
}

int main()
{
    precalc();
    printf("%lld",TotientSum(1000000000000LL,1));
    return 0;
}

Note: $\text{TotientSum}(n,k)=\Phi\left(\left\lfloor\dfrac{n}{k}\right\rfloor\right)$ in the code.

Just posting a code wouldn't do much good either for the blogger or the reader I think from the experience I've had myself. But am not gonna no through the complete code. Just a few details which we should be paying attention to.

First, the 'precalc' function. It's just a word-by-word translation of the algorithm that we saw on Optimized Precomputation.

Second is the 'SumNat' function. The important thing about this function is that we ought to keep the values from overflowing. For example, I first did $\text{(p%Mod * q%Mod)%Mod}$ which was perfectly okay looking except that it overflows.

Third, the 'TotientSum' function, is pure memoization. It'll take a separate post to completely explain it. Maybe I'll do that later. But for now I think this should suffice. I'll try to come up with the C-code for Merten's function in the next post.

Last is the limit you choose to sieve. That's the real trade-off you have to make between space and time. If you choose to sieve more, it's gonna take more space but the final result will be quick, whereas if you choose to memoize more you consume less space at cost of being time expensive. I made some trial and error to arrive at the above sieve limit which kinda minimized the computation time.

UPDATE(13 Mar 2017): I made an iterative version of the Totient summartory function in Python. You can check it out in this Online IDE.

from timeit import default_timer

start = default_timer()

n = 1000000000
s = int(n ** 0.5)
LowPhis = [k for k in range(1 + s)]
HiPhis = [0] * (1 + s)

for p in range(2, 1 + s):
    if p == LowPhis[p]:
        for k in range(p, 1 + s, p):
            LowPhis[k] = LowPhis[k] - LowPhis[k] // p
    LowPhis[p] += LowPhis[p - 1]

for k in range(s, 0, -1):
    tempn = n // k
    res = tempn * (1 + tempn) // 2
    u = int(tempn ** 0.5)
    v = tempn // (1 + u)
    for i in range(1, 1 + v):
        res -= (tempn // i) * (LowPhis[i] - LowPhis[i - 1])
    for i in range(2, 1 + u):
        res -= (HiPhis[i * k] if i * k <= s else LowPhis[tempn // i])
    i = 1 + u
    res += u * LowPhis[tempn // i]
    HiPhis[k] = res

print(HiPhis[1])
print(default_timer() - start)

Note: $\text{HiPhis}[k]=\Phi\left(\left\lfloor\dfrac{n}{k}\right\rfloor\right)$ in the code.

I like iterative versions in that they don't suffer from the annoying maximum-recursion-depth error. I think this is a little slower than the recursive version but looks well nonetheless in Python.


Until then
Yours Aye
Me

Wednesday, May 13, 2015

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

Thursday, May 7, 2015

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

Sunday, May 3, 2015

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

Wednesday, April 29, 2015

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

Tuesday, April 28, 2015

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




Totient Summatory Function


Applying Dirichlet's Hyberbola method for calculating the Totient summatory function exploits an interesting property of the Totient function. For integer $n$, we have

$\displaystyle\sum\limits_{d|n}\varphi(d)=n$ (or) $\varphi * 1 = n$

Let $f(n)=\varphi(n)$. Then

$F(n)=\displaystyle\sum\limits_{i=1}^n\varphi(i)=\Phi(n)$ and  $\hat{F}(n)=\displaystyle\sum\limits_{i=1}^n\sum_{d|i}\varphi(d)=\sum\limits_{i=1}^n i=\dfrac{n(n+1)}{2}$.

Using these values in Dirichlet's Hyberbola method, we get

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

Solving for the first term of the right summation,

$\Phi(n)=\dfrac{n(n+1)}{2}-\displaystyle\sum\limits_{i=1}^{n/(u+1)} \left\lfloor\frac{n}{i}\right\rfloor \varphi(i) - \sum_{d=2}^u \Phi\left(\left\lfloor\frac{n}{d}\right\rfloor\right) +u^{\text{ }}\Phi\left(\left\lfloor\frac{n}{u+1}\right\rfloor\right)$

So if we memoize our results, we have this beautiful method that enables us to calculate the Totient summatory function. A small hiccup in the above method for calculating $\Phi(n)$ is that you need to precompute the values of $\varphi(k)$ for $k \leq \sqrt{n}$. I have a Mathematica code later along with the timings for reference.

The code takes $70$ seconds to compute $\Phi(10^9)=303963551173008414$, $11$ seconds to compute $\Phi(10^8)=3039635516365908$ and $2$ seconds to compute $\Phi(10^7)=30396356427242$.

But we have another method which will be a bit slower at the expense of not having to precompute the any values. It exploits some of the properties of Euler's totient function and can be generalized to a certain extent. Meanwhile, can you see how we can apply Dirichlet's Hyperbola method to calculate the Mertens function?


Yours Aye
Me

Monday, April 27, 2015

Disclaimer

I'll tell this upfront. I'm no authority on algorithms or computer science to even remotely give you any idea of algorithms. That too, After started reading Donald Knuth's 'Art of Computer Programming', which I could not understand after a certain point, I definitely should not have made this post. You'll find most of these stuffs anyway in Wikipedia.

But should it stop me from sharing the knowledge I have? Blogging about things that I have learned which someone else may be trying to learn? I don't think so. Moreover, rather than seeing things as algorithms I see them as mathematical shortcuts which maybe as much interesting as much they are useful for students of computer science and connoisseurs of Mathematics.

Most of the posts are derived from a list of materials that I collected from the Internet over a period of time and my only intention in creating these posts is to keep all of the algorithms that I've learned in a specific location. I'm simply a student of Mathematics who has great taste for remarkable theorems and algorithms. Nothing more...

A special case of Dirichlet's Hyperbola Method


We saw how to quickly calculate the Divisor summatory function in the last post. I also said something about Dirichlet's Hyperbola method. We'll see how we can generalize the method to other functions.

Let $f(n)$ be a function defined for all integers $n$. Define

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

If you take $f(n)=n^k$, then $F(n)=S_k(n)=1^k+2^k+3^k+\cdots+n^k$ and $\hat{F}(n)=D_k(n)=\displaystyle\sum\limits_{i=1}^n\sigma_k(n)$. Proceeding with the general case as in the last post, we can write

$\hat{F}(n)=\displaystyle\sum\limits_{i=1}^n\left\lfloor\frac{n}{i}\right\rfloor f(i)$

Then using the constancy of the Floor function and the ideas explained in the previous post, it is easy to see that

$\hat{F}(n)=\displaystyle\sum\limits_{i=1}^{n/(u+1)} \left\lfloor\frac{n}{i}\right\rfloor f(i)  + \sum_{d=1}^u d \left( \sum_{1+\lfloor n/(d+1)\rfloor}^{\lfloor n/d \rfloor}f(k)\right)$, $u=\lfloor \sqrt{n}\rfloor$

Using the definition of $F(n)$ now, the equation can be succinctly written as,

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

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

The Divisor summatory function then becomes a direct application of the above formula with $f(n)=n^k$. One of my most favorite application of the above formulae is in finding a recursive relation of the Totient summatory function.

For the generalization, see Dirichlet's Hyperbola method.


Yours Aye
Me

Disclaimer

Divisor Summatory Function


Let $\sigma_k(n)$ denote the sum of the $k$th powers of the divisors of $n$. Or more formally,

$\sigma_k(n)=\sum\limits_{d|n}d^k$

So, $\sigma_0(n)$ will be the number of divisors of $n$, $\sigma_1(n)$ will be the sum of the divisors of $n$, $\sigma_2(n)$ will be the sum of the squares of the divisors of $n$ and so on. For example, the divisors of $15$ are $1,3,5,15$.

$\sigma_0(15)=4$
$\sigma_1(15)=1+3+5+15=24$
$\sigma_2(15)=1^2+3^2+5^2+15^2=260$.

Quoting Divisor function, the divisor function is multiplicative, but not completely multiplicative. It's important for one to understand the difference between multiplicative and completely multiplicative because we will be using these terms in the posts to follow. The Wikipedia page on Divisor function will also give you a formulae to calculate $\sigma_k(n)$ using the factorization of $n$. Now, lets define the Divisor summatory function.

$D_k(n)=\sum\limits_{i=1}^n\sigma_k(i)$

If you go with the straightforward computation of $D_k(n)$, you will find yourself in need the prime factorization of all the numbers below $n$. Now, factorization is one thing that you cannot do easily with a computer program. I mean you can do it, but if you want to do it for a long list of numbers, it's gonna take forever.

If you write a program to find $D_k(n)$, your program will have to iterate through all the numbers less than $n$, find the factorization of each number and then apply the formulae. So the question for the blog is, Is there an easier way of finding $D_k(n)$ without having to factorize all the numbers below $n$? Indeed there is.

Let's deal with $D_1(n)$ for simplicity. This being one of the most common divisor summatory functions, we will simply denote it as $D(n)$. Let's say we want find $D(10)$. Look at the divisors of all the numbers below $10$.

$1:1$
$2:1,2$
$3:1,3$
$4:1,2,4$
$5:1,5$
$6:1,2,3,6$
$7:1,7$
$8:1,2,4,8$
$9:1,3,9$
$10:1,2,5,10$

We see that $1$ is present all the numbers, $2$ is present in every second number, $3$ in every third number and so on. So $D(10)$ will have ten $(=\lfloor\frac{10}{1}\rfloor)$ $1$'s, five $(=\lfloor\frac{10}{2}\rfloor)$ $2$'s, three $(=\lfloor\frac{10}{3}\rfloor)$ $3$'s and so on. Note that, this will work for all integers $n>1$, not just $10$. So we can write,

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

Realigning the sum brings the floor function into the equation. The best thing about the Floor function is that it remains constant over a range of numbers before making a step. Let's take $100$ for example.

$\lfloor\frac{100}{k}\rfloor=\begin{cases}
1,&51\leq k \leq 100\\
2,&34\leq k \leq 50\\
3,&26\leq k \leq 33\\
\vdots\end{cases}$

In light of the above observation, we can write

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

The inner sum can easily be found with a simple formulae.

$\displaystyle\sum\limits_{\lfloor n/(d+1)\rfloor}^{\lfloor n/d \rfloor}k=S_1\left(\left\lfloor\frac{n}{d}\right\rfloor\right)-S_1\left(\left\lfloor\frac{n}{d+1}\right\rfloor\right)$, where $S_1(n)=\dfrac{n(n+1)}{2}$

Note that, $S_k(n)=1^k+2^k+\cdots+n^k$ which can be easily calculated using Faulhaber's formula. The same analysis applies for $D_k(n)$ as well. Following a similar procedure, we have

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

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

What we have used here is a special case of a technique called the Dirichlet's Hyperbola method. Though I will not be discussing the method deeply, much of the techniques that I'll be posting about are someway connected to this method. I'll try to post a Mathematica code here for finding $D_k(n)$ soon.

Until then
Yours Aye,
Me

Disclaimer