LG5325 模板Min_25筛

Posted 2020-11-27 autoint

Min_25筛

Min_25筛的目的是解决一类求和问题

[ sum_{i=1}^n f(i) ]

满足 (f) 是积性函数，且对于质数的幂 (f) 容易求出。比如 (f(p^k)=(p^k)^2)。

Min_25筛基于两个简单想法，模拟分解质因数和模拟埃氏筛。

模拟分解质因数

这一部分的想法是利用积性函数的性质，强行拆出最小质因数。再利用合数的最小质因数 (leq sqrt{n}) 的性质，强行分离合数和质数，来缩小最小质因数的枚举范围。

以下用 (operatorname{MPF}(x)) 函数表示 (x) 的最小质因数。

[ sum_{i=1}^n f(i)=1+sum_{i=2}^n f(i)=1+sum_{2 leq pleq sqrt{n},e geq 1,p^e leq n} f(p^e) left([e>1]+sum_{2 leq x leq lfloorfrac{n}{p^e} floor,operatorname{MPF}(x)>p}f(x) ight)+sum_{sqrt{n} < pleq n}f(p) ]

注意到 (sum_{i=2}^n f(i)) 和 (sum_{2 leq x leq lfloorfrac{n}{p^e} floor,operatorname{MPF}(x)>p}f(x)) 这两个求和式长得很像，只不过第二个求和式多了一个对最小质因数的限制。考虑用DP维护他们。

[ G(n,m)=sum_{2leq xleq n,operatorname{MPF}(x)>m} f(x)=sum_{m< pleqsqrt{n},egeq 1,p^eleq n} f(p^e)left([e>1]+sum_{2leq xleq lfloorfrac{n}{p^e} floor,operatorname{MPF}(x)>p}f(x) ight)+sum_{m< pleq n}f(p)=sum_{m< pleqsqrt{n},egeq 1,p^eleq n} f(p^e)([e>1]+G(lfloorfrac{n}{p^e} floor,p))+sum_{m< pleq n}f(p) ]

后面 (sum_{m< pleq n}f(p)) 这个针对质数的求和式是个经典问题。显然可以通过差分来计算。

[ H(n)=sum_{2leq pleq n}f(p)\sum_{m< pleq n}f(p)=H(n)-H(m) ]

而 (H) 的计算就是第二部分的内容。

模拟埃氏筛

这一部分的想法是仿造埃氏筛筛掉合数的方式从总和中减去合数的贡献。

所以我们令 (H′(i,j)) 表示埃式筛法枚举了前 (i) 个素数后，不超过 (j) 的所有还剩下的数（不包括 (1)）的平方和。

[ H'(i,j)=H'(i-1,j)-p_i^2(H'(i-1,lfloorfrac{j}{p_i} floor)-H'(i-1,p_i-1)) ]

这里因为 (f(p^k)) 是关于 (p^k) 的多项式，所以我们可以拆一个一次的质因数出来。

每一轮筛掉的数都是最小质因数为 (p_i) 且还含有另外一个 (geq p_i) 的质因数的数。所以如果 (j < p_i^2) 那么实际上没有筛掉任何数。

实现细节

[ H'(0,j)=sum_{i=1}^j i^2=frac{1}{6}j(j+1)(2j+1) ]

[ H'(i,j)=H'(i-1,j)-p_i^2(H'(i-1,lfloorfrac{j}{p_i} floor)-H'(i-1,p_i-1)) ]

[ H(n)=H'(max_{p_ileq sqrt{n}}{i},n) ]

[ G(n,max_{pleqsqrt{n}}{p})=H(n)-H(max_{pleqsqrt{n}}{p}) ]

[ G(n,m)=sum_{m< pleqsqrt{n},egeq 1,p^eleq n} f(p^e)([e>1]+G(lfloorfrac{n}{p^e} floor,p))+H(n)-H(m) ]

注意到 (G) 的第一维永远是 (lfloorfrac{n}{x} floor) 的形式，所以 (G) 和 (H') 都只用算 (sqrt{n}) 个值，并且可以用数组存储。

[ ans=1+G(n,0) ]

(G) 和 (H') 的空间复杂度都是 (O(sqrt{n})) 的。

我们要暴力把 (H) 在所有 (lfloorfrac{n}{x} floor) 处求出来，时间复杂度 (O(frac{n^{frac{3}{4}}}{log n}))。

如果直接递归求解 (G)，在 (nleq 10^{13}) 的情况下也是 (O(frac{n^{frac{3}{4}}}{log n})) 的。

LG5325 模板

(f(p^k)=(p^k)^2-p^k)，这时 (H) 可以通过两次求和算出。

代码参考：https://www.luogu.com.cn/blog/wucstdio/solution-p5325

CO int N=1e6+10;
int64 n;int S,pri[N],tot,sH1[N],sH2[N];
int64 w[N];int num,H1[N],H2[N],idx1[N],idx2[N];

int G(int64 x,int y){
    if(pri[y]>=x) return 0;
    int k=x<=S?idx1[x]:idx2[n/x];
    int ans=add(add(H2[k],mod-H1[k]),mod-add(sH2[y],mod-sH1[y]));
    for(int i=y+1;i<=tot and pri[i]<=x/pri[i];++i){
        int64 pe=pri[i];
        for(int e=1;pe<=x;++e,pe*=pri[i]){
            int f=mul(pe%mod,pe%mod-1);
            ans=add(ans,mul(f,G(x/pe,i)+(e>1)));
        }
    }
    return ans;
}

int main(){
    read(n),S=sqrt(n);
    for(int i=2;i<=S;++i){
        if(!pri[i]){
            pri[++tot]=i;
            sH1[tot]=add(sH1[tot-1],i);
            sH2[tot]=add(sH2[tot-1],mul(i,i));
        }
        for(int j=1;j<=tot and i*pri[j]<=S;++j){
            pri[i*pri[j]]=1;
            if(i%pri[j]==0) break;
        }
    }
    for(int64 l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        w[++num]=n/l;
        H1[num]=w[num]%mod;
        H2[num]=mul(H1[num],mul(H1[num]+1,mul(2*H1[num]+1,i6)));
        --H2[num];
        H1[num]=mul(H1[num],mul(H1[num]+1,i2));
        --H1[num];
        if(n/l<=S) idx1[n/l]=num;
        else idx2[n/(n/l)]=num;
    }
    for(int i=1;i<=tot;++i)for(int j=1;j<=num and pri[i]<=w[j]/pri[i];++j){
        int k=w[j]/pri[i]<=S?idx1[w[j]/pri[i]]:idx2[n/(w[j]/pri[i])];
        H1[j]=add(H1[j],mod-mul(pri[i],add(H1[k],mod-sH1[i-1])));
        H2[j]=add(H2[j],mod-mul(pri[i],mul(pri[i],add(H2[k],mod-sH2[i-1]))));
    }
    printf("%d
",G(n,0)+1);
    return 0;
}