莫比乌斯反演（持续开坑）

GDOI2018 day2 T1 被吊打后下定决心来学反演。。。

（图片来源于symbol的ppt）

例题

Bzoj2301
https://www.lydsy.com/JudgeOnline/problem.php?id=2301

50%

奇怪的方法就不讲了

所以可以直接求1~x、1~y的答案。

---

设$f(n)$表示1~x、1~y中gcd为n的对数(x< y)
$f(n)=\sum_{i=1}^{x} {\sum_{j=1}^{y} {[gcd(i,j)=n]}}$ ([]是艾佛森括号)

f(n)比较难算，所以再设一个g(n)表示gcd是n倍数的对数
$g(n)=\sum_{i=1}^{x} {\sum_{j=1}^{y} {[n\,|\,gcd(i,j)]}}$ （|是整除）
也可以这样写
$g(n)=\sum_{d=1}^{ ⌊ \frac{x}{n} ⌋ } { f(dn) }$

g(n)很容易算，就等于
$g(n)=⌊\frac{x}{n}⌋⌊\frac{y}{n}⌋$

然后用g乱搞就可以得到50分

70%

然而70分怎么搞。。。

莫比乌斯函数

好吧直接开始讲正题

假设有一个函数g，其表达式为
$g(n)=\sum_{d|n} {f(d)}$
并且g很好求，但f不好求
且最终答案跟f有关

所以要怎么求f？

观察一下g的值
$g(1)=f(1)$
$g(2)=f(1)+f(2)$
$g(3)=f(1)+f(3)$
$g(4)=f(1)+f(2)+f(4)$
$g(5)=f(1)+f(5)$
$g(6)=f(1)+f(2)+f(3)+f(6)$
$g(7)=f(1)+f(7)$
$g(8)=f(1)+f(2)+f(4)+f(8)$
$g(9)=f(1)+f(3)+f(9)$
$g(10)=f(1)+f(2)+f(5)+f(10)$
$g(11)=f(1)+f(11)$
$g(12)=f(1)+f(2)+f(3)+f(4)+f(6)+f(12)$

那么f可以用g求出来
$f(1)=g(1)$
$f(2)=g(2)-g(1)$
$f(3)=g(3)-g(1)$
$f(4)=g(4)-g(2)$
$f(5)=g(5)-g(1)$
$f(6)=g(6)-g(3)-g(2)+g(1)$
$f(7)=g(7)-g(1)$
$f(8)=g(8)-g(4)$
$f(9)=g(9)-g(3)$
$f(10)=g(10)-g(5)-g(2)+g(1)$
$f(11)=g(11)-g(1)$
$f(12)=g(12)-g(6)-g(4)+g(2)$

可以看出f与g之间是有某些规律的
那么可以写成这样
$f(n)=\sum_{d|n} {g(d)*?}$

可以发现?处是与$\frac{n}{d}$有关，将其记作$\mu(\frac{n}{d})$

$\mu$的取值如下：
$\mu(n)= \begin{cases} 1 ,\qquad \quad n=1 \\ (-1)^k \;\,\quad n=p_1*p_2*...*p_k \\ 0 \end{cases}$
其中p为互质的质数

如此便有
$g(n)=\sum_{d|n} {f(d)}$
$f(n)=\sum_{d|n} {g(d)*\mu(\frac{n}{d})}$

或

$g(n)=\sum_{d=1}^{ ⌊ \frac{x}{n} ⌋ } { f(dn) }$
$f(n)=\sum_{d=1}^{ ⌊ \frac{x}{n} ⌋ } { g(dn)*\mu(d) }$
（其实下面的式子原理跟上面类似，都是感性理解容斥原理）

还有一些奇怪的性质（都不会证）

求莫比乌斯函数

实际上可以根据它的性质来线性求出。
不会线筛看这里：https://blog.csdn.net/gmh77/article/details/77413368

线筛有一个性质：当它在$\text{i%p[j]=0}$时就会退出
所以可以以此来判断出每个数是否含有多个相同质因子

因为$\mu$是积性函数，所以$\mu(i)*\mu(p[j])=-\mu(i)$
大概就这样搞，可以线性处理出$\mu$的值

code

void init()
{
    int i,j;

    miu[1]=1;
    sum[1]=1;
    Len=0;

    fo(i,2,len)
    {
        if (!f[i])
        {
            p[++Len]=i;
            miu[i]=-1;
        }

        fo(j,1,Len)
        if (i*p[j]<=len)
        {
            f[i*p[j]]=1;

            if (!(i%p[j]))//如果含有多个相同质因子
            {
                miu[i*p[j]]=0;
                break;
            }

            miu[i*p[j]]=-miu[i];//求μ的值
        }
        else
        break;

        sum[i]=sum[i-1]+miu[i];//求μ的前缀和，很多题都要用到
    }
}

