莫比乌斯反演与狄利克雷卷积和杜教筛感性瞎扯

xtx1092515503 · 2020-02-01 21:52:44 · 个人记录

因为觉得省选将至不能再颓废了，便决定开始更加颓废地研究学也学不会，就算学会也不会用的莫比乌斯反演了。

但我觉得这东西比FFT要好van多了。

1.莫比乌斯函数与莫比乌斯反演

O.约定

主要是因为我太懒了，下取整符号的\LaTeX表达式太长了

I.作用

设有一函数g(n)=\sum\limits_{d|n}f(d)。

显然，如果我们知道每个f(d)的值，那么g(n)的值应该是很容易得到的废话，不就按照定义瞎递推一下吗

那么如果我们知道每个g(n)的值，又应该如何求出f(d)来呢？

II.定义

我们引出莫比乌斯函数\mu(x)。

定义：

设x可以被质因数分解成\prod\limits_{i=1}^n(a_i)^{p_i}，则

换句话说，如果这个x包含某一个质数的平方或更高次数项，则\mu(x)=0；否则，如果它包含偶数个质因数，\mu(x)=1；否则，即它包含奇数个质因数，\mu(x)=-1。

这个形式诡异的函数有什么用吗？

III.性质

1.\mu(x)是积性函数。即，

因此，我们可以采用欧拉筛线性求它。

就是按照它的定义来求，因为积性函数都可以在线性筛时顺便求出。

int mu[N],pri[N];
void getmu(int n){
    mu[1]=1;
    for(int i=2;i<=n;i++){
        if(!pri[i])pri[++pri[0]]=i,mu[i]=-1;
        for(int j=1;i*pri[j]<=n&&j<=pri[0];j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j]))break;
            mu[i*pri[j]]=-mu[i];
        }
    }
}

2.（最重要）

证：

当x=1时，由定义，显然成立。

否则，设x=\prod\limits_{i=1}^n(a_i)^{p_i}。

对于它的每一个因数d，只有当d=\prod\limits_{i=1}^n(a_i)^{p_i},p_i\in\{0,1\}时，才有\mu(d)\neq 0。

如果d包含m个质因数，则\mu(d)=(-1)^m。在所有x的因数中，共有C_n^m个这样的d。

也就是说，\sum\limits_{d|x}\mu(d)=\sum\limits_{m=0}^n(-1)^mC_n^m

依据某奇妙的二项式定理，这个东西有\sum\limits_{m=0}^n(-1)^mC_n^m=(1-1)^n=0。

证毕。

IV.反演

正片开始。

回到我们一开始的说法。莫比乌斯反演对于g(n)=\sum\limits_{d|n}f(d)的函数f和g，假使我们知道每个g(n)的值，我们就可以反推出f(d)的值。

反演定理I：

证：

如果f(x)=\sum\limits_{d|n}g(d)\mu(\dfrac{n}{d})，

则一定有f(x)=\sum\limits_{d|n}g(\dfrac{n}{d})\mu(d)

在\sum\limits_{d|n}g(\dfrac{n}{d})\mu(d)中代入g(x)的定义g(x)=\sum\limits_{d|n}f(d)，

得\sum\limits_{d|n}g(\dfrac{n}{d})\mu(d)=\sum\limits_{d|n}(\sum\limits_{i|\frac{n}{d}}f(i))\mu(d)=\sum\limits_{d|n}(\sum\limits_{i*d|n}f(i))\mu(d)

我们将\mu(d)乘进括号，得

然后发现这个实际上的意思就是枚举所有的i和d，求f(i)\mu(d)的和。

因此有\sum\limits_{d|n}g(\dfrac{n}{d})\mu(d)=\sum\limits_{i|n}(\sum\limits_{i*d|n}f(i)\mu(d))=\sum\limits_{i|n}f(i)(\sum\limits_{i*d|n}\mu(d))

看一下后面的\sum\limits_{i*d|n}\mu(d)，有\sum\limits_{i*d|n}\mu(d)=\sum\limits_{d|\frac{n}{i}}\mu(d)。

是不是很熟悉？当\frac{n}{i}=1，即i=n时，该式值为1；否则，值为0。

则\sum\limits_{i|n}f(i)(\sum\limits_{i*d|n}\mu(d))=f(n)。因为i\neq n时后面的东西都是0。

证毕。

反演定理II：

换成枚举\dfrac{d}{p}，得

代入g(x)=\left\lfloor\dfrac{n}{x}\right\rfloor\times\left\lfloor\dfrac{m}{x}\right\rfloor，得ans=\sum\limits_{p\ is\ prime}\sum\limits_{d=1}^{\min(\left\lfloor\frac{n}{p}\right\rfloor,\left\lfloor\frac{m}{p}\right\rfloor)}\left\lfloor\dfrac{n}{dp}\right\rfloor\times\left\lfloor\dfrac{m}{dp}\right\rfloor\mu(d)

令T=dp，则

则

然后就OK了。只需要将\mu(x)做一个前缀和，套上一个整除分块即可。

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll sum[10001001],ans;
int _1,n,m,pri[10001001],cnt,mu[10001001];
void prime(int ip){
    mu[1]=1,cnt=0;
    for(int i=2;i<=ip;i++){
        if(!pri[i])pri[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt&&i*pri[j]<=ip;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j]))break;
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<=cnt;i++)for(int j=1;j*pri[i]<=ip;j++)sum[j*pri[i]]+=mu[j];
    for(int i=1;i<=ip;i++)sum[i]+=sum[i-1];
}
signed main(){
    scanf("%d",&_1),prime(10000000);
    while(_1--){
        scanf("%d%d",&n,&m),ans=0;
        if(n>m)swap(n,m);
        for(int l=1,r;l<=n;l=r+1)r=min(n/(n/l),m/(m/l)),ans+=1ll*(n/l)*(m/l)*(sum[r]-sum[l-1]);
        printf("%lld\n",ans);
    }
    return 0;
} 

ii.[POI2007]ZAP-Queries

如果前一道题没有听懂的话，是我的锅。毕竟这道题应该放在第一道题，上一道题明显是这道题的升级版。

首先，观察一下题目，发现这道题让我们求的就是上一道题中的f(d)。

我们再来推一下f(d)：

设f(x)为\gcd(i,j)=x的个数，g(x)为\gcd(i,j)是x的倍数的个数。

则据定义，有g(n)=\sum\limits_{n|d}f(d)

据反演定理II，有f(n)=\sum\limits_{n|d}g(d)\mu(\dfrac{d}{n})。

换成枚举\dfrac{d}{n}，并设它为新的d，则有

又有g(x)=\left\lfloor\dfrac{a}{x}\right\rfloor\times\left\lfloor\dfrac{b}{x}\right\rfloor，

则只有当d\leq\min(\left\lfloor\frac{a}{x}\right\rfloor,\left\lfloor\frac{b}{x}\right\rfloor)时，才有g(dn)>0。

因此我们完全可以就求出这个式子的值即可，即

有了这个式子，我们就已经可以做到O(n)的单次询问。

而如果我们把g(dn)的定义代回去的话，就有

换句话说，两个数(i,j)积的因数个数，等于i的所有因数与j的所有因数中互质的对数。

简单证明（感谢suncongbo巨佬的讲解）：

显然各质因数之间相互独立（最后是个数乘在一起，只要每一项都相等，乘积也就相等），我们不如只考虑一个质因数p。设它在i中的次数是a，在j中的次数是b。那么，在乘积i\!j中，次数就是a+b。

考虑仅p一项为i\!j贡献了多少因数。显然，是a+b+1个，次数从0一直到i+j。

那么为等式右侧贡献了多少呢？当i中不取任何p时，j中可以取1到b个p，共b种；当j中不取任何p时，i中可以取1到a个p，共a种；除此之外，还有一种i和j都不选的，共1种。综上，为等式右侧它也贡献了a+b+1。

因为对于每一个单独的p，等式左右的贡献相等；所以综合起来考虑，等式左右贡献的积也相等。

