Skip to main content

NT Part 5: Mobius Function Part 2. Mobius Inversion, Dirichlet Convolution

 Theorem. Let $f(n), g(n)$ be arithmetic function such that $g(n)=\sum_{d|n} f(d)$. Then, $f(n) = \sum_{d|n} \mu (d) g(\frac{n}{d})$.


Proof.$$\sum_{d|n} \mu (d) g(\frac{n}{d}) = \sum_{d|n} \mu (d) \sum_{c|\frac{n}{d}} f(c) = \sum_{d|n} \sum_{c|\frac{n}{d}} \mu (d) f(c) = \sum_{c|n} \sum_{d| \frac{n}{c}} f(c) \mu (d) = \sum_{c|n} f(c) \sum_{d| \frac{n}{c}} \mu (d) = f(n)$$using $\sum_{d|n} \mu (d) = 0$ for all $n>1$ and $\mu (1)=1$.

The Dirichlet Convolution of two function $f, g$ is the following.
$$(f*g) (n) = \sum_{uv=n} f(u)g(v)$$
Denote $1$ as the constant function, $f(n)=1$.
$\epsilon (n)$ is a function which satisfies $\epsilon (1)=1$ and $\epsilon (n) = 0$ for $n \not= 1$.
$Id(n)$ is the identity function, $Id(n)=n$

The mobius inversion formula and the basic property of mobius functions give $1 * \mu = \epsilon$, and $g = f * 1 \iff f=g * \mu$

Okay. So what? How can we use this formula to improve our lives solve problems?

We can, change our "sum" or our answer using mobius inversion and mobius functions to calculate them fast.

I will give 3 examples which changes the desired sum by using number theory (like mobius function) to calculate them fast.

Example Problems

Problem 1. http://www.spoj.com/problems/LCMSUM/

$$\sum_{i=1}^n \text{lcm}(i,n) = n\sum_{i=1}^n \frac{i}{\text{gcd}(i,n)} = n+n \sum_{d|n, d\not= n} \frac{n \phi (\frac{n}{d})}{2d} = \frac{n}{2} + \frac{n}{2} \sum_{d|n} d \phi (d)$$
Let $res[x]$ be the answer when $n=x$.
Precalculate $\phi (d)$, and add $d \phi (d)$ to every $res[i \cdot d]$. Then multiply $\frac{x}{2}$ to $res[x]+1$.
This gives our solution in $O(n \log n + T)$. yay

#include <stdio.h>
#include <iostream>
using namespace std;
typedef long long int ll;
ll res[1000010];
ll phi[1000010];
void preprocess(void) 
{
       for(int i=1 ; i<=1000000 ; i++)
       {
           phi[i]=i;    
    }
    for(int i=2 ; i<=1000000 ; i++)
    {
        if(phi[i]==i) 
        {
            for(int j=1; i*j<=n; j++) 
            {
                phi[i*j]-=phi[i*j]/i;
            }
        }
    }
    for(int i=1 ; i<=1000000 ; i++)
    {
        for(int j=1; i*j<=1000000 ; j++) 
        {
            res[i*j]+=i*phi[i];
        }
    }
}
 
int main(void) 
{
    preprocess();
    int tc;
    cin>>tc;
    while (tc--)
    {
        ll n;
        scanf("%lld",&n);
        ll ans=res[n]+1;
        ans=(ans*n)/2;
        printf("%lld\n",ans);
    }
    return 0;
}


2. https://www.codechef.com/LTIME13/problems/COPRIME3

Here's the key sum transformation.
We want $\sum_{i=1}^n \sum_{j=i+1}^n \sum_{k=j+1}^n \epsilon ( \text{gcd} (a_i,a_j,a_k))$.
Using $1 * \mu = \epsilon$, we can transform this sum. I'll explain this step-by-step.

Start with $\epsilon ( \text{gcd} (a_i,a_j,a_k) )$ as $\sum_{d| \text{gcd} (a_i,a_j,a_k)} \mu (d)$.
Now let's decide how many $\mu (d)$ appears in this sum. Clearly, it is the number of 3-tuples $(i, j, k)$ such that $d|a_i, d|a_j, d|a_k$. Therefore, if $k$ is the number of $a_i$s which is a multiple of $d$, $\mu (d)$ appears $\binom{k}{3}$ times.

Therefore, our result is $\sum_{k=1}^{MAX} \mu (k) \cdot \binom{\text{divcnt}[k]}{3}$, where $\text{divcnt}[k]$ is the number of $a_i$s which are a multiple of $k$.

We can precalculate binomial numbers, $\mu (k)$s, and $\text{divcnt}[k]$ easily in $O(MAX \log MAX)$, where $MAX$ is the maximum of $a_i$s.

#include <stdio.h>
#include <algorithm>
#include <vector>
#include <string.h>
#include <string>
#include <stack>
#include <queue>
#include <iostream>
#include <assert.h>
#include <math.h>
using namespace std;
typedef long long int ll;
ll primechk[111111];
ll mu[111111];
ll divcnt[111111];
ll cnt[111111];
ll com[333][333];
ll ans, mod=1e7+3;
int n;
 
void input(void)
{
    int i, j, x;
    cin>>n;
    for(i=1 ; i<=n ; i++)
    {
        scanf("%d",&x);
        cnt[x]++;
    }
}
 
