Min25-筛

给定积性函数 \(f(x)\)，求 \(F(n)=\sum_{i=1}^nf(i)\)。

如果没有特殊说明，忽略掉所有函数在 \(1\) 处的取值。

第二部分

设 \(G(n)=\sum_{i=1}^n f(i)[i\in Prime]\)，假设我们现在对于所有 \(i\)，求出了 \(G(\left\lfloor\frac{n}{i}\right\rfloor)\) 的值，考虑如何计算 \(F(n)\)。

考虑设 \(F(i,j)\) 表示 \(\sum_{k=1}^j f(k)[\mathrm{min}_k\ge p_i]\)，其中 \(\mathrm{min}_k\) 表示 \(k\) 的最小质因子，而 \(p_i\) 表示第 \(i\) 个质数。

我们考虑如何计算 \(F(i,j)\)。

我们考虑枚举这个最小的质因子，设为 \(k\)，显然 \(k\ge i\)，接下来我们考虑计数满足 \(\mathrm{min}_j=p_k\) 的 \(j\)，显然有 \(p_k^2\le j\)，否则的话 \(j\) 一定有更小的质因子（我们先不考虑质数）。接下来要考虑的是 \(j\) 里面 \(p_k\) 出现了几次，不难发现答案应该是 \(F(k+1,\left\lfloor\frac{j}{p_k}\right\rfloor)+F(k+1,\left\lfloor\frac{j}{p_k^2}\right\rfloor)+\cdots\)。注意我们的 \(F(i,j)\) 里面是不考虑 \(1\) 的，所以 \(p_k^e\) 都没有被考虑，我们考虑把它们加上，这样我们有：\(\sum_{p_{k}^{e+1}\le j,e\ge 1} f(p_k^{e+1})+F(k+1,\frac{j}{p_k^e})\)。我们这里先不考虑 \(f(p_k)\)。不难发现我们还有一些包括 \(f_{p_k}\) 在内的一些质数的取值没有计算，所以我们需要加上它们，于是我们有：

\[F(i,j)=G(j)-\sum\limits_{k=1}^{i-1}f(p_k)+\sum\limits_{p_k^2\le j,k\ge i}\sum\limits_{p_k^{e+1}\le j,e\ge 1} f(p_k^{e+1})+F(k+1,\left\lfloor\frac{j}{p_k^e}\right\rfloor)*f(p_k^e) \]

综上，我们可以通过爆搜来计算 \(F(i,j)\)，最后答案就是 \(F(1,n)\)。

这段的复杂度并没有被证明，但是实测跑的并不慢。复杂度只被证明是 \(O(n^{1-\epsilon})\)

注意，上面的部分不需要记忆化，而且只能求出 \(F(n)\)，而不能求出所有的 \(F(\left\lfloor\frac{n}{i}\right\rfloor)\)。

由 \(\sum\) 的限制可以知道，我们预处理质数需要预处理 \(\sqrt{n}\) 以内的质数和第一个比 \(\sqrt{n}\) 大的质数，这是因为我们递归的停止条件是 \(p_k^2\le j\)，所以是需要最后一个质数满足 \(p^2>n\) 的。

第一部分

关键在于如何求出我们的 \(G\)。我们找一个完全积性函数 \(h\) 使得对于所有的 \(p\in Prime\)，有 \(h(p)=f(p)\)，且 \(h\) 便于计算前缀和。设 \(G(i,n)=\sum_{k=1}^n h(k)[k\in Prime \lor min_k> p_i]\)，这里我们递推考虑 \(p_i\) 筛掉了哪些数，就是那些最小质因子为 \(p_i\) 的合数。则有：

\[G(i,n)=G(i-1,n)-\sum\limits_{k=1}^n h(k)[k\notin Prime \land min_k=p_i] \]

假设 \(p_i^2\le n\)

考虑 \(G(i-1,\lfloor\frac{n}{p_i}\rfloor)\) 里面包含了哪些数，首先它包含 \(p_1,p_2,\cdots,p_{i-1}\)，其次，它包含所有最小质因子 \(\ge p_i\) 的数给它乘上一个 \(h(p_i)\)，那么我们所有最小质因子大于等于 \(p_i\) 的就变成了所有最小质因子为 \(p_i\) 的数，但是前面那部分 \(p_1p_i,\cdots,p_{i-1}p_{i}\) 我们是不要的，所以再把它们加回来，于是我们有：

