求和技巧（min_25筛和被吊打的杜教筛）

吉老师的博客

第一类：因为 $g$ 是完全积性，故

((f\cdot g)*g)(n)=\sum_{d|n}f(d)g(d)g(\frac{n}{d})=g(n)\sum_{d|n}f(d)=(g\cdot(f*1))(n)

于是套了一发杜教筛的公式，相当于给原函数卷了一个 $g$ 。
若 $g\cdot(f*1)$ 非常优秀且 $g$ 是完全积性就可以套公式啦（比如求 $\sum \varphi\cdot id^k$ ）
第二类：

\sum_{i=1}^n(f*g*1)(i)=\sum_{i=1}^n\sum_{d|i}g(d)(f*1)(\frac{i}{d})=\sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}(f*1)(i)

故若我们能快速求 $f*1$ 的前缀和及 $g$ 的点值我们就能套杜教筛了，相当于给原函数卷了一个 $1$ 。

这里推式子的方法确实是条条大路通罗马，一般用 $\varphi$ 和用 $\mu$ 都能成，用 $\varphi$ 推的一般可能会用一些骚操作，而用 $\mu$ 的却往往可以非常暴力顺畅地推...

这个是复杂度分析，这种分析方式很重要，不仅限于杜教筛。有时有一些除法分块套除法分块直接做是 $O(n^{\frac{3}{4}})$ 的，也可以通过线性筛一部分的方式得到 $O(n^{\frac{2}{3}})$ ..
min25筛
这个算法就是把“一个数最多只有一个大于 $\sqrt n$ 的质因子”和“ $\lfloor\frac{n}{x}\rfloor$ 的值个数 $O(\sqrt n)$ ”运用到了极致。

考虑把所有数分成质数和合数。

第一步要求质数的函数值求和。
考虑筛法，拿质数 $p_i$ 去筛掉mindiv为 $p_i$ 的数。因为 $n$ 内不存在一个数的 mindiv 大于 $\sqrt n$ ，于是只要用小于等于 $\sqrt n$ 的质数去筛就好了。

第二步是合数，对于合数枚举它的mindiv和它的指数，即 $p^e$ ，再除掉递归求解就ok。

实际上min-25筛和洲阁筛都是对埃氏筛的扩展。

#include<bits/stdc++.h>
#define ULL unsigned long long
using namespace std;
const int N=2e6+50;
int t,sq,p[N],tot,sp[N],tp;ULL a[N],mx,n,k,d0[N],st[N],b[N];bool mark[N];
int getid(ULL x){return x<=sq?x:sq+n/x;}//经典trick求id
ULL S(ULL n,int y){//函数意思是2到n mindiv>p[y]的函数值和
    if(n<=p[y])return 0;
    int id=getid(n);ULL ret=(k+1)*(d0[id]-sp[y]);//这里是质数的函数值
    for(int i=y+1;i<=tot&&1ll*p[i]*p[i]<=n;i++){//合数的mindiv一定<=sqrt(n)。这里和下面的剪枝保证了复杂度，一定要加上。
        for(ULL j=p[i],d=k+1;j<=n;j*=p[i],d+=k)
            ret+=d*(S(n/j,i)+(j!=p[i]));//p^k本身也是合数，也要加上。
    }
    return ret;
}
void solve(){
    sq=sqrt(n);tp=0;
    for(ULL l=1,r;l<=n;l=r+1){
        ULL d=n/l;r=n/d;
        st[++tp]=d;d0[getid(d)]=d-1;//先把1减掉，之后好处理。（否则会出现本来想构造合数但是一个质数*1还是质数的情况）
    }
    for(int i=1;i<=tot;i++){
        ULL sq=1ll*p[i]*p[i];if(sq>n)break;
        for(int j=1;st[j]>=sq;j++){//st数组里的数是递减的。被筛掉的数一定>=sq。
            int id=getid(st[j]),iid=getid(st[j]/p[i]);
            d0[id]-=d0[iid]-sp[i-1];//这里减去筛掉的部分，因为数组的值还包括已经是质数的，所以还要减掉。
        }
    }
    printf("%llu\n",S(n,0)+1);//最后把1加上
}
int main(){
    scanf("%d",&t);
    for(int i=1;i<=t;i++)scanf("%llu%llu",&a[i],&b[i]),mx=max(mx,a[i]);
    sq=sqrt(mx);
    for(int i=2;i<=sq;i++){//这里是线性筛预处理。
        if(!mark[i])p[++tot]=i;
        for(int j=1,y;j<=tot&&(y=p[j]*i)<=sq&&(mark[y]=1);j++)
            if(i%p[j]==0)break;
    }
    for(int i=1;i<=tot;i++)sp[i]=i;//质数的函数值前缀和
    for(int i=1;i<=t;i++)n=a[i],k=b[i],solve();
    return 0;
}