证毕。

然后我们看一看套上这个东西后我们要求什么：

发现这个d可以是任何x的倍数。因此有g(x)=\sum\limits_{i=1}^n\sum\limits_{j=1}^m\left\lfloor\dfrac{n}{i}\right\rfloor\left\lfloor\dfrac{m}{j}\right\rfloor[x|\gcd(i,j)]

我们将x除出来，有

我们令一个sum(x)=\sum\limits_{i=1}^{x}\left\lfloor\dfrac{x}{i}\right\rfloor，这个东西可以在O(n\sqrt{n})时间内通过整除分块写出来。

则g(x)=sum(\dfrac{n}{x})*sum(\dfrac{m}{x})。

又有ans=\sum\limits_{x=1}^{\infty}g(d)\mu(d)，这东西也可以整除分块，因为里面有一个可以整除分块的g(d)和一个可以求前缀和的\mu(d)。

然后就可以了。

代码：

#include<bits/stdc++.h>
using namespace std;
#define int long long
int T,n,m,mu[50100],pri[50100],sum[50100],res;
void getmu(int N){
    mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!pri[i])pri[++pri[0]]=i,mu[i]=-1;
        for(int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j]))break;
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<=N;i++)mu[i]+=mu[i-1];
}
void getsum(int N){
    for(int i=1;i<=N;i++)for(int l=1,r;l<=i;l=r+1)r=i/(i/l),sum[i]+=(r-l+1)*(i/l);
}
signed main(){
    scanf("%lld",&T),getmu(50000),getsum(50000);
    while(T--){
        scanf("%lld%lld",&n,&m);
        if(n>m)swap(n,m);
        res=0;
        for(int l=1,r;l<=n;l=r+1)r=min(n/(n/l),m/(m/l)),res+=(mu[r]-mu[l-1])*sum[n/l]*sum[m/l];
        printf("%lld\n",res);
    }
    return 0;
}

v.[NOI2010]能量采集

真正自己做出来的第一道莫反题祭~~~~

题意：

求\sum_{i=1}^n\sum_{j=1}^m(2\gcd(i,j)-1)。

开始推式子：

到这里为止，就已经可以做了。

复杂度为（设n与m同级）O(\sum_{i=1}^n\sqrt{\dfrac{n}{i}})\approx O(n)。

但实际上还有O(\sqrt{n})的做法（虽然预处理是O(n)的）：

设T=dx，则

这东西不还是O(n)的吗？

我们看一下后面的东西：\sum_{d|T}d\mu(\dfrac{T}{d})，

我们规定一个函数id(x)=x，则这个东西=\sum_{d|T}id(d)\mu(\dfrac{T}{d})。

这个运算，有一个具体的名字，叫做狄利克雷卷积。它是这样的运算：

它的更多性质我们将在接下来一一提到。不过，在这道题中它最大的作用是帮助我们将上面的东西转成狄利克雷卷积的形式，即

好好好，我们搞出了一个新的运算。但是这有什么用吗？

是积性函数就意味着我们可以用欧拉筛去筛它！

而我们已经知道了\mu函数是积性函数。用脑子稍微想想也知道，id函数肯定是积性函数。

这意味着，*$id\mu$也是一个积性函数**，可以用欧拉筛。

而这个东西前一半的东西可以整除分块。

那么我们就丧心病狂地把这个东西优化到了O(\sqrt{n})。

尽管预处理还是O(n)的，但是已经很优秀了。感觉可以出成一道很毒瘤的题

先给出普通的O(n)的代码，而O(\sqrt{n})的我们将在下一题中用到：

#include<bits/stdc++.h>
using namespace std;
#define int long long
int n,m,mu[100100],pri[100100],res;
void getmu(int N){
    mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!pri[i])pri[++pri[0]]=i,mu[i]=-1;
        for(int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j]))break;
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<=N;i++)mu[i]+=mu[i-1]; 
}
signed main(){
    scanf("%lld%lld",&n,&m),getmu(min(n,m));
    for(int i=1;i<=min(n,m);i++){
        int x=n/i,y=m/i,lim=min(x,y),sum=0;
        for(int l=1,r;l<=lim;l=r+1)r=min(x/(x/l),y/(y/l)),sum+=(x/l)*(y/l)*(mu[r]-mu[l-1]);
        res+=2*i*sum;
    }
    res-=n*m;
    printf("%lld\n",res);
    return 0;
}

vi.于神之怒加强版

在这之前，我们引出一个数论函数id^k(x)=x^k。这个函数就是整数域上的k次函数。很明显，它是积性函数，准确地说，是完全积性函数。

它的两个特例，一是k=1，就是我们之前提到的id函数。二是k=0，即id^0函数。因为\forall x\in\mathbb{N}^+,id^0(x)=1，因此这个函数有一个形象的名字：“1函数”，即1(x)。

在上一题中，我们最大的任务是求出\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)，而这个东西，我们有\sqrt{n}复杂度的\sum\limits_{T=1}^{\min(n,m)}\left\lfloor\dfrac{n}{T}\right\rfloor\left\lfloor\dfrac{m}{T}\right\rfloor((id*\mu)(T))。

在这道题中，\gcd(i,j)变成了\gcd^k(i,j)。相应地，如果你把它再证明解一遍的话，就会发现这道题的答案为\sum\limits_{T=1}^{\min(n,m)}\left\lfloor\dfrac{n}{T}\right\rfloor\left\lfloor\dfrac{m}{T}\right\rfloor((id^k*\mu)(T))。

显然，后面的东西是积性函数。我们设g=id^k*\mu。

因为g是积性函数，所以如果有x=\prod\limits_{i=1}^n(P_i)^{a_i}，那么g(x)=\prod\limits_{i=1}^ng((P_i)^{a_i})。

代入g的定义：g(T)=\sum_{d|T}id^k(d)\mu({\dfrac{T}{d}})，得g(x)=\prod\limits_{i=1}^n\sum_{d|(P_i)^{a_i}}id^k(d)\mu({\dfrac{(P_i)^{a_i}}{d}})。

当d<P_i^{a_i-1}时，\mu({\dfrac{(P_i)^{a_i}}{d}})=0。因此我们只用考虑d=P_i^{a_i-1}和P_i^{a_i}。

则g(x)=\prod\limits_{i=1}^nid^k(P_i^{a_i-1})\mu(P_i)+id^k(P_i^{a_i})\mu(1)。

代入id^k和\mu的定义，得g(x)=\prod\limits_{i=1}^n(P_i^{(a_i-1)*k})*(-1)+(P_i^{a_i*k})*1

则有g(x)=\prod\limits_{i=1}^n(P_i^{(a_i-1)*k})*(P_i^{a_i}-1)。依照这个定义，我们就可以写出线性筛的程序了。

代码：

#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
int T,n,m,k,g[5001000],f[5001000],pri[5001000],res;
int ksm(int x,int y){
    int res=1;
    for(;y;x=(1ll*x*x)%mod,y>>=1)if(y&1)res=(1ll*res*x)%mod;
    return res;
}
void getg(int N){
    g[1]=1;
    for(int i=2;i<=N;i++){
        if(!pri[i])pri[++pri[0]]=i,f[pri[0]]=ksm(i,k),g[i]=(f[pri[0]]-1+mod)%mod;
        for(int j=1;j<=pri[0]&&1ll*i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j])){
                g[i*pri[j]]=1ll*g[i]*f[j]%mod;
                break;
            }
            g[i*pri[j]]=1ll*g[pri[j]]*g[i]%mod;
        }
    }
    for(int i=1;i<=N;i++)g[i]=(g[i-1]+g[i])%mod;
}
int main(){
    scanf("%d%d",&T,&k),getg(5000000);
    while(T--){
        scanf("%d%d",&n,&m),res=0;
        for(int l=1,r;l<=min(n,m);l=r+1)r=min(n/(n/l),m/(m/l)),(res+=1ll*(g[r]-g[l-1]+mod)%mod*(n/l)%mod*(m/l)%mod)%=mod;
        printf("%d\n",res);
    }
    return 0;
}