\[G(i,n)=G(i-1,n)-G(i-1,\lfloor\frac{n}{p_i}\rfloor)h(p_i)+(\sum_{j=1}^{i-1}h(p_j))h(p_i) \]

那么最后的答案为 \(G(n)=G(m,n)\)，其中 \(m\) 可以设为 \(n\) 内的质数个数。

考虑这部分的复杂度，这部分实际上是一个可以进行 \(O(1)\) 转移的 DP，因此只需要考虑其状态数，第一维显然是 \(O(\sqrt n)\) 的，第二维需要满足 \(n\ge p_i^2\)，所以状态数等同于下面的式子：

\[\sum\limits_{i=1}^{\sqrt{n}}[i\in Prime]\sum\limits_{j,\not \exists j'<j,s.t.\lfloor\frac{n}{j'}\rfloor=\lfloor\frac{n}{j}\rfloor}[\lfloor\frac{n}{j}\rfloor\ge i^2] \]

注意第二个 \(\sum\)，\(\lfloor\frac{n}{j}\rfloor\) 相同的 \(j\) 我们只会计算一次。

其中，我们认为当 \(i\le n^\frac{1}{4}\) 的时候代价为 \(\sqrt{n}\)，其余情况认为是 \(\frac{n}{i^2}\)，前面直接加起来，后面用积分积起来，而质数的限制我们直接除以 \(\ln n\) 即可，这是因为 \(\pi(n)=\frac{n}{n\ln n}\)，等价于给质数每 \(\frac{1}{\ln n}\) 出现一次，而相邻 \(\frac{1}{\ln n}\) 个数以上函数取值差异有限，所以可以认为一段 \(\ln n\) 中的数函数值相等，于是我们就可以这样计算：

\[\frac{1}{\ln n}(\sum\limits_{i=1}^{n^{\frac{1}{4}}}\sqrt{n}+\int_{n^{\frac{1}{4}}}^{\sqrt n}\frac{n}{x^2}dx) \]

前面 \(\sum\) 显然是 \(n^{\frac{3}{4}}\)，后面是 \(n(n^{-\frac{1}{4}}-n^{-\frac{1}{2}})\)，可以认为是 \(n^{\frac{3}{4}}\)，所以总复杂度为 \(\frac{n^{\frac{3}{4}}}{\ln n}\)。

注意我们暂且认为 \(f(1)=h(1)=0\)

同时，第一部分启发我们，枚举所有数除了最大质因子以外的部分是很快的。

整个算法在数据范围是 \(10^{10}\) 时消耗的时间约为 \(2s\) 左右。

例题

洛谷模板题

链接

如果给定的函数 \(f\) 找不到对应的完全积性函数 \(h\) 满足要求的话，可以尝试把 \(f\) 拆成 \(a_1f_1+a_2f_2+\cdots a_nf_n\) 的形式，然后如果后者每个都可以找到对应的 \(h\)，然后对后面每个 \(f_i\) 都做一遍第一部分，然后就可以得到 \(G\)。

这提示我们，\(G\) 是可以拆开来计算的。

代码：

#include<bits/stdc++.h>
#define mset(a,b) memset((a),(b),sizeof((a)))
#define rep(i,l,r) for(int i=(l);i<=(r);i++)
#define dec(i,l,r) for(int i=(r);i>=(l);i--)
#define cmax(a,b) (((a)<(b))?(a=b):(a))
#define cmin(a,b) (((a)>(b))?(a=b):(a))
#define Next(k) for(int x=head[k];x;x=li[x].next)
#define vc vector
#define ar array
#define pi pair
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define N 1000010
#define M number
using namespace std;

typedef double dd;
typedef long double ld;
typedef long long ll;
typedef unsigned int uint;
typedef unsigned long long ull;
#define int long long
typedef pair<int,int> P;
typedef vector<int> vi;

const int INF=0x3f3f3f3f;
const dd eps=1e-9;
const int Len=100010;
const int mod=1e9+7;