例题1

没错就是刚刚那题……

根据上面50分的做法，
$g(n)=\sum_{d=1}^{ ⌊ \frac{x}{n} ⌋ } { f(dn) }$
可得
$f(n)=\sum_{d=1}^{⌊\frac{x}{n}⌋} {g(dn)*\mu(d)}$
$f(n)=\sum_{d=1}^{⌊\frac{x}{n}⌋} {⌊\frac{x}{dn}⌋⌊\frac{y}{dn}⌋*\mu(d)}$

于是暴力求就有70分了……

然后通过打表观察可以发现
$⌊\frac{n}{d}⌋$最多有$\sqrt{n}$级别种取值
①d<$\sqrt{n}$
因为d最多只有$\sqrt{n}$种取值，所以$⌊\frac{n}{d}⌋$有$\sqrt{n}$种取值
②d>$\sqrt{n}$
因为$⌊\frac{n}{d}⌋$<$\sqrt{n}$，所以共有$\sqrt{n}$种取值

于是分别按照$⌊\frac{x}{dn}⌋$和$⌊\frac{y}{dn}⌋$分成$\sqrt{x}+\sqrt{y}$段，每段的取值都相同，所以可以用$\mu$的前缀和直接求

时间复杂度：$O(\sqrt{x}+\sqrt{y})$

code

//Bzoj2301 - Proble b
#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define min(a,b) (a<b?a:b)
#define len 50000
using namespace std;

int miu[len+1]; //μ
bool f[len+1];
int p[len+1];
int sum[len+1];
int A,B,C,D,K,i,j,k,l,Len,T;
double X,Y;

void swap(int &a,int &b) {int c=a;a=b;b=c;}

void init()
{
    int i,j;

    miu[1]=1;
    sum[1]=1;
    Len=0;

    fo(i,2,len)
    {
        if (!f[i])
        {
            p[++Len]=i;
            miu[i]=-1;
        }

        fo(j,1,Len)
        if (i*p[j]<=len)
        {
            f[i*p[j]]=1;

            if (!(i%p[j]))
            {
                miu[i*p[j]]=0;
                break;
            }

            miu[i*p[j]]=-miu[i];
        }
        else
        break;

        sum[i]=sum[i-1]+miu[i];
    }
}

int Ans(int x,int y)
{
    int i,j,k,ans;

    if (x>y) swap(x,y);

    k=x/K;

    for (ans=0,i=0,j=0; i<k;)
    {
        j=min(x/(x/((i+1)*K))/K,y/(y/((i+1)*K))/K);
        ans+=(sum[j]-sum[i])*(x/(j*K))*(y/(j*K));

        i=j;
    }

    return ans;
}

int main()
{
    init();

    for (scanf("%d",&T); T; T--)
    {
        scanf("%d%d%d%d%d",&A,&B,&C,&D,&K);
        printf("%d\n",Ans(B,D)-Ans(A-1,D)-Ans(B,C-1)+Ans(A-1,C-1));
    }

    return 0;
}

例题2

JZOJ5701. 【gdoi2018 day2T1】谈笑风生
http://172.16.0.132/senior/#contest/show/2372/4

60（80?）%

因为$min(Ax,Ay)<=1000$ （比赛时的题面）
所以可以以小的那边来做，另一边直接分段等差数列求和
（考试时数组开小写炸了）

不难也就想了两个小时

100%

两端的数为n和m，f(d)表示
$f(d)=\sum_{i=1}^{n} {\sum_{j=1}^{m} {(i+j)[gcd(i,j)=d]}}$
f直接求比较蛋疼，考虑设g来辅助
$g(d)=\sum_{i=1}^{n} {\sum_{j=1}^{m} {(i+j)[d\,|\,gcd(i,j)]}}$
$g(d)=\sum_{i=1}^{⌊\frac{n}{d}⌋} {f(di)}$
$f(d)=\sum_{i=1}^{⌊\frac{n}{d}⌋} {g(di)*\mu(i)}$

显然 g可以直接求
$g(d)=\frac{d⌊\frac{n}{d}⌋⌊\frac{m}{d}⌋(2+⌊\frac{n}{d}⌋+⌊\frac{m}{d}⌋)}{2}$
（大概就是两个等差数列乱搞）