vii.[SDOI2014]数表

仍然是线性筛筛各种东西。我们引出一个东西\sigma(x)=\sum\limits_{d|n}d，也就是x的约数和。这个东西明显是积性函数。设x=\prod\limits_{i=1}^n(P_i)^{a_i}，则\sigma(x)=\prod\limits_{i=1}^n(\sum\limits_{j=0}^{a_i}(P_i)^j)。

如果我们忽略a的限制，则我们要求的就是\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sigma(\gcd(i,j))。

开始推式子。

假使我们忽略a的限制，到现在我们已经做完了。但是现在有a的限制，怎么办呢？

设g=\mu*\sigma。

将询问按照a排序，同时用树状数组动态地对于不同的a更新g函数。比如说，假设上一个a是5，这一个a是6，那对于所有\sigma(x)=6的x，枚举在10^5范围内所有x的倍数T，更新g(T)。

而我们最后统计时，就用树状数组维护g函数的区间和。

哦，对了，关于取模，直接自然溢出即可。

关于线性筛筛各种东西，我推荐一篇blog。里面就有线性筛筛\sigma(x)的讲解。

代码：

#include<bits/stdc++.h>
using namespace std;
const int N=100000;
int T,pri[100100],mu[100100],sigma[100100],low[100100],sum[100100],t[100100],res[100100];;
pair<int,int>p[100100];
void sieve(){
    mu[1]=sigma[1]=1;
    for(int i=2;i<=N;i++){
        if(!pri[i])low[i]=pri[++pri[0]]=i,mu[i]=-1,sum[i]=sigma[i]=i+1;
        for(int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j])){
                mu[i*pri[j]]=0;
                low[i*pri[j]]=low[i]*pri[j];
                sum[i*pri[j]]=sum[i]+low[i*pri[j]];
                sigma[i*pri[j]]=sigma[i]/sum[i]*sum[i*pri[j]];
                break;
            }
            mu[i*pri[j]]=-mu[i];
            low[i*pri[j]]=pri[j];
            sum[i*pri[j]]=pri[j]+1;
            sigma[i*pri[j]]=sigma[i]*sigma[pri[j]];
        }
    }
    for(int i=1;i<=N;i++)p[i]=make_pair(sigma[i],i);
    sort(p+1,p+N+1);
}
void add(int x,int y){
    while(x<=N)t[x]+=y,x+=x&-x;
}
int ask(int x){
    int res=0;
    while(x)res+=t[x],x-=x&-x;
    return res;
}
struct query{
    int n,m,a,id;
    friend bool operator <(const query &x,const query &y){
        return x.a<y.a;
    }
}q[100100];
int solve(int n,int m){
    if(n>m)swap(n,m);
    int ans=0;
    for(int l=1,r;l<=n;l=r+1)r=min(n/(n/l),m/(m/l)),ans+=(n/l)*(m/l)*(ask(r)-ask(l-1));
    return ans;
}
int main(){
    scanf("%d",&T),sieve();
    for(int i=1;i<=T;i++)scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].a),q[i].id=i;
    sort(q+1,q+T+1);
    for(int i=1,j=1;i<=T;i++){
        while(j<=N&&p[j].first<=q[i].a){
            for(int k=p[j].second;k<=N;k+=p[j].second)add(k,sigma[p[j].second]*mu[k/p[j].second]);
            j++;
        }
        res[q[i].id]=solve(q[i].n,q[i].m);
    }
    for(int i=1;i<=T;i++)printf("%d\n",res[i]&(~(1<<31)));
    return 0;
} 

viii.[SDOI2017]数字表格

题意：求出

观察它的指数\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]，正是我们这几题来最常见的反演形式。它等于\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}\sum_{x|i,x|j}\mu(x)=\sum_{x=1}^{\min(n/d,m/d)}\mu(x)\left\lfloor\dfrac{n}{dx}\right\rfloor\left\lfloor\dfrac{m}{dx}\right\rfloor。

也就是说，

设\prod_{d|T}(f_d)^{\mu(T/d)}=F(T)，

则有\prod\limits_{i=1}^n\prod\limits_{j=1}^mf_{\gcd(i,j)}=\prod\limits_{T=1}^{\min(n,m)}F(T)^{\left\lfloor n/T\right\rfloor\left\lfloor m/T\right\rfloor}。

我们还有其它性质，例如：

提到狄利克雷卷积，我们就不得不提它所作用着的载体：数论函数。

数论函数是指定义域为\mathbb{N}^+，值域为\mathbb{Z}的全体函数。

常见数论函数包括（设x=\prod_{i=1}^n(P_i)^{a_i}）：

（符号不唯一）

这些函数至少都是积性函数，部分是完全积性函数。

我们之前接触过\mu的一个性质：

或者说，

用上我们新学的狄利克雷卷积的知识，这就是：

换句话说，\mu和1互为逆元。

这样，我们就可以通过另一个途径证出莫比乌斯反演：

证毕。

而欧拉函数\varphi也有如此性质：

换句话说，

我们发现，\mu和\varphi在卷上1后都会得到一些奇妙的性质。那么这两者之间有何关联呢？

瞧瞧我们得到了什么好van的东西！

虽然这东西也没啥用，只是很美妙罢了。

3.杜教筛

之前在做莫反的题时，有很多题都需要用到杜教筛，因而我非常不爽。因此便来研究杜教筛了。

杜教筛可以干什么？

在非线性时间内（准确说，O(n^{\frac{2}{3}})）求出某些积性函数的前缀和。例如，\sum_{i=1}^n\mu(i)。

怎么办呢？

假设我们要求S(n)=\sum_{i=1}^nf(i)。这时，我们随便找出一个函数g。我们计算\sum_{i=1}^n(f*g)(i)。

则g(1)S(n)=\sum_{d=1}^ng(d)S(\dfrac{n}{d})-\sum_{d=2}^ng(d)S(\dfrac{n}{d})=\sum_{d=1}^n(f*g)(d)-\sum_{d=2}^ng(d)S(\dfrac{n}{d})

后面的那一大坨可以整除分块，递归地跑下去。这样，如果我们可以找出一个合适的g，可以O(\sqrt{n})以内地算出\sum_{d=1}^n(f*g)(d)，就能以一个较快的速度跑出S(n)。

回到我们的例子。\sum_{i=1}^n\mu(i)。

回忆起1*\mu=\epsilon，我们可以尝试使g=1。

则g(1)S(n)=\sum_{d=1}^n\epsilon(d)-\sum_{d=2}^ng(d)S(\dfrac{n}{d})

又\forall (x),g(x)=1，因此S(n)=1-\sum_{d=2}^nS(\dfrac{n}{d})。

这样我们就找出了一种好方法。

考虑当n较小时，递归计算的复杂度甚至还要劣于线性筛。可以证明，当我们预处理出n^{\frac{2}{3}}以内的S(i)时，复杂度达到最优，为O(n^{\frac{2}{3}})。

我们可以用手写哈希表，或者unordered_map，来实现记忆化搜索。

再考虑\sum_{i=1}^n\varphi(i)。

这次，我们还是选择g=1。因为1*\varphi=id，所以\sum_{i=1}^n(f*g)(i)=\sum_{i=1}^nid(i)=\sum_{i=1}^ni=\dfrac{n(n+1)}{2}。

这样我们就可以通过【模板】杜教筛（Sum）了。

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=5000000;
int T,n,pri[N+5],mu[N+5];
ll phi[N+5];
void init(){
    phi[1]=mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!pri[i])pri[++pri[0]]=i,mu[i]=-1,phi[i]=i-1;
        for(int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j])){
                phi[i*pri[j]]=phi[i]*pri[j];
                break;
            }else phi[i*pri[j]]=phi[i]*phi[pri[j]],mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<=N;i++)mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}