关于它的复杂度，除了搜索以外的部分的复杂度是 $n^{\frac{3}{4}}{\log n}$ 的，搜索的部分如果写成记忆化或者dp也是这个复杂度，因为它们基本是同构的..然而根据玄学分析发现不记忆化时虽然理论复杂度是很差的 $O(\frac{n}{\log \log n})$ ，但在它可以用的范围 $n\leq 10^13$ 内，运算次数实际上也是 $n^{\frac{3}{4}}{\log n}$ 的。而实际表现更优秀。关于min25筛和洲阁筛的复杂度分析也很重要，有时虽然不是数论函数求和，但也要用到这种dp。

zzt论文中的这部分复杂度分析引用了这篇博客

而min25筛也可以进一步优化..

不过比较神仙。这个被人戏称为min26筛..
min25筛非常好用。只要发现函数值是一个积性函数， $p^k$ 处的函数值容易计算， $p$ 处的函数值是低次多项式即可快速计算。

刚刚的代码是求 $\sum_{i=1}^n\sigma_0(i^k)$ 。

当 $k=2$ 时有一个疯狂推式子的做法。简单提一提

对 $n^2$ 约数的每个质因子进行观察可以发现它们的次数要么是奇数要么是偶数（。。。），可以考虑构造映射，对奇数次数的质因子的次数都加1，就变成了一个平方数 $d^2$ 。 $d^2$ 一定是 $n^2$ 的约数，于是 $d$ 是 $n$ 的约数。所以我们可以枚举 $n$ 的所有约数 $d$ ，再对 $d$ 的所有不同的质因子选择是否给它的次数 $-1$ 。我们有了一个式子

\sigma_0(n^2)=\sum_{d|n}2^{w(d)}

其中 $w(n)$ 表示 $n$ 的不同质因子个数。
再从刚才提到的“对 $d$ 的所有不同的质因子选择是否给它的次数 $-1$ ”可以发现我们相当于是对 $n$ 的所有没有相同质因子的因数计了个数，即

2^{w(n)}=\sum_{d|n}\mu^2(d)

我们现在求前缀和。

\sum_{i=1}^n\sigma_0(i^2)=\sum_{d=1}^n2^{w(d)}\lfloor\frac{n}{d}\rfloor

我们除法分块以后要求 $2^{w(n)}$ 的前缀和。

\sum_{i=1}^n2^{w(i)}=\sum_{d=1}^n\mu^2(d)\lfloor\frac{n}{d}\rfloor

又要求 $\mu^2(n)$ 的前缀和，而这是一个经典的容斥。

\sum_{i=1}^n\mu^2(i)=\sum_{i=1}^{\sqrt n}\mu(i)\lfloor\frac{n}{i^2}\rfloor

发现这三个式子都是根号计算，于是统统筛出 $n^{\frac{2}{3}}$ 的值之后暴力计算时空复杂度都为 $O(n^{\frac{2}{3}})$ 。这题太草了。。还卡空间。
而用min25筛呢？发现要求的函数显然是个积性函数，在 p 处的函数值为 $3$ ，在 $p^e$ 处函数值为 $2e+1$ 。直接套就AC了。。。min_25nb！

还要注意min25其实不要局限于求积性函数前缀和。它可以魔改的。注意我们搜索的部分其实已经得到了所有素数幂的信息，然后我们在适当的时候统计上贡献就可以了。比如这么一个函数的前缀和： $f(x)^k$ ，其中 $f(x)$ 表示它的第二大质因子，可重，比如 $f(4)=2,f(12)=2,f(p)=1,f(1)=0$ 。这时注意到素数的贡献是 $1$ ，合数的次大质因子也比 $\sqrt n$ 小，这样我们可以在搜索时统计一下多少数的次大值因子是某个质数，累加一下贡献即可。

而如果我们要用min25筛求出所有 $\sqrt n$ 个位置的值，那么这个搜索就不太有用了...也许我们会想记忆化，但这个空间好像也有点承受不住...

