• [基本操作] Mobius 反演, Dirichlet 卷积和杜教筛


    Dirichlet 卷积是两个定义域在正整数上的函数的如下运算,符号为 $*$

    $(f * g)(n) = sum_{d|n}f(d)g(frac{n}{d})$

    如果不强调 $n$ 可简写为 $f * g$

    常用:

    $mu * 1 = epsilon$

    $phi * 1 = id$

    $epsilon(n) = [n=1]$

    $id(n)=n$

    Mobius 反演是基于 Dirichlet 卷积的一种....化简式子的方法?

    比较有用的结论就是 $mu * 1 = [n=1]$

    由这个可以引出两个式子

    1.如果 $$F(n) = sum_{n|d}f(d)$$

    则 $$f(n) = sum_{n|d} F(d)mu(lfloor frac{d}{n} floor)$$

    2.如果 $$F(n) = sum_{d|n}f(d)$$

    则 $$f(n) = sum_{d|n} mu(lfloor frac{n}{d} floor)F(d)$$

    还有一个很好用的东西叫做数论分块,即 $lfloor frac{n}{i} floor$ 只有 $sqrt{n}$ 种取值

    知道这些就可以做题了

    bzoj1101 Zap

    求$$sum_{i=1}^n sum_{j=1}^m[gcd(i,j)==k]$$

    sol:先除以 $k$ ,转化为 $$sum_{i=1}^x sum_{j=1}^y[gcd(i,j)==1] space space (x = lfloor frac{n}{k} floor,y=lfloor frac{m}{k} floor)$$

    然后发现 $[gcd(i,j)==1]$ 是一个 $[n=1]$ 的形式,我们把它转化成 $mu * 1$

    得到

    $$sum_{i=1}^x sum_{i=1}^y sum_{d|gcd(i,j)} mu(d)$$

    因为 $d|(gcd(x,y)$等价于 $d|x$ & $d|y$

    而且 $1 hicksim n$ 中满足 $d|x$ 的 $x$ 数量为 $lfloor frac{n}{x} floor$

    所以原式为 $$sum_{d=1}^x mu(d) lfloor frac{x}{d} floor lfloor frac{y}{d} floor$$

    预处理 $mu$ 的前缀和,数论分块即可

    #include<bits/stdc++.h>
    #define LL long long
    using namespace std;
    inline int read()
    {
        int x = 0,f = 1;char ch = getchar();
        for(;!isdigit(ch);ch = getchar())if(ch == '-')f = -f;
        for(;isdigit(ch);ch = getchar())x = 10 * x + ch - '0';
        return x * f;
    }
    const int maxn = 50010;
    int ntp[maxn],pri[maxn],mu[maxn],tot;
    int sum[maxn];
    void getmu()
    {
        mu[1]=1;
        for(int i=2;i<=50000;i++)
        {
            if(!ntp[i])pri[++tot]=i,mu[i]=-1;
            for(int j=1;j<=tot&&pri[j]*i<=50000;j++)
            {
                ntp[i*pri[j]]=1;
                if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
                else mu[i*pri[j]]=-mu[i];
            }
        }
        for(int i=1;i<=50000;i++)sum[i]=sum[i-1]+mu[i];
    } 
    int cal(int x,int y)
    {
        int ret = 0;
        if(x > y)swap(x,y);
        for(int L=1,R=0;L<=x;L=R+1) 
        {
            R = min(x/(x/L),y/(y/L));
            ret += (x/L) * (y/L) * (sum[R] - sum[L - 1]);
        }return ret;
    }
    int main()
    {
        int T = read();
        getmu();
        while(T--)
        {
            int a = read(),b = read(),d = read();
            printf("%d
    ",cal(a / d,b / d));
        }
    }
    View Code

    bzoj2154 Crash的数字表格 && bzoj2693 jzptab

    求 $$sum_{i=1}^n sum_{j=1}^m lcm(i,j)$$

    $n,m leq 10^7$ 对于 bzoj2693 ,有 10000 组询问

    sol:

    首先我们知道 $lcm(i,j) = frac {i imes j}{gcd(i,j)} $

    于是转化成 $$sum_{i=1}^n sum_{j=1}^m frac {i imes j}{gcd(i,j)}$$

    枚举 gcd $$sum_{d=1}^n sum_{i=1}^n sum_{j=1}^m [gcd(i,j) == d]frac {i imes j}{d}$$

    把 d 挪上去 $$sum_{d=1}^n sum_{i=1}^{lfloor n/d floor} sum_{j=1}^{lfloor m/d floor} [gcd(i,j) == 1]i imes j imes d$$

    然后发现和式只有 1 个 d,我们把它拿出来 $$sum_{d=1}^n d imes sum_{i=1}^{lfloor n/d floor} sum_{j=1}^{lfloor m/d floor} [gcd(i,j) == 1]i imes j$$

    换元一下$$sum_{d=1}^n d imes sum_{i=1}^{x} sum_{j=1}^{y} [gcd(i,j) == 1]i imes j space space (x = lfloor n/d floor,y = lfloor m/d floor)$$

    发现后面是个板题,可以枚举 $gcd(i,j)$ 的因数然后反演一波

    $$sum_{d=1}^n d imes sum_{i=1}^{x} sum_{j=1}^{y} i imes j imes sum_{p|i,p| j} mu(p)$$

    把 $p|i$ 和 $p|j$ 这两个条件拆开

    $$sum_{d=1}^{n}d imes sum_{p=1}^{x} mu(p) imes sum_{i=1}^{lfloor frac{x}{p} floor} i imes p imes sum_{j=1}^{lfloor frac{y}{p} floor} j imes p$$

    记 $sum(n) = sum_{i=1}^n i$,用这个推一步就是

    $$sum_{d=1}^{n}d imes sum_{p=1}^{x} mu(p) imes p^2 imes sum(lfloor frac{x}{p} floor) imes sum(lfloor frac{y}{p} floor)$$

    然后把元换回来

    $$sum_{d=1}^{n}d imes sum_{p=1}^{lfloor frac{n}{d} floor} mu(p) imes p^2 imes sum(lfloor frac{n}{d imes p} floor) imes sum(lfloor frac{m}{d imes p} floor)$$

    发现 $lfloor frac{n}{d} floor$ 和 $lfloor frac{m}{d} floor$ 这两个东西也是可以数论分块的,预处理一下 $mu(i) imes i^2$ 就是数论分块套数论分块,时间复杂度是 $O(sqrt{n} imes sqrt{n}) = O(n)$ 的

    做到这就可以做 bzoj 2154 了,但 bzoj 2693 还要再搞一点,毒瘤 bzoj 2693

    #include<bits/stdc++.h>
    using namespace std;
    const int MOD = 20101009,maxn = 1e7 + 20;
    #define LL long long
    inline int read()
    {
        int x = 0,f = 1;char ch = getchar();
        for(;!isdigit(ch);ch = getchar())if(ch == '-') f = -f;
        for(;isdigit(ch);ch = getchar())x = 10 * x + ch - '0';
        return x * f;
    }
    int mu[maxn],pri[maxn],tot,ntp[maxn];
    int n,m;
    int G[maxn],ans;
    int smu[maxn],sqr[maxn];
    void getmu()
    {
        ntp[1]=1;mu[1]=1;
        for(int i=2;i<=n;i++)
        {
            if(!ntp[i]){pri[++tot]=i;mu[i]=-1;}
            for(int j=1;j<=tot&&i*pri[j]<=n;++j)
            {
                ntp[i*pri[j]]=1;
                if(i%pri[j])mu[i*pri[j]]=-mu[i];
                else{mu[i*pri[j]]=0;break;}
            }
        }
        for(int i=1;i<=n;i++)smu[i]=(smu[i-1]+mu[i]+MOD)%MOD;
    }
    int cal(int xx,int yy)
    {
        LL ans=0;
        for(int L=1,R=0;L<=xx;L=R+1)
        {
            R=min(xx/(xx/L),yy/(yy/L));
            int cur = 1LL * (1LL * (1+xx/L) * (xx/L) / 2 % MOD) * (1LL * (1+yy/L) * (yy/L) / 2 % MOD) % MOD;
            (ans += 1LL * (sqr[R] - sqr[L-1]) % MOD * cur % MOD) %= MOD;
        }
        return (ans+MOD)%MOD;
    }
    int main()
    {
        n=read(),m=read();
        if(n>m)swap(n,m);
        getmu();
        for(int i=1;i<=n;i++)sqr[i]=(sqr[i-1]+1ll*i*i%MOD*mu[i]%MOD+MOD)%MOD;
        for(int L=1,R=0;L<=n;L=R+1)
        {
            R=min(n/(n/L),m/(m/L));
            int cur = 1LL * (L + R) * (R - L + 1) / 2 % MOD;
            (ans += 1LL * cal(n / L,m / L) * cur % MOD) %= MOD;
        }
        printf("%d
    ",ans);
        return 0;
    }
    bzoj2154

     令 $t=p imes d$ ,这样就可以把 $p imes d$ 捆在一起枚举,式子变成

    $$sum_{t=1}^n sum(lfloor frac{n}{t} floor) imes sum(lfloor frac{m}{t} floor) imes t imes sum_{d|t} d imes mu(d) $$

    前缀和不太兹磁化简,我们考虑这个式子后面的部分,我们看一看

    $$f(x)=x imes sum_{d|x} d imes mu(d)$$

    经过仔(cha)细(kan)推(ti)理(jie),发现 $f(x)$ 是 $id(x)$ 和 $g(x)=x^2 imes mu(x)$ 的 Dirichlet 卷积,而 $g(x)$ 和 $id(x)$ 都是积性函数,则他们的 Dirichlet 卷积也是积性函数

    考虑线性筛出 $f(x)$

    线性筛的话要考虑两个东西

    1.$f(p)$ p 为质数

    2.$f(p imes q)$ p为质数,p|q

    对于 1. 我们可以观察/打表得出 $f(p) = p - p^2$

    对于 2. 我们考虑 $f(n)$ 与 $f(q)$ 的差别,可以发现 $n$ 比 $q$ 多的因数一定都包含 $p^2$ ,因为 $n=p imes q$,所以后面 $sum_{d|n}d imes mu(d)$ 的值跟 $sum_{d|q}d imes mu(d)$ 的值是一样的,把前面的 $q$ 换成 $n$ 即可

    于是 $f(n) = frac{f(q)}{q} imes n = f(q) imes p$

    将 $f$ 带入原式可以得到答案就是

    $sum_{t=1}^n sum(lfloor frac{n}{t} floor) imes sum(lfloor frac{m}{t} floor) imes f(t)$

    预处理 $f$ 的单点值和前缀和,数论分块即可,单组询问 $O(sqrt{n})$

    #include<bits/stdc++.h>
    #define int long long 
    #define LL long long
    using namespace std;
    inline int read()
    {
        int x = 0,f = 1;char ch = getchar();
        for(;!isdigit(ch);ch = getchar())if(ch == '-') f = -f;
        for(;isdigit(ch);ch = getchar())x = 10 * x + ch - '0';
        return x * f;
    }
    const int mod = 100000009,maxn = 10000050,n = 10000000;
    int ntp[maxn],tot,pri[maxn],f[maxn],sum[maxn];
    void shaker()
    {
        ntp[1] = f[1] = 1;
        for(int i=2;i<=n;i++)
        {
            if(!ntp[i])pri[++tot] = i,f[i] = (((i - (i * i)) % mod) + mod) % mod;
            for(int j=1;j<=tot && i * pri[j] <= n;j++)
            {
                ntp[i * pri[j]] = 1;
                if(i % pri[j] == 0)
                {
                    f[i * pri[j]] = f[i] * pri[j] % mod;
                    break;
                }
                else f[i * pri[j]] = f[i] * f[pri[j]] % mod;
            }
        }
        for(int i=1;i<=n;i++)sum[i] = (sum[i - 1] + f[i]) % mod;
    }
    int sig(int x){return (x * (x + 1) / 2) % mod;}
    int cal(int a,int b)
    {
        int ans = 0;if(a > b)swap(a,b);
        for(int L=1,R=0;L<=a;L=R+1)
        {
            R = min(a / (a / L),b / (b / L));
            ans = (ans + (sum[R] - sum[L - 1] + mod) % mod * sig(a / L) % mod * sig(b / L) % mod) % mod;
        }return ans;
    }
    signed main()
    {
        shaker();
        int T = read();
        while(T--)
        {
            int a = read(),b = read();
            printf("%lld
    ",cal(a,b));
        }
    }
    bzoj2673

    杜教筛用到了 Dirichlet 卷积的性质:如果 $f$ 是积性函数,$g$ 是积性函数,则 $f *g$ 是积性函数

    杜教筛主要用来解决积性函数前缀和的问题,具体方法是给积性函数前缀和卷上一个积性函数 $g$

    快速求出 $g$ 的前缀和和 $f *g$ 的前缀和

    之后分块求原积性函数前缀和

    推导有点复杂,不过也是基础 Mobius 反演

    假设我们需要计算 $S(n) = sum_{i=1}^n f(i)$

    先大力推一波式子

    $$sum_{i=1}^n(f*g)(i) = sum_{i=1}^n sum_{d|i}g(d) imes f(frac{i}{d})$$

    把两个 $sum$ 拆开·,把 $d$ 放到指标上

    $$sum_{d=1}^n g(d) imes sum_{i=1}^{lfloor frac{n}{d} floor}f(i)$$

    发现后面可以用 $S$ 表示
    $$sum_{d=1}^n g(d) imes S(lfloor frac{n}{d} floor)$$

    我们要求 $S(n)$ 也就是 $frac{g(1) imes S(n)}{g(1)}$ ,我们可以先求出 $g(1) imes S(n)$

    也就是$sum_{i=1}^n(f*g)(i) - sum_{i=2}^n g(i) imes S(lfloor frac{n}{i} floor)$

    后面那项可以递归下去,预处理 $S$ 的前 $O(n^{frac{2}{3}})$ 项就可以 $O(n^{frac{2}{3}})$ 的时间筛出 $S(n)$

    bzoj3944 Sum

    求 $phi$ 的前缀和

    求 $mu$ 的前缀和

    10 组询问,每组 n 不超过 INT_MAX

    sol:

    已经知道

    $mu * 1 = epsilon$

    $phi * 1 = id$

    先来算 $phi$ 的前缀和,取 $g = 1$,记 $S(i) = sum_{i=1}^n phi(i)$,则$g(1)S(n) = sum_{i=1}^n i - sum_{i=2}^nS(lfloor frac{n}{i} floor)$,后面的数论分块一下

    $mu$ 的前缀和的话,跟 $phi$ 一样,把$ sum_{i=1}^n i $ 换成 $1$ 就可以了(因为 $epsilon$ 的前缀和是 $1$)

    #include<bits/stdc++.h>
    #define LL long long
    using namespace std;
    inline int read()
    {
        int x = 0,f = 1;char ch = getchar();
        for(;!isdigit(ch);ch = getchar())if(ch == '-') f = -f;
        for(;isdigit(ch);ch = getchar())x = 10 * x + ch - '0';
        return x * f;
    }
    const LL _Bound = 2500000;
    int pri[_Bound + 10],tot;
    LL mu[_Bound + 10];
    LL phi[_Bound + 10];
    char ntp[_Bound + 10];
    map<LL,LL> _phi,_mu;
    void sieve()
    {
        ntp[1] = phi[1] = mu[1] = 1;phi[0] = mu[0] = 0;
        for(LL i=2;i<=_Bound;i++)
        {
            if(!ntp[i])pri[++tot] = i,phi[i] = i - 1,mu[i] = -1;
            for(LL j=1;j<=tot && (LL)pri[j] * i <= (LL)_Bound;j++)
            {
                ntp[i * pri[j]] = 1;
                if(i % pri[j] == 0)
                {
                    mu[i * pri[j]] = 0;
                    phi[i * pri[j]] = phi[i] * pri[j];
                    break;
                }
                else
                {
                    mu[i * pri[j]] = -mu[i];
                    phi[i * pri[j]] = phi[i] * (pri[j] - 1); 
                }
            }
        }
        for(int i=2;i<=_Bound;i++)
        {
            phi[i] += phi[i - 1];
            mu[i] += mu[i - 1];
        }
    }
    inline pair<LL,LL> getans(LL n)
    {
        if(n < _Bound)return make_pair(phi[n],mu[n]);
        map<LL,LL>::iterator it = _phi.find(n);
        if(it != _phi.end())return make_pair(_phi[n],_mu[n]);
        LL ans_phi = n * (n + 1) / 2,ans_mu = 1;
        pair<LL,LL> cur;
        for(LL L=2,R=0;L<=n;L=R+1)
        {
            R = n / (n / L);cur = getans(n / L);
            ans_phi -= (R - L + 1) * cur.first;
            ans_mu  -= (R - L + 1) * cur.second;
        }
        _phi[n] = ans_phi,_mu[n] = ans_mu;
        return make_pair(_phi[n],_mu[n]);
    }/*
    inline LL getmu(LL n)
    {
        if(n < _Bound)return mu[n];
        map<LL,LL>::iterator it = _mu.find(n);
        if(it != _mu.end())return _mu[n];
        LL ans = 1;
        for(LL L=2,R=0;L<=n;L=R+1)
        {
            R = n / (n / L);
            ans -= (R - L + 1) * getmu(n / L);
        }return _mu[n] = ans;
    }*/
    int main()
    {
        int T = read();
        sieve();
        while(T--)
        {
            LL n = read();pair<LL,LL> ans = getans(n);
            printf("%lld %lld
    ",ans.first,ans.second);
        }
    }
    View Code

    bzoj4652 NOI2016 循环之美

    求$$sum_{i=1}^n sum_{j=1}^m [gcd(i,j) == 1][gcd(j,k)==1]$$

    $n,m leq 10^9,k leq 2000$

    (当然原题并没有给出这个式子,手推一会就能推出来的嘛

    sol:

    化简一下式子

    $$sum_{i=1}^n sum_{j=1}^m [gcd(i,j) == 1][gcd(j,k)==1]$$

    把 k 放下来

    $$sum_{i=1}^n sum_{j=1,gcd(j,k)==1}^m [gcd(i,j) == 1]$$

    后面反演一波

    $$sum_{i=1}^n sum_{j=1,gcd(j,k)==1}^m sum{x|gcd(i,j)} mu(x)$$

    转换枚举顺序

    $$sum_{x=1,gcd(x,k)==1}^nmu(x)sum_{i=1}^{lfloor frac{n}{x} floor}sum_{j=1}^{lfloor frac{m}{x} floor} gcd(j,k) == 1$$

    发现有一层的枚举是没啥必要的

    $$sum_{x=1,gcd(x,k)==1}^nmu(x)lfloor frac{n}{x} floorsum_{j=1}^{lfloor frac{m}{x} floor} gcd(j,k) == 1$$

    现在把这个问题转换成了两个前缀和即

    $$f(x)=sum_{i=1}^n[gcd(i,k)==1]$$

    $$g(x)=sum_{i=1}^nmu(i)[gcd(i,k)==1]$$

    第一个式子很好求,根据辗转相除法,$f(n) = f(nspace mod space k) + lfloor frac{n}{k} floor f(k)$

    预处理 $f$ 的前 $k$ 项就可以 $O(1)$ 求了

    对于 $g(x)$ ,我们记它对 k 的值为 $G(n,k)$ ,则有

    $$G(n,k)=sum_{i=1}^nmu(i)[gcd(i,k)==1]$$

    反演

    $$G(n,k)=sum_{i=1}^nmu(i) sum_{x|i,x|k} mu(x)$$

    调一下求和顺序

    $$G(n,k)=sum_{x|k} mu(x) sum_{i=1}^{lfloor frac{n}{x} floor} mu(x imes i)$$

    发现 $mu(x imes i)$ 这一项很有趣,如果 $gcd(i,j)==1$ ,它就是 $mu(i) imes mu(x)$ ,否则它是 $0$ 

    于是就可以用类似“反莫比乌斯反演”(莫比乌斯正演?)的方法把它搞成

    $$G(n,k)=sum_{x|k} [mu(x)]^2 sum_{i=1}^{lfloor frac{n}{x} floor} mu(i)[gcd(i,x)==1]$$

    发现 $sum_{i=1}^{lfloor frac{n}{x} floor} mu(i)[gcd(i,x)==1]$ 很有趣,实际上它就是我们定义的 $G(lfloor frac{n}{x} floor,x)$

    注意到,当 $n = 1$ 的时候它是 $mu$ 的前缀和,这一步可以杜教筛

    然后呢,暴力就可以了,复杂度 $O(飞快)$

    // luogu-judger-enable-o2
    #include<bits/stdc++.h>
    #define LL long long
    using namespace std;
    inline int read()
    {
        int x = 0,f = 1;char ch = getchar();
        for(;!isdigit(ch);ch = getchar())if(ch == '-') f = -f;
        for(;isdigit(ch);ch = getchar())x = 10 * x + ch - '0';
        return x * f;
    }
    int n,m,k;
    const int maxn = 2500001;
    int pri[maxn],pre_f[5010],tot;
    char ntp[maxn];
    LL mu[maxn];
    map<int,LL> _mu;
    map< pair<int,int>,LL > G;
    int ps[50],ToT; 
    void sieve()
    {
        mu[1] = 1;
        for(int i=2;i<maxn;i++)
        {
            if(!ntp[i]){pri[++tot] = i,mu[i] = -1;}
            for(int j=1;j<=tot && pri[j] * i < maxn;j++)
            {
                ntp[i * pri[j]] = 1;
                if(i % pri[j] == 0)
                {
                    mu[i * pri[j]] = 0;
                    break;
                }
                else mu[i * pri[j]] = -mu[i];
            }
        }
        for(int i=2;i<maxn;i++)mu[i] += mu[i - 1];
        for(int i=1;i<=k;i++)pre_f[i] = pre_f[i - 1] + (__gcd(i,k) == 1);
    }
    inline int get_f(int x){return pre_f[x % k] + (x / k) * pre_f[k];}
    inline LL get_mu(int x)
    {
        if(x < maxn)return mu[x];
        if(_mu.find(x) != _mu.end())return _mu[x];
        LL ans = 1;
        for(int L=2,R=0;L<=x;L=R+1)
        {
            R = x / (x / L);
            ans -= (R - L + 1) * get_mu(x / L);
        }return _mu[x] = ans;
    }
    inline LL get_G(int x,int y)
    {
        if(!x)return get_mu(y);
        if(y <= 1)return y;
        if(G.find(make_pair(x,y)) != G.end())return G[make_pair(x,y)];
        return G[make_pair(x,y)] = get_G(x - 1,y) + get_G(x,y / ps[x]);
    }
    int main()
    {
        n = read(),m = read(),k = read();
        sieve();for(int i=1;pri[i]<=k;i++)if(k % pri[i] == 0)ps[++ToT] = pri[i];
        LL ans = 0;
        for(int L=1,R=0;L<=min(n,m);L=R+1)
        {
            R = min(n / (n / L),m / (m / L));
            ans += 1LL * (get_G(ToT,R) - get_G(ToT,L - 1)) * (n / L) * get_f(m / L);  
        }cout<<ans;
    }
    View Code
  • 相关阅读:
    win7 64位安装mongodb及管理工具mongoVUE1.6.9.0
    常见共识算法
    Go语言学习笔记(5)——集合Map
    UPUPW Apache5.5系列本地开发环境配置
    TCP/IP协议
    HTTP协议
    Gossip协议
    《CAP定理》
    比特币双花攻击
    Fabric中的节点类型
  • 原文地址:https://www.cnblogs.com/Kong-Ruo/p/10074901.html
Copyright © 2020-2023  润新知