unordered_map<int,int>mm;
unordered_map<int,ll>mp;
int getmu(int x){
    if(x<=N)return mu[x];
    if(mm[x])return mm[x];
    int res=1;
    for(int l=2,r;l<=x;l=r+1)r=x/(x/l),res-=getmu(x/l)*(r-l+1);
    return mm[x]=res;
}
ll getphi(int x){
    if(x<=N)return phi[x];
    if(mp[x])return mp[x];
    ll res=1ll*x*(x+1)/2;
    for(int l=2,r;l<=x;l=r+1)r=x/(x/l),res-=getphi(x/l)*(r-l+1);
    return mp[x]=res;
}
int main(){
    init(),scanf("%d",&T);
    while(T--)scanf("%d",&n),printf("%lld %d\n",getphi(n),getmu(n));
    return 0;
}

I.简单的数学题

在做这题之前，我们先来见一位老朋友：

我们在1.v.[NOI2010]能量采集中就已经接触到了这道题。当时我们运用了\sum\limits_{d|n}\mu(d)=[n=1]的公式。现在，我们要运用另一个公式\sum\limits_{d|n}\varphi(d)=n。

前面的东西可以线性筛前缀和，后面的东西可以整除分块，这样我们也得到一种和之前运用的\mu复杂度一致的运算。

而现在这道题呢？

我们要求\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\gcd(i,j)。

再来！

后面的(\sum_{i=1}^{n/x}i)^2可以直接套等差数列求和公式并整除分块。前面的我们考虑杜教筛。

设f(x)=\varphi(x)x^2，我们要求S(n)=\sum_{i=1}^n\varphi(i)i^2。

我们想找一个合适的g。

考虑(f*g)(n)=\sum\limits_{d|n}\varphi(d)d^2g(\dfrac{n}{d})。

我们想把这个恶心的d^2消掉。这样的话，我们不如设g=id^2。这样的话，(f*g)(n)=\sum\limits_{d|n}\varphi(d)d^2(\dfrac{n}{d})^2=\sum\limits_{d|n}\varphi(d)n^2=n^3！！！

我们好像找对了。

再来看式子：

代入，得S(n)=\sum_{i=1}^ni^3-\sum_{d=2}^ni^2S(\dfrac{n}{d})

我们有公式：

这样我们就可以杜教筛了。

我们回到一开始的式子：\sum_{x=1}^n\varphi(x)x^2(\sum_{i=1}^{n/x}i)^2。后半部分是O(\sqrt{n})的整除分块，前半部分要用杜教筛，这样，看起来复杂度为O(n^{\frac{2}{3}}*\sqrt{n})=O(n^{\frac{7}{6}})。

Emm?这东西一看就不像能过的样子呀？

首先，因为杜教筛是基于记忆化搜索实现的，多次使用复杂度明显是低于O(n^{\frac{2}{3}})的。至于低多少，就要看出题人是否是用脚造数据了。

同时，因为整除分块的过程中，杜教筛要筛的n不会总是n级别的。总复杂度应该低于O(\sum_{i=1}^n(\dfrac{n}{i})^{\frac{2}{3}})

并且，因为杜教筛是基于分块思想的，我们一开始预处理的部分也可以大于n^{\frac{2}{3}}，这样单次杜教筛复杂度就会低于O(n^{\frac{2}{3}})。这样一番优化下来，这个算法复杂度应该是亚线性，甚至直接是O(n^{\frac{2}{3}})的（其实是我不会证）

或者换一种说法，它的复杂度是O(\text{能过})。

代码：

#pragma GCC optimize(3)
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=4000000;
int p,phi[N+10],inv6,pri[N+10],n;
inline int ksm(int x,int y){
    register int res=1;
    for(;y;x=(1ll*x*x)%p,y>>=1)if(y&1)res=(1ll*res*x)%p;
    return res;
}
inline void init(){
    inv6=ksm(6,p-2);
    phi[1]=1;
    for(register int i=2;i<=N;i++){
        if(!pri[i])pri[++pri[0]]=i,phi[i]=i-1;
        for(register int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j])){phi[i*pri[j]]=(1ll*phi[i]*pri[j])%p;break;}
            phi[i*pri[j]]=phi[i]*phi[pri[j]];
        }
    }
    for(register int i=1;i<=N;i++)phi[i]=(phi[i-1]+1ll*phi[i]*i%p*i%p)%p;
}
inline int sqrsum(int x){
    x%=p;
    return 1ll*x*(x+1)%p*(2*x+1)%p*inv6%p;
}
inline int cubsum(int x){
    x%=p;
    return (1ll*((1ll*x*(x+1)/2)%p)*((1ll*x*(x+1)/2)%p))%p;
}
unordered_map<int,int>mp;
inline int djs(int x){
    if(x<=N)return phi[x];
    if(mp[x])return mp[x];
    register int res=cubsum(x);
    for(register int l=2,r;l<=x;l=r+1)r=x/(x/l),res=(res-(1ll*(sqrsum(r)-sqrsum(l-1)+p)%p*djs(x/l)%p)+p)%p;
    return mp[x]=res;
}
inline int solve(){
    register int res=0;
    for(register int l=1,r;l<=n;l=r+1)r=n/(n/l),res=(1ll*(djs(r)-djs(l-1)+p)%p*cubsum(n/l)%p+res)%p;
    return res;
}
signed main(){
    scanf("%lld%lld",&p,&n),init();
    printf("%lld\n",solve());
    return 0;
} 

II.[CQOI2015]选数

我们要求这个东西：

老套路，除一下，得到

把x拖出去，得到

后面的东西整除分块掉，前面的东西杜教筛掉。

然后就可以了。

复杂度O(\text{能过})。

代码：

#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
const int N=5000000;
int n,k,L,R,mu[N+10],pri[N+10],ans;
void getmu(){
    mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!pri[i])pri[++pri[0]]=i,mu[i]=-1;
        for(int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j]))break;
            mu[i*pri[j]]=-mu[i];
        }
    }
//  for(int i=1;i<=10;i++)printf("%d\n",mu[i]);
    for(int i=1;i<=N;i++)mu[i]=(0ll+mu[i]+mu[i-1]+mod)%mod;
//  for(int i=1;i<=10;i++)printf("%d\n",mu[i]);
}
unordered_map<int,int>mp;
int djs(int x){
    if(x<=N)return mu[x];
    if(mp[x])return mp[x];
    int res=1;
    for(int l=2,r;l<=x;l=r+1)r=x/(x/l),res=(0ll+res-(1ll*(r-l+1)*djs(x/l)%mod)+mod)%mod;
    return mp[x]=res;
}
int ksm(int x,int y){
    int res=1;
    for(;y;x=(1ll*x*x)%mod,y>>=1)if(y&1)res=(1ll*res*x)%mod;
    return res;
}
int main(){
    scanf("%d%d%d%d",&n,&k,&L,&R),L--,L/=k,R/=k,getmu();
//  printf("%d %d\n",L,R);
    for(int l=1,r;l<=R;l=r+1){
        if(!(L/l))r=R/(R/l);
        else r=min(L/(L/l),R/(R/l));
        ans=(ans+(1ll*(djs(r)-djs(l-1)+mod)%mod*ksm((R/l)-(L/l),n)%mod))%mod;
    }
    printf("%d\n",ans);
    return 0;
}

III.DZY Loves Math

题意：求

但是，稍微想想就发现，f函数并不是积性函数。例如，f(5)=1,f(6)=1,\gcd(5,6)=1,f(5*6)=f(30)=1。

那怎么办呢？

考虑(f*\mu)(x)=\sum\limits_{d|x}\mu(d)f(\dfrac{x}{d})。

我们设x=\prod_{i=1}^k(P_i)^{a_i}。

显然，只有当d=\prod_{i=1}^k(P_i)^{b_i},b_i\in\{0,1\}时，才有\mu(d)\neq 0。

我们设某个具有\min(a_i)的P_i为p。