于是我们只能写dp的做法了。dp部分的状态仍然是这个状态。我们从大到小枚举质数，每次加上 mindiv为这个质数的数的贡献。不过注意我们复杂度优化的关键在 $x\leq p^2$ 的时候，我们就不能访问 $x$ 了。否则复杂度会退化为 $\frac{n}{\log n}$ 。当 $p\leq x$ 时 $p$ 对 $x$ 就有贡献，这有一点麻烦。不过注意到这时只有质数的贡献，我们已经把质数的贡献预处理了。于是为了方便递推，我们可以把状态改一下，改成只有合数的贡献。质数的贡献我们现算即可。其实发现第一部分的dp也是这个思路。
而复杂度分析我们只要注意到 $\sum\limits_{p^k\leq n}1$ 是 $O(\frac{n}{\ln n})$ 就能得到它和前一部分的复杂度是相同的...其实这种式子中把 $p$ 换成 $p^k$ 大部分时候是不影响它的量级的...
我们发现它确实跑得比搜索慢，因为搜索访问不到所有的状态。

#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=1e5+50;
int T,n,phi[N],p[N],tot,st[N],tp,sq,d0[N],tt,cnt1,cnt2;bool vis[N];LL sp[N],ans[N],d1[N];
int getid(int d){return d<=sq?d:sq+n/d;}
LL S(int n,int f){
    if(n<p[f])return 0;
    int id=getid(n);LL ret=d1[id]-sp[f-1]-(d0[id]-(f-1));
    for(int i=f;1ll*p[i]*p[i]<=n;i++)
        for(LL j=p[i],k=p[i]-1;j<=n;j*=p[i],k*=p[i])
            cnt2++,ret+=k*(S(n/j,i+1)+(j!=p[i]));
    return ret;
}
void solve(){
    scanf("%d",&n);tp=0;sq=sqrt(n);
    for(LL l=1,r,d,id;l<=n;l=r+1){
        d=n/l;r=n/d;st[++tp]=d;
        id=getid(d);d0[id]=d-1;
        d1[id]=1ll*d*(d+1)/2-1;
    }
    for(int i=1;i<=tot;i++){
        if(1ll*p[i]*p[i]>n)break;tt=i;
        int sq=p[i]*p[i];
        for(int j=1;sq<=st[j];j++){
            int id=getid(st[j]),iid=getid(st[j]/p[i]);
            d0[id]-=d0[iid]-(i-1);
            d1[id]-=1ll*p[i]*(d1[iid]-sp[i-1]);
        }
    }
    for(int i=tt;i;i--){
        int sq=p[i]*p[i];
        for(int j=1;sq<=st[j];j++){
            int m=st[j],id=getid(m);
            for(LL k=p[i],v=p[i]-1,iid;k<=m;k*=p[i],v*=p[i]){
                iid=getid(m/k);cnt1++;
                ans[id]+=v*(ans[iid]+max(0ll,d1[iid]-sp[i]-(d0[iid]-i))+(k!=p[i]));
            }
        }
    }
    int id=getid(n);
    cout<<ans[id]+d1[id]-d0[id]+1<<" "<<S(n,1)+1<<endl;
    cout<<cnt1<<" "<<cnt2<<endl;
}
int main(){
    for(int i=2;i<N;i++){
        if(!vis[i])p[++tot]=i,phi[i]=i-1;
        for(int j=1,y;j<=tot&&(y=p[j]*i)<N;j++){
            vis[y]=1;
            if(i%p[j])phi[y]=phi[i]*(p[j]-1);
            else{phi[y]=phi[i]*p[j];break;}
        }
    }
    for(int i=1;i<=tot;i++)sp[i]=sp[i-1]+p[i];
    scanf("%d",&T);
    while(T--)solve();
    return 0;
}
// 1
// 2147483647

powerful number

这个技巧相当于是杜教筛的一个扩展。我们要求一个积性函数 $F$ 的前缀和，它可能长得比较奇怪，不过没有关系。我们观察它在质数处的取值，找到一个积性函数 $G$ 满足在质数处（也可以有常数个值例外，这样不影响复杂度）取值与 $F$ 相同，且容易求前缀和。那么函数 $H=F/G$ 在质数处的取值就为 $0$ 。其中 $/$ 为狄利克雷除法，即卷积的逆运算。那么