那么
$f(d)=\sum_{i=1}^{⌊\frac{n}{d}⌋} {\frac{di⌊\frac{n}{di}⌋⌊\frac{m}{di}⌋(2+⌊\frac{n}{di}⌋+⌊\frac{m}{di}⌋)\mu(i)}{2}}$
因为最后求的是互质情况下的和，所以d=1
所以
$Ans=\sum_{i=1}^{n} {\frac{i⌊\frac{n}{i}⌋⌊\frac{m}{i}⌋(2+⌊\frac{n}{i}⌋+⌊\frac{m}{i}⌋)\mu(i)}{2}}$

同理，分块就可以$O(\sqrt{n}+\sqrt{m})$解决一次询问
（注意维护的是$i*\mu(i)$的前缀和）

code

求边权（只有反演部分）

#include <iostream>
#include <cstdlib>
#include <cstdio>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define min(a,b) (a<b?a:b)
#define Len 100000
using namespace std;

int miu[Len+1];
bool f[Len+1];
int p[Len+1];
int sum[Len+1];
int i,j,n,m,len;
long long ans,x,y;

void init()
{
    int i,j;

    miu[1]=1;sum[1]=1;

    fo(i,2,Len)
    {
        if (!f[i])
        {
            p[++len]=i;
            miu[i]=-1;
        }

        fo(j,1,len)
        if (i*p[j]<=Len)
        {
            f[i*p[j]]=1;

            if (!(i%p[j]))
            {
                miu[i*p[j]]=0;
                break;
            }

            miu[i*p[j]]=-miu[i];
        }
        else
        break;

        sum[i]=sum[i-1]+miu[i]*i;
    }
}

int main()
{
//  freopen("magic.in","r",stdin);
//  freopen("magic.out","w",stdout);

    init();

    scanf("%d%d",&n,&m);
    if (n>m) n^=m,m^=n,n^=m;

    for (ans=0,i=0,j=0; i<n;)
    {
        j=min(n/(n/(i+1)),m/(m/(i+1)));

        x=n/j;
        y=m/j;

        ans=ans+(sum[j]-sum[i])*x*y*(2+x+y)/2;

        i=j;
    }

    printf("%lld\n",ans);

    fclose(stdin);
    fclose(stdout);

    return 0;
}

真·code

完整版

#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cstring>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define min(a,b) (a<b?a:b)
#define max(a,b) (a>b?a:b)
#define Len 100000
using namespace std;

long long A[10001];
long long a[40001][2];
long long ls[10001];
long long b[20001];
long long miu[Len+1];
bool f[Len+1];
long long p[Len+1];
long long sum[Len+1];
long long d[1000001];
long long F[10001];
long long i,j,k,n,m,len,len2,h,t;
long long T,s,l,r,mid;

void New(long long x,long long y,long long s)
{
    len2++;
    a[len2][0]=y;
    a[len2][1]=ls[x];
    b[len2]=s;
    ls[x]=len2;
}

void init()
{
    long long i,j;

    miu[1]=1;sum[1]=1;

    fo(i,2,Len)
    {
        if (!f[i])
        {
            p[++len]=i;
            miu[i]=-1;
        }

        fo(j,1,len)
        if (i*p[j]<=Len)
        {
            f[i*p[j]]=1;

            if (!(i%p[j]))
            {
                miu[i*p[j]]=0;
                break;
            }

            miu[i*p[j]]=-miu[i];
        }
        else
        break;

        sum[i]=sum[i-1]+miu[i]*i;
    }
}

long long get(long long n,long long m)
{
    long long ans,x,y;
    long long i,j;

    if (n>m) n^=m,m^=n,n^=m;

    for (ans=0,i=0,j=0; i<n;)
    {
        j=min(n/(n/(i+1)),m/(m/(i+1)));

        x=n/j;
        y=m/j;

        ans=ans+(sum[j]-sum[i])*x*y*(2+x+y)/2;

        i=j;
    }

    return ans;
}

bool pd(long long s)
{
    long long i,S;

    memset(F,127,sizeof(F));

    h=0;
    t=1;
    d[1]=1;
    F[1]=0;

    while (h<t)
    {
        for (i=ls[d[++h]]; i; i=a[i][1])
        {
            S=max(0,b[i]-s);

            if (F[d[h]]+S<F[a[i][0]])
            {
                if (a[i][0]<n)
                d[++t]=a[i][0];
                F[a[i][0]]=F[d[h]]+S;
            }
        }
    }

    return F[n]<=T;
}