再设一个d|x但\gcd(d,p)=1的d。

则，\mu(d)=-\mu(dp),f(\dfrac{x}{d})=f(\dfrac{x}{dp})

则\mu(d)f(\dfrac{x}{d})+\mu(dp)f(\dfrac{x}{dp})=0

因此，大部分情况下，(f*\mu)(x)=0

唯一的特例就是所有的a_i全部相等时。这时，就算是\min(a_i)也会对f有影响。

当然，除了d=1和d=\prod_{i=1}^k(P_i)两个值以外，其他的\mu(d)f(\dfrac{x}{d})还是会互相抵消掉。

当d=1时，\mu(d)f(\dfrac{x}{d})=a_i；当d=\prod_{i=1}^k(P_i)时，\mu(d)f(\dfrac{x}{d})=(a_i-1)*\mu(d)

这两个求和的话，我们最终得到了-\mu(d)。

我们设f*\mu=g。然后就只要枚举每个\mu(x)\neq 0的x，将它的所有次方的g设为-\mu(x)即可。这个复杂度为O(n)，因为每一个g(x)只会被访问一次。

虽然这是莫反题，但本题最大的难点是线性时间内求出f*\mu的值。数学素养是很重要的。（就像我，根本求不出来f*\mu到底是个什么东西）

代码：

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=1e7;
int n,m,pri[N+5],mu[N+5],g[N+5],T,res;
void init(){
    mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!pri[i])pri[++pri[0]]=i,mu[i]=-1;
        for(int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j]))break;
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=2;i<=N;i++)if(mu[i])for(int j=i;j<=N;j*=i)g[j]=-mu[i];
    for(int i=1;i<=N;i++)g[i]+=g[i-1];
}
signed main(){
    scanf("%lld",&T),init();
    while(T--){
        scanf("%lld%lld",&n,&m),res=0;
        for(int l=1,r;l<=min(n,m);l=r+1)r=min(n/(n/l),m/(m/l)),res+=(g[r]-g[l-1])*(n/l)*(m/l);
        printf("%lld\n",res);
    }
    return 0;
}

IV.Lcm

既然上一道题中的DZY都能自定义函数，那我们为什么不能呢？

定义f(x)为x中是否含有平方项。没有则为1，有则为0。显然，它是积性函数。而我们要求的，就是

都是老一套了。做多了就发现都是这个德行。

我们设sum(x)=\dfrac{x(x+1)}{2}。

我们设g(T)=T\sum_{d|T}f(d)(\dfrac{T}{d})\mu(\dfrac{T}{d})。因为f是积性函数，\mu(x)*id(x)（注意这里是数乘不是卷积）也是积性函数，因此g是积性函数。

我们考虑线性筛g。

设p为一质数,则g(p)=p*(f(1)*p*\mu(p)+f(p)*1*\mu(1))=p*(1-p)

考虑如何求出g(xp)，其中p为一质数，x为任意数（保证大于p）。

当\gcd(x,p)=1，即x,p互质时，我们按照积性函数性质，有g(xp)=g(x)*g(p)。

否则，如果xp中含有三个及以上的因数p，则f(xp)=0，可以直接得出g(xp)=0。

否则，我们将xp分割成p^2和\dfrac{x}{p}两个必然互质的部分，运用积性函数性质，则有g(xp)=g(\dfrac{x}{p})*g(p^2)。

考虑如何求出g(p^2)。按照性质递推，会发现它等于-p^3。

当然咯，为了减少计算量，因为g(T)=T\sum_{d|T}f(d)(\dfrac{T}{d})\mu(\dfrac{T}{d})，因此这个因子T我们统一留到最后加上。这时，有g(p)=(1-p),g(p^2)=-p（复辟）。

然后就直接整除分块即可。

代码：

#include<bits/stdc++.h>
using namespace std;
const int N=4e6;
int n,m,T,pri[N+5],g[N+5],res;
void init(){
    g[1]=1;
    for(int i=2;i<=N;i++){
        if(!pri[i])pri[++pri[0]]=i,g[i]=1-i;
        for(int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j])){
                int r=i/pri[j];
                if(r%pri[j])g[i*pri[j]]=-pri[j]*g[r];
                break;
            }
            g[i*pri[j]]=g[i]*g[pri[j]];
        }
    }
    for(int i=1;i<=N;i++)g[i]=g[i-1]+g[i]*i;
}
int main(){
    scanf("%d",&T),init();
    while(T--){
        scanf("%d%d",&n,&m),res=0;
        for(int l=1,r;l<=min(n,m);l=r+1)r=min(n/(n/l),m/(m/l)),res+=((1ll*(n/l)*(n/l+1)/2)*(1ll*(m/l)*(m/l+1)/2)*(g[r]-g[l-1]));
        printf("%d\n",res&(~((1<<30)|(1<<31))));
    }
    return 0;
}

V.Product

要求这个东西：

开始推式子。

上下分别考虑。首先，我们有\prod_{i=1}^n\prod_{j=1}^nij=(n!)^{2n}。

考虑下面的东西。

然后我们最后有

下面的东西，是可以\sqrt{n}\log n或者类似的亚线性时间里，用两次整除分块加快速幂做出来的。但是，如果这样的话，就需要预处理出阶乘数组。但是毒瘤出题人卡空间，8M几乎开不了什么东西。先把一份空间8.4M，刚好被卡掉的正确代码奉上：

#include<bits/stdc++.h>
using namespace std;
const int N=1e6;
const int mod=104857601;
int n,pri[N+2],mu[N+2],ans=1;
void init(){
    mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!pri[i])pri[++pri[0]]=i,mu[i]=-1;
        for(int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j]))break;
            mu[i*pri[j]]=-mu[i];
        }
    }
    pri[0]=1;
    for(int i=1;i<=N;i++)pri[i]=1ll*pri[i-1]*i%mod,mu[i]=(0ll+mu[i]+mu[i-1]+mod-1)%(mod-1);
}
int ksm(int x,int y){
    int res=1;
    for(;y;x=(1ll*x*x)%mod,y>>=1)if(y&1)res=(1ll*res*x)%mod;
    return res;
}
int calc(int x){
    int res=0;
    for(int l=1,r;l<=x;l=r+1)r=x/(x/l),res=(res+1ll*(mu[r]-mu[l-1]+(mod-1))%(mod-1)*(x/l)%(mod-1)*(x/l)%(mod-1))%(mod-1);
    return res;
}
int main(){
    scanf("%d",&n),init();
    for(int l=1,r;l<=n;l=r+1)r=n/(n/l),ans=(1ll*ans*ksm(1ll*pri[r]*ksm(pri[l-1],mod-2)%mod,calc(n/l)))%mod;
    ans=(1ll*ans*ans)%mod;
    ans=ksm(ans,mod-2);
    ans=1ll*ans*ksm(pri[n],2*n)%mod;
    printf("%d\n",ans);
    return 0;
} 

为了避免MLE，我们只能放弃外层的分块，只保留内层的分块。则复杂度为\sum\limits_{i=1}^n\sqrt{\dfrac{n}{i}}=n\log n。

代码：

#include<bits/stdc++.h>
using namespace std;
const int N=1e6;
const int mod=104857601;
int n,pri[80100],mu[N+2],ans=1,ans2=1;
bool vis[N+2];
void init(){
    mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!vis[i])pri[++pri[0]]=i,mu[i]=-1;
        for(int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            vis[i*pri[j]]=true;
            if(!(i%pri[j]))break;
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<=N;i++)mu[i]=(0ll+mu[i]+mu[i-1]+(mod-1))%(mod-1);
}
int ksm(int x,int y){
    int res=1;
    for(;y;x=(1ll*x*x)%mod,y>>=1)if(y&1)res=(1ll*res*x)%mod;
    return res;
}
int calc(int x){
    int res=0;
    for(int l=1,r;l<=x;l=r+1)r=x/(x/l),res=(1ll*(mu[r]-mu[l-1]+(mod-1))%(mod-1)*(x/l)%(mod-1)*(x/l)%(mod-1)+res)%(mod-1);
    return res;
}
int main(){
    scanf("%d",&n),init();
    for(int i=1;i<=n;i++)ans=1ll*ans*i%mod;
    ans=ksm(ans,2*n);
    for(int i=1;i<=n;i++)ans2=(1ll*ans2*ksm(i,calc(n/i)))%mod;
    ans2=(1ll*ans2*ans2)%mod;
    ans2=ksm(ans2,mod-2);
    printf("%d\n",1ll*ans*ans2%mod);
    return 0;
} 