\sum\limits_{i=1}^n F(i)=\sum\limits_{i=1}^n\sum\limits_{d|i}H(d)G(\frac{i}{d})=\sum\limits_{d=1}^nH(d)\sum\limits_{i=1}^{\lfloor\frac{n}{i}\rfloor}G(i)

而因为 $H$ 在质数处为 $0$ 且是积性函数，那么只有每个质因子次数都大于 $1$ 的数的 $H$ 才可能有值，我们把这种数叫做powerful number。它一定能被表示成 $a^2b^3$ ，其中 $b$ 没有平方因子。那么 $n$ 以内这种数的个数不超过 $\sum\limits_{i=1}^{\sqrt n}(n/i^2)^{1/3}$ ，积分一下可以得到是 $O(\sqrt n)$ 的。于是我们可以爆搜出所有的powerful number，算出对应的 $H$ ，再算一下 $G$ 的前缀和（且只需要所有 $\lfloor\frac{n}{x}\rfloor$ 处的前缀和，所以也可以杜教筛等来处理），就可以得到答案了。
下面这个代码是 $f(1)=1,f(p^e)=p$ xor $e$ 的积性函数的前缀和，发现 $p>2$ 时 $f(p)=p-1$ ，那么用 $\varphi$ 即可。注意杜教筛线性筛的部分略比 $n^{\frac{2}{3}}$ 小一点可以跑得更快，也更省空间。这题也可以用min25筛做，也是注意到 $f(p)=p-1$ ，是多项式，特判一下 $2$ 即可。跑得比这个还要快。

#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=2e6,mod=1e9+7,M=1e5,P=1e4;//不超过1e5的质数只有1e4个
int p[N/10],tot,phi[N],sq,val[M],ans,h[P][40];bool vis[N];LL n;
int S(LL n){
    LL x=n&1?n:n/2,y=n&1?(n+1)/2:(n+1);
    return 1ll*(x%mod)*(y%mod)%mod;
}
void dfs(int x,LL v,int d){
    ans=(ans+1ll*d*(n/v<N?phi[n/v]:val[n/(n/v)]))%mod;
    for(int i=x;i<=tot;i++){
        if(i>1&&v*p[i]*p[i]>n)return;assert(i<P);
        for(LL j=1,k=p[i];k*v<=n;j++,k*=p[i]){
            if(h[i][j]==-1){
                int H=p[i]^j,G=p[i]-1;//因为这里值可能很大，我们并没有筛到，所以现算
                for(int l=1;l<=j;l++,G=1ll*G*p[i]%mod)
                    H=(H-1ll*h[i][j-l]*G)%mod;//求狄利克雷除法，因为是积性函数所以只在质数幂求即可
                h[i][j]=H;
            }
            if(h[i][j])dfs(i+1,v*k,1ll*d*h[i][j]%mod);
        }
    }
}
int main(){//loj6053
    scanf("%lld",&n);sq=sqrt(n);phi[1]=1;
    for(int i=2;i<N;i++){
        if(!vis[i])p[++tot]=i,phi[i]=i-1;
        for(int j=1,y;j<=tot&&(y=p[j]*i)<N;j++){
            vis[y]=1;
            if(i%p[j])phi[y]=phi[i]*(p[j]-1);
            else {phi[y]=phi[i]*p[j];break;}
        }
    }
    for(int i=1;i<N;i++)(phi[i]+=phi[i-1])%=mod;
    int s=1;while(n/s>=N)s++;
    for(int i=s;i;i--){
        LL m=n/i,v=S(m);
        for(LL l=2,r,d;l<=m;l=r+1){
            d=m/l;r=m/d;
            v=(v-(r-l+1)%mod*(d<N?phi[d]:val[i*r]))%mod;//这里推一下n/(n/i/d)可以发现是i*r。
        }
        val[i]=v;
    }
    memset(h,-1,sizeof(h));
    for(int i=1;i<P;i++)h[i][0]=1;
    dfs(1,1,1);cout<<(ans+mod)%mod;
    return 0;
}

而在冬令营课件中有一个运用了该技巧的特别特别神仙的 $O(\frac{n^{\frac{2}{3}}}{\log n})$ 对积性函数求和的技巧......这辈子都不会碰的。

Prutekoi's blog

求和技巧（min_25筛和被吊打的杜教筛）

powerful number

【清华集训2015】V