杜教筛入门

当学 Min 25 的一个前置知识。

---

算法内容。

定义 $S(n)=\sum _{i=1}^{n}f(i)$。对于一个函数 $g$，有：

$$\begin{array}{rl}\sum _{i=1}^n(f\ast g)(i)& =\sum _{i=1}^n\sum _{d|i}f(\frac{i}{d})g(d)\\ & =\sum _{d=1}^{n}g(d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }f(i)\\ & =\sum _{d=1}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )\\ ⟹g(1)S(n)& =\sum _{i=1}^n(f\ast g)(i)-\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )\\ ⟹S(n)& =\frac{\sum _{i=1}^n(f\ast g)(i)-\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )}{g(1)}\end{array}$$

所以如果存在函数 $g$，满足：

1. $g(1)\ne 0$
2. $\sum _{i=1}^n(f\ast g)(i)$ 可以快速计算
3. $g(i)$ 可以快速计算

可以通过记忆化搜索+数论分块快速计算 $S(n)$。可以用 unordered_map存储结果。

直接计算复杂度为 $O(n^{\frac{3}{4}})$。更好的计算是预处理前 $O(n^{\frac{2}{3}})$ 的 $S$ 值，可以做到 $O(n^{\frac{2}{3}})$ 的复杂度。

具体证明可以参考 link。证明本质是积分，简单的。

常用构建函数 $g$ 技巧：

1. $\sum _{d|n}^n\mu (d)=[n=1]$
2. $\sum _{d|n}^{n}u(d)\frac{n}{d}=\phi (d)$
3. $\sum _{d|n}^n\phi (d)=n$
4. $i^k·(\frac{n}{i}{)}^k=n^k$

例子：

---

$S(n)=\sum _{i=1}^n\mu (i)$。

令 $g(n)=1$，也即常函数。

则 $\sum _{i=1}^n(g\ast \mu )(i)=[n\ge 1]$

则

$$S(n)=[n\ge 1]-\sum _{i=2}^{n}S(\lfloor \frac{n}{i}\rfloor )$$

---

$S(n)=\sum _{i=1}^n\phi (i)$

令 $g(n)=1$

则 $\sum _{i=1}^n(g\ast \phi )(i)=\frac{n(n+1)}{2}$。

则

$$S(n)=\frac{n(n+1)}{2}-\sum _{i=2}^{n}S(\lfloor \frac{n}{i}\rfloor )$$

---

实操

计算 $\sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j)$

先化简：

$$\begin{array}{rl}\sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j)& =\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)=d]ijd\\ & =\sum _{d=1}^n\sum _{i=1}^{id\le n}\sum _{j=1}^{jd\le n}[gcd(i,j)=1]ij{d}^3\\ & =\sum _{d=1}^n{d}^3\sum _{i=1}^{id\le n}\sum _{j=1}^{jd\le n}\sum _{t|i,t|j}\mu (t)ij\\ & =\sum _{x=1}^n\sum _{t=1}^{xt\le n}{x}^3\mu (t)\sum _{i=1}^{i\le \lfloor \frac{n}{xt}\rfloor }\sum _{j=1}^{j\le \lfloor \frac{n}{xt}\rfloor }ij{t}^2\\ & =\sum _{x=1}^n\sum _{t=1}^{xt\le n}{x}^3\mu (t)t^2{\left(\frac{\lfloor \frac{n}{xt}\rfloor (\lfloor \frac{n}{xt}\rfloor +1)}{2}\right)}^2\\ & =\sum _{T=1}^n{T}^2\sum _{d|T}\mu (d)\frac{T}{d}{\left(\frac{\lfloor \frac{n}{T}\rfloor (\lfloor \frac{n}{T}\rfloor +1)}{2}\right)}^2\\ & =\sum _{T=1}^n{T}^2\phi (T){\left(\frac{\lfloor \frac{n}{T}\rfloor (\lfloor \frac{n}{T}\rfloor +1)}{2}\right)}^2\end{array}$$

令 $h(n)=\frac{n^2(n+1{)}^2}{4}$，$f(n)=n^2\phi (n)$。

对于上式，后者可以数论分块，问题化为求解 $f$ 的前缀和。

令 $g(n)=n^2$

则 $\sum _{i=1}^n(f\ast g)(i=\sum _{i=1}^n\sum _{d|i}{d}^2\phi (d)\ast \frac{i^2}{d^2}=\sum _{i=1}^n{i}^3=\frac{n^2(n+1{)}^2}{4}$

则 $g(1)=1$

则 $g(i)$ 可以快速计算。

所以 $S_f(n)=\frac{n^2(n+1{)}^2}{4}-\sum _{i=2}^n{i}^2{S}_f(\lfloor \frac{n}{i}\rfloor )$

数论分块即可。

模版代码一份。

#include<bits/stdc++.h>
using namespace std;
#define N 1050500
#define int long long 
const int it=1e6+7;
int v[N],pri[N],tot,p,phi[N],n,m,s[N],inv2,inv6,inv4;
int power(int a,int b){
    int ans=1;
    while(b){
        if(b&1)ans=ans*a%p;
        a=a*a%p;b>>=1;
    }
    return ans;
}
unordered_map<int,int>sit;
void init(){
    phi[1]=1;
    for(int i=2;i<it;i++){
        if(!v[i]){
            pri[++tot]=i;phi[i]=i-1;
        }
        for(int j=1;j<=tot&&i*pri[j]<it;++j){
            v[pri[j]*i]=1;
            if(i%pri[j]==0){
                phi[pri[j]*i]=pri[j]*phi[i];break;
            }
            phi[pri[j]*i]=phi[pri[j]]*phi[i];
        }
    }
    for(int i=1;i<it;++i)s[i]=(s[i-1]+i*i%p*phi[i]%p)%p;
}
int h(int n){
    n%=p;
    return n*n%p*(n+1)%p*(n+1)%p*inv4%p;
}
int pfs(int n){
    n%=p;
    return n*(n+1)%p*(n+n+1)%p*inv6%p;    
}
int calc(int n){
    if(n<it)return s[n];
    if(sit[n]!=0)return sit[n];
    int res=h(n);int lst=1,cur=0;
    for(int l=2,r;l<=n;l=r+1){
        r=min(n,n/(n/l));cur=pfs(r);
        res=(res+p-(cur-lst)%p*calc(n/l)%p)%p;lst=cur;
    }res=(res%p+p)%p;
    return sit[n]=res;
}
signed main(){
    ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
    cin>>p;cin>>n;inv2=power(2,p-2);inv4=power(4,p-2);inv6=power(6,p-2);
    init();
    int ans=0,lst=0;
    for(int l=1,r;l<=n;l=r+1){
        r=min(n,n/(n/l));int cur=calc(r);
        ans+=(cur-lst)*h(n/l)%p;ans%=p;lst=cur;
    }
    ans=(ans%p+p)%p;
    cout<<ans<<"\n";
}