VI.LJJ爱数数

题目给出要求这样的东西

开始胡搞

设(a-c)=x,(b-c)=y。如果\gcd(a,b,c)> 1的话，则一定满足\gcd(a,b,c)|a,\gcd(a,b,c)|b,\gcd(a,b,c)|c，则又有\gcd(a,b,c)|(a-c),\gcd(a,b,c)|(b-c)，最终有\gcd(x,y)>1。反之亦然。

这样，只要保证\gcd(x,y)=1即可。

这样的话，如果将c^2的全部质因子分配给a-c和b-c的话，相同的质因子必然全部分给了一个数，不然必有\gcd(x,y)>1。

因此，必有x,y为完全平方数。

我们就设i=\sqrt{x},j=\sqrt{y}。则有c=ij。

因为有n这个限制，我们必有\max(a,b)\leq n，不妨设a\geq b，则必有x\geq y,i\geq j。

如果我们枚举i的话，j就可以是所有满足c\leq n,j\leq i，且\gcd(i,j)=1的值。

我们看j最多可以在什么范围内取。必有x+c\leq n，即i^2+ij\leq n，即j\leq \dfrac{n}{i}-i。

同时，为了避免重复计算，我们还有限制j<i。

这样的话，j的上限就是\min(\dfrac{n}{i}-i,i-1)。设为b_i。

则我们要求\sum\limits_{i=1}^{\sqrt{n}}\sum\limits_{j=1}^{b_i}[\gcd(i,j)=1]

老套路了吧。不过注意，这回j的上限b_i与i有关，因此不能直接把它变成下取整的形式，必须保留这个东西。

直接没有任何技术含量地暴力即可。复杂度为\sum_{i=1}^{\sqrt{n}}(\dfrac{\sqrt{n}}{i})=\sqrt{n}\log n

代码：

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=1e6;
int n,mu[N+5],pri[N+5],ans;
void init(){
    mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!pri[i])pri[++pri[0]]=i,mu[i]=-1;
        for(int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j]))break;
            mu[i*pri[j]]=-mu[i];
        }
    }
}
signed main(){
    scanf("%lld",&n),init();
    for(int i=1;i*i<=n;i++)if(mu[i])for(int j=i;j*j<=n;j+=i)ans+=mu[i]*((min(j-1,(n/j)-j))/i);
    printf("%lld\n",ans*2+1);
    return 0;
}

VII.[NOI2016] 循环之美

依据小学数论知识，我们要求

因为后面的 k 是个常数，所以我们就想把它搞出来。

推到这步时，我们发现它好像跟我们最开头的式子莫名相似。

因此我们设原式为 f(n,m,k)，则我们这里后面的式子好像就是 f(m/d,n,d)。

因为每次都会除掉一个 d，而 n 和 m 是轮流颠倒着除的，所以复杂度应该有个 \log m\log n。又因为我们要暴力枚举 k 的所有因数，所以又乘上一个 \sqrt k。跑到 k=1 的层之后又要求 \sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=1]=\sum\limits_{d=1}^{\min(n,m)}\mu(d)\left\lfloor\dfrac{n}{d}\right\rfloor\left\lfloor\dfrac{m}{d}\right\rfloor，因此还要上杜教筛。

但实际上，这个算法复杂度不是很优，2s内卡过不去，权当参考。

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=20000000;
int pri[20010000],mu[20010000];
void sieve(){
    mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!pri[i])pri[++pri[0]]=i,mu[i]=-1;
        for(int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(i%pri[j])mu[i*pri[j]]=-mu[i];
            else break;
        }
    }
    for(int i=1;i<=N;i++)mu[i]+=mu[i-1];
}
unordered_map<int,int>mp;
int MU(int n){
    if(n<=N)return mu[n];
    if(mp.find(n)!=mp.end())return mp[n];
    int ret=1;
    for(int l=2,r;l<=n;l=r+1)r=n/(n/l),ret-=(r-l+1)*MU(n/l);
    return mp[n]=ret;
}
ll func(int n,int m,int k){
    if(!n||!m)return 0;
    if(k==1){
        ll ret=0;
        for(int l=1,r;l<=min(n,m);l=r+1)r=min(n/(n/l),m/(m/l)),ret+=1ll*(n/l)*(m/l)*(MU(r)-MU(l-1));
        return ret;
    }
    ll ret=0;
    for(int d=1;d*d<=k;d++){
        if(k%d)continue;
        ret+=func(m/d,n,d)*(mu[d]-mu[d-1]);
        if(d*d!=k)ret+=func(m/(k/d),n,k/d)*(mu[k/d]-mu[k/d-1]);
    }
    return ret;
}
int n,m,k;
int main(){
    sieve();
    scanf("%d%d%d",&n,&m,&k);
    printf("%lld\n",func(n,m,k));
    return 0;
}

接下来我们反其道而行之，保留 k，处理 i。

我们令 f(n) 表示 \sum\limits_{j=1}^n\Big[\gcd(j,k)=1\Big]。当 j\leq k 时，可以预处理；当 j>k 时，每 k 位一个循环节，且循环节内部大小为 \varphi(k)。

于是有式子

则现在 f 即可 O(1) 计算。我们即可将原式改写为 \sum\limits_{d=1}^{\min(n,m)}\mu(d)\left\lfloor\dfrac{n}{d}\right\rfloor\Big[\gcd(d,k)=1\Big]f\Big(\left\lfloor m/d\right\rfloor\Big)。

因为 \mu 是积性函数，且有平方项的数的 \mu 值为 0，故 \mu(dp)=\mu(d)\mu(p)\Big[\gcd(d,p)=1\Big]，于是上式改写成 \sum\limits_{p|k}\mu(p)\sum\limits_{d=1}^{\min(n,m)/p}\left\lfloor\dfrac{n}{dp}\right\rfloor f\Big(\left\lfloor m/dp\right\rfloor\Big)\mu(d)\mu(p)\Big[\gcd(d,p)=1\Big]。

发现了什么？后面一大坨跟我们开始的式子 \sum\limits_{d=1}^{\min(n,m)}\mu(d)\left\lfloor\dfrac{n}{d}\right\rfloor\Big[\gcd(d,k)=1\Big]f\Big(\left\lfloor m/d\right\rfloor\Big) 极度相似。于是我们设 g(n,m,k)=\sum\limits_{d=1}^{\min(n,m)}\mu(d)\left\lfloor\dfrac{n}{d}\right\rfloor\Big[\gcd(d,k)=1\Big]f\Big(\left\lfloor m/d\right\rfloor\Big)，则将最后的式子代进去，发现就是