template<typename T> inline void read(T &x) {
    x=0; int f=1;
    char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c == '-') f=-f;
    for(;isdigit(c);c=getchar()) x=x*10+c-'0';
    x*=f;
}

int pr[N],ta,h[N][2],sumh[N][2],id[N],idtot,n,g[N][2],inv6,inv2,G[N],pp[N][41],sumf[N];
bool no[N];
unordered_map<int,int> rk;

inline int ksm(int a,int b,int mod){
    int res=1;while(b){if(b&1)res=1ll*res*a%mod;a=1ll*a*a%mod;b>>=1;}return res;
}
inline int inv(int x){return ksm(x,mod-2,mod);}
inline int f(int p,int k){
    return pp[p][k]*(pp[p][k]-1)%mod;
}
inline void sieve(){
    no[1]=1;
    rep(i,2,Len){
        if(!no[i]) pr[++ta]=i;
        for(int j=1;j<=ta&&1ll*pr[j]*i<=Len;j++){
            no[pr[j]*i]=1;if(i%pr[j]==0) break;
        }
    }
    rep(i,1,ta){
        h[i][0]=pr[i]*pr[i]%mod;sumh[i][0]=(sumh[i-1][0]+h[i][0])%mod;
        h[i][1]=pr[i];sumh[i][1]=(sumh[i-1][1]+h[i][1])%mod;
    }
    rep(i,1,ta) pp[i][0]=1;
    rep(i,1,ta){
        rep(j,1,40) pp[i][j]=pp[i][j-1]*pr[i]%mod;
    }
    rep(i,1,ta){
        sumf[i]=(sumf[i-1]+f(i,1))%mod;
    }
}
inline int SumH(int n,int op){
    n%=mod;
    if(op==0) return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
    else return n*(n+1)%mod*inv2%mod;
}
inline void Init(){
    sieve();
    for(int l=1,r;l<=n;l=r+1){
        r=n/(n/l);id[++idtot]=n/l;
    }
    reverse(id+1,id+idtot+1);
    rep(i,1,idtot) rk[id[i]]=i;
    inv6=inv(6);inv2=inv(2);
}
inline void GetG(int op){
    rep(i,1,idtot) g[i][op]=SumH(id[i],op)-1;
    rep(i,1,ta){
        dec(j,1,idtot){
            if(pr[i]*pr[i]>id[j]) break;
            // if(pr[i-1]*pr[i-1]>id[j]/pr[i]) break;
            g[j][op]=((g[j][op]-g[rk[id[j]/pr[i]]][op]*h[i][op]%mod+sumh[i-1][op]*h[i][op]%mod)%mod+mod)%mod;
        }
    }
}
inline int dfs(int i,int j){
    if(i>ta) return 0;
    if(pr[i]>id[j]) return 0;
    int ans=((G[j]-sumf[i-1])%mod+mod)%mod;
    for(int k=i;k<=ta&&pr[k]*pr[k]<=id[j];k++){
        int now=pr[k];
        for(int e=1;now<=id[j]/pr[k];e++,now*=pr[k]){
            ans=((ans+f(k,e+1))%mod+dfs(k+1,rk[id[j]/(now)])*f(k,e)%mod)%mod;
        }
    }
    return ans;
}
inline int GetF(){return dfs(1,idtot);}

signed main(){
    // freopen("my.in","r",stdin);
    // freopen("std.out","w",stdout);
    read(n);Init();GetG(0);GetG(1);
    rep(i,1,idtot) G[i]=((g[i][0]-g[i][1])%mod+mod)%mod;
    int Ans=GetF();
    printf("%lld\n",((Ans+1)%mod+mod)%mod);
    return 0;
}

犯的错误：

• 在预处理 \(G\) 的时候注意不要加上限制 \(p_{i-1}^2\le \lfloor\frac{j}{p_i}\rfloor\)，这是因为即使是 \(p_i^2<j\) 的时候，\(G(i,j)\) 也是有意义的，只不过其没必要再继续计算，直接从 \(G(i-1,j)\) 传递过来即可。

求莫比乌斯函数前缀和