int main()
{
    freopen("magic.in","r",stdin);
    freopen("magic.out","w",stdout);

    init();

    scanf("%lld%lld%lld",&n,&m,&T);
    fo(i,1,n)
    scanf("%lld",&A[i]);

    fo(i,1,m)
    {
        scanf("%lld%lld",&j,&k);
        s=get(A[j],A[k]);

        New(j,k,s);
        New(k,j,s);
    }

    l=0;
    r=607931674418546;

    while (l<r)
    {
        mid=(l+r)/2;

        if (!pd(mid))
        l=mid+1;
        else
        r=mid;
    }

    pd(l);
    printf("%lld %lld\n",l,F[n]);

    fclose(stdin);
    fclose(stdout);

    return 0;
}

例题3

JZOJ4161 于神之怒
http://172.16.0.132/senior/#main/show/4161

O(x)

同样的套路（x< y）
设f(d)表示
$f(d)=\sum_{i=1}^{⌊\frac{x}{d}⌋} {\sum_{j=1}^{⌊\frac{y}{d}⌋} {[gcd(i,j)=1]}}$
设g(d)表示
$g(d)=\sum_{i=1}^{⌊\frac{x}{d}⌋} {\sum_{j=1}^{⌊\frac{y}{d}⌋} {1\,|\,[gcd(i,j)]}}$
$g(d)=\sum_{i=1}^{⌊\frac{x}{d}⌋}{f(di)}$
$f(d)=\sum_{i=1}^{⌊\frac{x}{d}⌋}{g(di)\mu(i)}$

无比显然地得到
$g(d)=⌊\frac{x}{d}⌋⌊\frac{y}{d}⌋$
$f(d)=\sum_{i=1}^{⌊\frac{x}{d}⌋}{⌊\frac{x}{di}⌋⌊\frac{y}{di}⌋\mu(i)}$

那么枚举d，便可得到答案
$Ans=\sum_{d=1}^{x} {d^k\sum_{i=1}^{⌊\frac{x}{d}⌋} {⌊\frac{x}{di}⌋⌊\frac{y}{di}⌋\mu(i)}}$
$Ans=\sum_{d=1}^{x} {d^k\sum_{i=1}^{⌊\frac{x}{d}⌋} {⌊\frac{⌊\frac{x}{d}⌋⌊\frac{y}{d}⌋}{i}⌋\mu(i)}}$

把前后都分一下块来搞，最终复杂度：$O(\sqrt{x}*\sqrt{x})=O(x)$

code

//Jzoj4161 - 于神之怒
#include <iostream>
#include <cstdlib>
#include <cstdio>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define min(a,b) (a<b?a:b)
#define Len 5000000
#define mod 1000000007
using namespace std;

long long miu[Len+1];
long long sum[Len+1];
long long Sum[Len+1];
bool f[Len+1];
long long p[Len+1];
long long N,M,K,len,i,j,k,l;
long long ans;

void swap(long long &x,long long &y) {long long z=x;x=y;y=z;}

long long qpower(long long a,long long b)
{
    long long ans=1;

    while (b)
    {
        if (b&1)
        ans=ans*a%mod;

        a=a*a%mod;
        b>>=1;
    }

    return ans;
}

void init()
{
    long long i,j;

    miu[1]=1;sum[1]=1;Sum[1]=1;
    fo(i,2,N)
    {
        if (!f[i])
        {
            Sum[i]=qpower(i,K);
            p[++len]=i;
            miu[i]=-1;
        }

        long long k=N/i;

        fo(j,1,len)
        if (p[j]<=k)
        {
            f[i*p[j]]=1;
            Sum[i*p[j]]=(Sum[p[j]]*Sum[i])%mod;

            if (!(i%p[j]))
            {
                miu[i*p[j]]=0;
                break;
            }

            miu[i*p[j]]=-miu[i];
        }
        else break;

        sum[i]=sum[i-1]+miu[i];
    }
}

long long Ans(long long d)
{
    long long i,j,k,ans;

    for (ans=0,i=0,j=0,k=N/d; i<k;)
    {
        j=min(N/(N/(d*(i+1)))/d,M/(M/(d*(i+1)))/d);
        ans=(ans+(sum[j]-sum[i])*(N/(d*j))*(M/(d*j)))%mod;

        i=j;
    }

    return ans;
}

int main()
{
    scanf("%lld%lld%lld",&N,&M,&K);
    init();

    if (N>M) swap(N,M);
    fo(i,1,N) Sum[i]=(Sum[i]+Sum[i-1])%mod;

    for (ans=0,i=0,j=0; i<N;)
    {
        j=min(N/(N/(i+1)),M/(M/(i+1)));
        ans=(ans+(Sum[j]-Sum[i])*Ans(j))%mod;

        i=j;
    }

    printf("%lld\n",(ans+mod)%mod);
}

O(sqrt(x))

待补