边界就是到 k=1 时，有 g(n,m,1)=\sum\limits_{d=1}^{\min(n,m)}\mu(d)\left\lfloor\dfrac{n}{d}\right\rfloor f\Big(\left\lfloor m/d\right\rfloor\Big)。\mu 上杜教筛处理，剩下两个可以简单整除分块。

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=20000000;
int pri[20010000],mu[20010000];
void sieve(){
    mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!pri[i])pri[++pri[0]]=i,mu[i]=-1;
        for(int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
            pri[i*pri[j]]=true;
            if(i%pri[j])mu[i*pri[j]]=-mu[i];
            else break;
        }
    }
    for(int i=1;i<=N;i++)mu[i]+=mu[i-1];
}
unordered_map<int,int>mp;
int MU(int n){
    if(n<=N)return mu[n];
    if(mp.find(n)!=mp.end())return mp[n];
    int ret=1;
    for(int l=2,r;l<=n;l=r+1)r=n/(n/l),ret-=(r-l+1)*MU(n/l);
    return mp[n]=ret;
}
int F(int n); 
ll func(int n,int m,int k){
    if(!n||!m)return 0;
    if(k==1){
        ll ret=0;
        for(int l=1,r;l<=min(n,m);l=r+1)r=min(n/(n/l),m/(m/l)),ret+=1ll*(n/l)*F(m/l)*(MU(r)-MU(l-1));
        return ret;
    }
    ll ret=0;
    for(int d=1;d*d<=k;d++){
        if(k%d)continue;
        if(mu[d]!=mu[d-1])ret+=func(n/d,m/d,d);
        if(d*d!=k&&mu[k/d]!=mu[k/d-1])ret+=func(n/(k/d),m/(k/d),k/d);
    }
    return ret;
}
int k,f[2010];
int F(int n){return (n/k)*f[k]+f[n%k];}
int n,m;
int main(){
    sieve();
    scanf("%d%d%d",&n,&m,&k);
    for(int i=1;i<=k;i++)f[i]=f[i-1]+(__gcd(i,k)==1);
    printf("%lld\n",func(n,m,k));
    return 0;
}

本题是道难得的莫反好题，有多种思路，其中有三点我认为最为精妙：一是 [\gcd(i,j)=1][\gcd(j,k)=1]\Leftrightarrow[\gcd(j,ik)=1]，二是 \mu(dp)=\mu(d)\mu(p)[\gcd(d,p)=1]，这两点在乘积项和项的乘积间架起了桥梁；还有一点就是函数式思考，遇见形式相似的东西就尝试设函数递归处理，在某些时候不失为一种手段。

4.min25筛

听说这玩意能干杜教筛干不了的事？

同杜教筛一样，这也是用来求积性函数前缀和的东西。其复杂度为 O(\dfrac{n^{0.75}}{\log n})，大部分时候要略优于杜教筛。

min25筛作用的积性函数，应保证对于一切质数 p，f(p) 均是有关 p 的低阶多项式。同时，质数的正整数次幂的位置的 f 值要能够简单求出。

既然是多项式，就可以计算令 f(p)=p^k 时函数的前缀和，然后分别加在一起即可。

我们考虑将所有 i 按质数与合数分类：

接着，考虑枚举合数的最小质因子（该数一定 \leq\sqrt p）及其次数。

因为保证 f 是积性函数，所以可以变为

现在，我们考虑令 g(n,j)=\sum\limits_{i=1}^n[i\in\mathtt{prime}\lor i\text{的最小质因子大于第}j\text{个质数}]i^k，其中 k 是一常数。明显，i^k 这个函数是完全积性的，但它仅在质数处与我们强制令 f(p)=p^k 得到的函数值相等。

考虑由 g(n,j-1) 转移至 g(n,j)。显然，此时有些原本计入的 i 不能再被计入了。这样的 i 是所有最小质因子恰为第 j 个质数（设为 p_j）的合数。则我们显然可以提出一个 p_j 来。于是

其中，最后又减掉的那个 g 是为了避免重复计算前 j-1 个质数的影响（它们会使得最小质因子并非 p_j）

因为 i^k 是完全积性函数，所以里面套着的 g(n/p_j,j-1) 可能还有的 p_j 因子也可继续相乘。

因为 \left\lfloor\dfrac{\left\lfloor n/a\right\rfloor}{b}\right\rfloor=\left\lfloor\dfrac{n}{ab}\right\rfloor，所以对于某个 n，其所会访问到的状态一共只有所有的 \left\lfloor\dfrac{n}{x}\right\rfloor，共 \sqrt n 个。可以对这 \sqrt n 个数建立值与下标的双射，这样就压缩了 g 所需要的下标范围。

我们此时来设 h(n,j)=\sum\limits_{i=1}^n[i\text{的最小质因子大于第}j\text{个质数}]f(i)。则答案即为 h(n,0)。

考虑将 h(n,j) 划作两部分，即 >p_j 的质数的贡献以及最小质因子 >p_j 的合数的贡献。

质数的贡献是 g(n,|\mathtt{prime}|)-\sum\limits_{i=1}^{j}f(p_i)。

合数的贡献是 \sum\limits_{k=j+1}^{|\mathtt{prime}|}\sum\limits_{(p_k)^a\leq n}f\Big((p_k)^a\Big)\sum\limits_{i=1}^{\left\lfloor n/(p_k)^a\right\rfloor}[i\text{的最小质因子大于}p_k]f(i)。后面一坨的定义刚好就是 h\Big(n/(p_k)^a,k\Big)。于是可以化为 \sum\limits_{k=j+1}^{|\mathtt{prime}|}\sum\limits_{(p_k)^a\leq n}f\Big((p_k)^a\Big)\Bigg(h\Big(n/(p_k)^a,k\Big)+[a\neq1]\Bigg)。

注意到后面还有一个 [a\neq1]，因为我们在min25筛的时候只考虑了质数与合数，没有考虑 1，所以这里要补上；但是，当 a=1 时，如果考虑了 1 就会成为质数，已经在前面处理过了。这也意味着在min25筛执行完毕后，还要加上 1 的答案。

总递归式为

质数筛到 \sqrt n 就行了。

I.【模板】Min_25筛

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=1e9+7;
const int inv6=166666668;
ll n;
int m,pri[1001000],sp[2][1001000];//sp is the prefix sum of prime numbers' f, when 0 and 1 are the linear and quardratic terms respectively
void sieve(){
    for(int i=2;i<=m;i++){
        if(!pri[i])pri[++pri[0]]=i,sp[0][pri[0]]=(sp[0][pri[0]-1]+i)%mod,sp[1][pri[0]]=(sp[1][pri[0]-1]+1ll*i*i)%mod;
        for(int j=1;j<=pri[0]&&i*pri[j]<=m;j++){pri[i*pri[j]]=true;if(!(i%pri[j]))break;} 
    }
}
ll kth[1001000];//the kth of n/i
int tot,sml[1001000],lar[1001000];//sml is the bucket of elements no greater than m, while lar is for those greater than m.
int g[2][1001000];//the g array can DP directly.
int H(ll x,int y){
    if(pri[y]>=x)return 0;
    int X=(x<=m?sml[x]:lar[n/x]);
    int ret=(0ll+g[1][X]-g[0][X]+mod-(sp[1][y]-sp[0][y])+mod)%mod;
    for(int i=y+1;i<=pri[0]&&1ll*pri[i]*pri[i]<=x;i++){
        ll pa=pri[i];
        for(int a=1,tmp;pa<=x;a++,pa*=pri[i])tmp=pa%mod,(ret+=1ll*tmp*(tmp-1)%mod*(H(x/pa,i)+(a!=1))%mod)%=mod;
    }
//  printf("%d,%d:%d\n",x,y,ret);
    return ret;
}
int main(){
    scanf("%lld",&n),m=sqrt(n),sieve();
//  for(int i=1;i<=pri[0];i++)printf("%d ",pri[i]);puts("");
//  for(int i=1;i<=pri[0];i++)printf("%d ",sp[0][i]);puts("");
//  for(int i=1;i<=pri[0];i++)printf("%d ",sp[1][i]);puts("");
    for(ll i=1;i<=n;i=n/(n/i)+1){
        kth[++tot]=n/i;
        int I=kth[tot]%mod;//use the prefix sum formula of x and x^2 as the DP initial values of g(don't forget to ignore 1)
        g[0][tot]=(1ll*I*(I+1)/2+mod-1)%mod,g[1][tot]=(1ll*I*(I+1)%mod*(2*I+1)%mod*inv6+mod-1)%mod;
        if(kth[tot]<=m)sml[kth[tot]]=tot;else lar[n/kth[tot]]=tot;//because here we need to store values in buckets, so we can't use the moduloed values here.
    }
    for(int i=1;i<=pri[0];i++)for(int j=1;j<=tot&&1ll*pri[i]*pri[i]<=kth[j];j++){//cyclicize array g.
        ll nexi=kth[j]/pri[i];int I=(nexi<=m?sml[nexi]:lar[n/nexi]);//get the id of n/p[i]
        (g[0][j]+=mod-1ll*pri[i]*(g[0][I]+mod-sp[0][i-1])%mod)%=mod;
        (g[1][j]+=mod-1ll*pri[i]*pri[i]%mod*(g[1][I]+mod-sp[1][i-1])%mod)%=mod;
    }
    printf("%d\n",(H(n,0)+1)%mod);
    return 0;
}