void divpproc(void)
{
    int i, j, x;
    for(i=1 ; i<=110000 ; i++)
    {
        for(j=1 ; i*j<=110000 ; j++)
        {
            divcnt[i]+=cnt[i*j];
        }
    }
}
 
void compproc(void)
{
    com[0][0]=1;
    int i, j;
    for(i=1 ; i<=320 ; i++)
    {
        com[i][0]=1;
        com[i][i]=1;
    }
    for(i=2 ; i<=320 ; i++)
    {
        for(j=1 ; j<=i-1 ; j++)
        {
            com[i][j]=(com[i-1][j-1]+com[i-1][j]);
        }
    }
}
 
void preprocess(void)
{
    int i, j;
    for(i=1 ; i<=111100 ; i++)
    {
        mu[i]=1;
        primechk[i]=1;
    }
    primechk[1]=0;
    for(i=2 ; i<=111100 ; i++)
    {
        if(primechk[i]==1)
        {
            mu[i]=-mu[i];
            for(j=2 ; i*j<=111100 ; j++)
            {
                primechk[i*j]=0;
                if(j%i==0)
                {
                    mu[i*j]=0;  
                }   
                else
                {
                    mu[i*j]=-mu[i*j];
                }
            }   
        }   
    }
}
 
void calc(void)
{
    int i, j;
    for(i=1 ; i<=105000 ; i++)
    {
            ans=ans+mu[i]*com[divcnt[i]][3];
    }
}
 
int main(void)
{
    preprocess();
    input();
    divpproc();
    compproc();
    calc();
    cout<<ans;
}


3. https://www.codechef.com/COOK29/problems/EXGCD

Here's a sketch. Calculating denominator + calculating inverse of it is just modular inverse calculation.

The main problem is calculating the sum of all gcds. Using $\sum_{d|n} \phi (d) = n$, we can change the sum.

The sum is pretty much $\sum_{g} g \cdot \text{number of tuples that have gcd } g$.

Now change $g$ to $\sum_{d|g} \phi (d)$. Therefore, we can change the sum to $\sum_{d} \phi (d) \cdot \text{number of tuples that have gcd a multiple } d$.

But this is way easier to calculate! For the $gcd$ to be a multiple of $d$, we can just use $\prod_{i=0}^{k-1} \text{number of multiple of } d \text{ in the interval } [A_i,B_i]$. Done!

In the more harder problems, the result of mobius inversion gets more complex, and in some problems we also have to keep track of some prefix sums and use the fact that $\lfloor \frac{n}{i} \rfloor$ takes $O(\sqrt{n})$ values. But that is for later when I can actually solve those problems .

Comments

Post a Comment

Popular posts from this blog

NT Part 2: Generating Primes, Prime Test, Prime Factorization

  Generating primes fast is very important in some problems. Let's cut to the chase and introduce Eratosthenes's Sieve. The main idea is the following. Suppose we want to find all primes between 2 and 50. Iterate from 2 to 50. We start with 2. Since it is not checked, it is a prime number. Now check all numbers that are multiple of    except  2. Now we move on, to number 3. It's not checked, so it is a prime number. Now check all numbers that are multiple of  ,   except  3. Now move on to 4. We see that this is checked - this is a multiple of 2! So 4 is not a prime. We continue doing this. Here's the implementation. #include <stdio.h> int primechk [ 21000 ] ;   void preprocess ( void ) { int i, j ; for ( i = 2 ; i <= 20000 ; i ++ ) { primechk [ i ] = 1 ; } for ( i = 2 ; i <= 20000 ; i ++ ) { if ( primechk [ i ] == 1 ) { for ( j = 2 ; i * j <= 20000...
Post a Comment
Read more

NT Part 6: Sieve of Eratosthenes(Details)

  About 2300 years ago from today, famous Greek mathematician Euclid proved that there are an infinite number of prime numbers. Since then people have been searching for these prime numbers. In 1849, one of the greatest mathematician of all time, Carl Fredrick Gauss, had identified almost all of the prime numbers within the first 3 hundred thousand whole numbers. In the age of computers, we can find large prime numbers in the blink of an eye. But to do that, we need to know a bit of programming and a 2000 year old algorithm. By the end of this tutorial, you will be able to figure out a solution on your own to Gauss’s problem. What is a Prime Number? A prime number is an integer number that is divisible by 1 and the number itself only. For example, 7 is divisible by 1 and 7 only. But 6 is not a prime number because 6 is be divisible by 2 and 3 as well. It is worth mentioning that 1 itself is not a prime number. Now if you are asked to determine if a number is a prime number, you can...
Post a Comment
Read more

NT Part 1:GCD, LCM, Euclidean Algorithm

  The definitions of GCD and LCM are well-known, (and taught in middle school I think) I will skip the definitions. Also, since  ,  calculating GCD is equivalent to calculating LCM. Now, how do we calculate the GCD of two numbers? A naive solution would be iterating over all positive integers no more than  . This will get the GCD in  ,  very very slow. We can calculate the GCD of   in   using Euclidean Algorithm. This algorithm uses the easy-to-prove fact  ,  where   is the remainder when   is divided by  ,  or just a%b. We can now use the following code. #include <iostream> using namespace std ; int gcd ( int u, int v ) { return u % v == 0 ? v : gcd ( v,u % v ) ; } int main ( void ) { int x, y ; cin >> x >> y ; cout << gcd ( x,y ) ; } How do we prove that this algorithm is  ?  Well, let's suppose that we started with  . Then...
Post a Comment
Read more