\[G(n)=\sum\limits_{i=1}^n-[i\in Prime]=-\sum\limits_{i=1}^n [i\in Prime] \]

取 \(h=I\) 即可。

欧拉函数前缀和

\[G(n)=\sum\limits_{i=1}^n(i-1)[i\in Prime]=\sum\limits_{i=1}^ni[i\in Prime]+\sum\limits_{i=1}^n 1[i\in Prime] \]

取 \(h_1(x)=x,h_2=I\) 即可。

LOJ6053

有 \(f(p)=p\oplus1\)，实际上除了 \(p=2\) 的情况其它情况都是 \(f(p)=p-1\)，我们暂且认为 \(f(p)=p-1\)，这可以直接用上面的做法做，然后特殊处理 \(2\)，即把所有 \(G\) 的位置加上 \(2\)。然后做第二部分即可。

代码：

int pr[N],ta,h[N][2],sumh[N][2],id[N],idtot,n,g[N][2],inv6,inv2,G[N],sumf[N];
bool no[N];
unordered_map<int,int> rk;

inline int ksm(int a,int b,int mod){
    int res=1;while(b){if(b&1)res=1ll*res*a%mod;a=1ll*a*a%mod;b>>=1;}return res;
}
inline int inv(int x){return ksm(x,mod-2,mod);}
inline int f(int p,int k){return pr[p]^k;}
inline void sieve(){
    no[1]=1;
    rep(i,2,Len){
        if(!no[i]) pr[++ta]=i;
        for(int j=1;j<=ta&&1ll*pr[j]*i<=Len;j++){
            no[pr[j]*i]=1;if(i%pr[j]==0) break;
        }
    }
    rep(i,1,ta){
        h[i][0]=pr[i]%mod;sumh[i][0]=(sumh[i-1][0]+h[i][0])%mod;
        h[i][1]=1;sumh[i][1]=(sumh[i-1][1]+h[i][1])%mod;
    }
    rep(i,1,ta){
        sumf[i]=(sumf[i-1]+f(i,1))%mod;
    }
}
inline int SumH(int n,int op){
    n%=mod;
    if(op==0) return n*(n+1)%mod*inv2%mod;
    else return n;
}
inline void Init(){
    sieve();
    for(int l=1,r;l<=n;l=r+1){
        r=n/(n/l);id[++idtot]=n/l;
    }
    reverse(id+1,id+idtot+1);
    rep(i,1,idtot) rk[id[i]]=i;
    inv6=inv(6);inv2=inv(2);
}
inline void GetG(int op){
    rep(i,1,idtot) g[i][op]=SumH(id[i],op)-1;
    rep(i,1,ta){
        dec(j,1,idtot){
            if(pr[i]*pr[i]>id[j]) break;
            // if(pr[i-1]*pr[i-1]>id[j]/pr[i]) break;
            g[j][op]=((g[j][op]-g[rk[id[j]/pr[i]]][op]*h[i][op]%mod+sumh[i-1][op]*h[i][op]%mod)%mod+mod)%mod;
        }
    }
}
inline int dfs(int i,int j){
    if(i>ta) return 0;
    if(pr[i]>id[j]) return 0;
    int ans=((G[j]-sumf[i-1])%mod+mod)%mod;
    for(int k=i;k<=ta&&pr[k]*pr[k]<=id[j];k++){
        int now=pr[k];
        for(int e=1;now<=id[j]/pr[k];e++,now*=pr[k]){
            ans=((ans+f(k,e+1))%mod+dfs(k+1,rk[id[j]/(now)])*f(k,e)%mod)%mod;
        }
    }
    return ans;
}
inline int GetF(){return dfs(1,idtot);}

signed main(){
    // freopen("my.in","r",stdin);
    // freopen("std.out","w",stdout);
    read(n);Init();GetG(0);GetG(1);
    rep(i,1,idtot) G[i]=((g[i][0]-g[i][1])%mod+mod)%mod;
    rep(i,2,idtot) G[i]=(G[i]+2)%mod;
    // rep(i,1,idtot) printf("G[%d]=%d\n",i,G[i]);
    int Ans=GetF();
    printf("%lld\n",((Ans+1)%mod+mod)%mod);
    return 0;
}