II.LOJ#6053. 简单的函数

重申一下min25筛应用的条件：

• 是积性函数。
• 质数处取值是低阶多项式。
• 质数次幂处取值可以快速求出。

满足以上三点的任意函数均可以min25筛。

现在看到这题。乍一看 p\operatorname{xor}a 这种东西看上去一脸非多项式的样子；但是因为我们只关心质数处取值是否是多项式，因为除了 2 以外的质数均为奇，且质数处的 a 必为 1，则除了 f(2)=3 以外，其余的 f(p)=p-1，是低阶多项式。故我们只需特判掉 n\geq 2 时手动给筛出来的结果加上 2 即可。

代码：

#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
typedef long long ll;
ll n;
int m,pri[1001000],sp[2][1001000],g[2][1001000];
void sieve(){
    for(int i=2;i<=m;i++){
        if(!pri[i])pri[++pri[0]]=i,sp[0][pri[0]]=sp[0][pri[0]-1]+1,sp[1][pri[0]]=(sp[1][pri[0]-1]+i)%mod;
        for(int j=1;j<=pri[0]&&i*pri[j]<=m;j++){
            pri[i*pri[j]]=true;
            if(!(i%pri[j]))break;
        }
    }
}
ll kth[1001000];
int sml[1001000],lar[1001000],tot;
int H(ll x,int y){
    if(pri[y]>=x)return 0;
    int X=(x<=m?sml[x]:lar[n/x]);
    int ret=(2ll*mod+g[1][X]-g[0][X]-(sp[1][y]-sp[0][y]))%mod;
    for(int i=y+1;i<=pri[0]&&1ll*pri[i]*pri[i]<=x;i++){
        ll pa=pri[i];
        for(int a=1;pa<=x;a++,pa*=pri[i])(ret+=1ll*(pri[i]^a)*(H(x/pa,i)+(a!=1))%mod)%=mod;
    }
    return ret;
}
int main(){
    scanf("%lld",&n),m=sqrt(n),sieve();
//  for(int i=1;i<=pri[0];i++)printf("%d ",pri[i]);puts("");
//  for(int i=1;i<=pri[0];i++)printf("%d ",sp[0][i]);puts("");
//  for(int i=1;i<=pri[0];i++)printf("%d ",sp[1][i]);puts("");
    for(ll i=1;i<=n;i=n/(n/i)+1){
        kth[++tot]=n/i;
        int I=kth[tot]%mod;
        g[0][tot]=(I+mod-1)%mod,g[1][tot]=(1ll*I*(I+1)/2+mod-1)%mod;
        if(kth[tot]<=m)sml[kth[tot]]=tot;else lar[n/kth[tot]]=tot;
    }
//  for(int i=1;i<=tot;i++)printf("%d ",kth[i]);puts("");
    for(int i=1;i<=pri[0];i++)for(int j=1;j<=tot&&1ll*pri[i]*pri[i]<=kth[j];j++){
        ll nexi=kth[j]/pri[i];int I=(nexi<=m?sml[nexi]:lar[n/nexi]);
        (g[0][j]+=(mod-g[0][I]+sp[0][i-1])%mod)%=mod;
        (g[1][j]+=mod-1ll*pri[i]*(g[1][I]-sp[1][i-1]+mod)%mod)%=mod;
    }
    printf("%d\n",H(n,0)+1+2*(n>=2));
    return 0;
} 

III.UOJ#188. 【UR #13】Sanrd

题意：求 \sum\limits_{i=l}^rf(i)，其中 f(i) 为 i 的次大质因子。

显然其可以被转为两个前缀和相减的形式。

明显 f(i) 并非积性函数，所以常规min25筛处理不了。但是我们可以用非常规min25筛。

考虑仍令 h(n,j)=\sum\limits_{i=1}^n[i\text{的最小质因子大于第}j\text{个质数}]f(i)。

考虑仍然分质数和合数计算贡献。质数贡献为 0，故我们只需考虑合数。

考虑枚举最小质因子 p 及其在 i 中次数。则贡献为 \sum\limits_{k=j+1}^{|\mathtt{prime}|}\sum\limits_{(p_k)^a\leq n}h\Big(n/(p_k)^a,k+1\Big)+p_k\sum\limits_{i=p_k}^{\left\lfloor n/(p_k)^a\right\rfloor}[i\in\mathtt{prime}]。

其中，加号前面的是 f(i)\neq p_k 的答案，后面的是 f(i)=p_k 时的答案，此时必有且仅有一个 \geq p_k 的因数作为 i 的因子才会使 p_k 成为其次大质因子。

后面的东西可以用常规min25筛中的 g 简单得出。

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll n,ans;
int m,pri[1001000];
ll g[1001000],kth[1001000];
int tot,sml[1001000],lar[1001000];
void sieve(){
    memset(pri,0,sizeof(pri));
    for(int i=2;i<=m;i++){
        if(!pri[i])pri[++pri[0]]=i;
        for(int j=1;j<=pri[0]&&i*pri[j]<=m;j++){pri[i*pri[j]]=true;if(!(i%pri[j]))break;}
    }
}
ll H(ll x,int y){
    if(x<=pri[y])return 0;ll ret=0;
    for(int i=y+1;i<=pri[0]&&1ll*pri[i]*pri[i]<=x;i++)for(ll X=x/pri[i];X;X/=pri[i])ret+=H(X,i),ret+=(X>=pri[i]?pri[i]*(g[(X<=m?sml[X]:lar[n/X])]-i+1):0);
    return ret;
}
ll solve(){
    m=sqrt(n),sieve();
    tot=0,memset(kth,0,sizeof(kth)),memset(sml,0,sizeof(sml)),memset(lar,0,sizeof(lar)),memset(g,0,sizeof(g));
    for(ll i=1,I;i<=n;i=n/(n/i)+1)I=kth[++tot]=n/i,g[tot]=I-1,(I<=m?sml[I]:lar[n/I])=tot;
    for(int i=1;i<=pri[0];i++)for(int j=1;j<=tot&&1ll*pri[i]*pri[i]<=kth[j];j++){
        ll nexi=kth[j]/pri[i];int I=(nexi<=m?sml[nexi]:lar[n/nexi]);
        g[j]-=g[I]-(i-1);
    }
//  for(int j=1;j<=tot;j++)printf("%d ",kth[j]);puts("");
//  for(int j=1;j<=tot;j++)printf("%d ",g[j]);puts("");
    return H(n,0);
}
int main(){
    scanf("%lld",&n),n--,ans-=solve();
    scanf("%lld",&n),ans+=solve();
    printf("%lld\n",ans);
    return 0;
}

IV.LOJ#572. 「LibreOJ Round #11」Misaka Network 与求和

首先先考虑莫反推一波式子。

考虑整除分块掉 T，问题转换为求 f^k*\mu 的前缀和，其中 * 是狄利克雷卷积符号。

明显 f^k 不会是积性函数，故 f^k*\mu 亦不可能是。但我们的确可以类杜教筛地计算其前缀和。具体而言，即使 f^k 不是积性函数，但它仍然是数论函数，则狄利克雷卷积中什么单位元与逆元的理论仍然可以应用，因而在杜教筛中我们选取的函数是 1。于是就有