【BIJI】MIN25筛笔记

虽然好早之前就学了，但是感觉每次要用到的时候都要搞忘orz】 并不准备讲太多原理。。

这是啥子

名称:min25筛 时间复杂度:$O(n^{\frac{3/4}{log2(n)} })$ 用处:求积性函数$f(x)$的前缀和(以及一些特殊的非积性函数的前缀和) 前提要求:对于质数$p$,$f(p)$是关于的p的多项式，$f(p^c)$可以比较方便地计算出来。

使用步骤

分为两步。

STEP1

求出对于任何$ x = \lfloor \frac{n}{i} \rfloor $ 的 $ \sum\limits_{i=2}^x f(i) * [i是质数] $ 我们先利用线性筛求出$ \sqrt{n}　$以内的所有质数之后，我们设$P_j$为从小到大第j个质数 设定$ g(n,j)=\sum_{i=1}^{n}[i \in P \ or\ \min(p)>P_j]f(i) $ 即n以内所有的f(i),这个i要满足是质数或者最小质因数$> P_j$。 关于g(n,0),就是所有数都当做质数来考虑的总和。 那么有如下的公式。 简单理解一下如果$P_j > n $，那么肯定无法筛掉任何数。 否则，我们考虑一下哪些会新被$P_j$新筛去。那么首先一定其最小质因数为$ P_j$,并且其次小质因数一定大于等于$P_j$。也就是我们删去的 $f(P_j)g(\frac{n}{P_j},j-1) $,然后我们由于会在$g(\frac{n}{P_j},j-1)$中多删去$P_1 $到$P_{j-1} $的质数的$f(x) $,那么减去就好了。 $$ g(n,j)= \begin{cases} g(n,j-1)&P_j^2\gt n \ g(n,j-1)-f(P_j)[g(\frac{n}{P_j},j-1)-\sum_{i=1}^{j-1}f(P_i)]&P_j^2\le n\end{cases} $$ 比较方便地我们在考虑的时候一般都从2开始(这个没有任何质数的1太奇葩了)，最后再加上f(1)即可

STEP2

对于step2分为两种。递归版本相对较快。而递推版本可以计算出对于任何$ x = \lfloor \frac{n}{i} \rfloor $ 的 $ \sum\limits_{i=1}^x f(i) $

递归版

其实就是我们在预处理完了之后再进行暴算。我们定义$S(n,j)=\sum_{i=2}^n[\min(p)\ge P_j]f(i)$。 那么有; $$ S(n,i)=g(n,|P|)-\sum_{j=1}^{i-1}f(P_j)+\sum_{k\ge i}^{|P|}\sum_{e}(f(P_k^e)S(\lfloor\frac{n}{P_k^e}\rfloor,k+1)+f(P_{k}^{e+1})) $$ 最后答案就是$S(n,1)$,直接递归下去搞就是了(没有啥初始化的问题)。 考虑一下就是所有的大于$P_{i-1}$的质数的$f(x)$之和，以及考虑暴力枚举所有最小质因数和最小质因数的幂次对应和对应的合数们。

递推版

这个可以算出所有$sum(n/1)$ , $sum(n/2) $,$sum(n/3) \dots $. 改变S定义我们设定$　S(i,j) $表示　　$\sum_{i=1}^{n}[i \in P \ or\ \min(p) \ge P_j]f(i) $看起来是不是和g的定义很像？ $$ S(n,j)=S(n,j+1)+\sum_{e|e>=1和p_j^{e+1}<=n}~f(p_j^e) * ( S(n/p_j^e,j+1)-sumf(j) )+f(p_j^{e+1}) $$ 这样我们在枚举递归的时候要从大往小枚举。 初始化$ S(n,MAX) $就是$ g(n,|P|) $ code: 递归版本 loj6063 简单的函数:

#include<stdio.h>
#include<cstdio>
#include<algorithm>
#include<cstdlib>
#include<iostream>
#include<cmath>
#define int long long 
using namespace std;
typedef int oo;

const int maxn = 2e5+5;
const int mod = 1e9+7;
bool mk[maxn];
int sq,w[maxn],pri[maxn],cnt,id1[maxn],id2[maxn],h[maxn],g[maxn];
int sm[maxn],tot,inv2;
int ksm(int a,int b) {
    int ans = 1; 
    for(;b;b>>=1,a=a*a%mod)
        if(b&1) ans = ans*a%mod;
    return ans;
}
//ｈ质数个数　ｇ质数和
void oula(int FF) {
    mk[1] = mk[0] = 1;
    for(int i=2;i<=FF;i++) {
        if(!mk[i]) { pri[++cnt]=i; sm[cnt]=(sm[cnt-1]+i)%mod; }
        for(int j=1;j<=cnt&&pri[j]*i<=FF;j++) {
            int k = pri[j]*i; mk[k] = 1;
            if(i%pri[j]==0) break;
        }
    }
}
int n;
int S(int x,int y) {
    if(x<=1||pri[y]>x) return 0;
    int k = (x<=sq)?id1[x]:id2[n/x];
    int ans = ( ( (g[k]-h[k])%mod -  ( sm[y-1] - (y-1) )%mod ) %mod + mod )%mod;
    if(y==1) ans+=2;
    for(int i=y;i<=cnt&&pri[i]*pri[i]<=x;i++) {
        int t1 = pri[i],t2=pri[i]*pri[i];
        for(int e=1;t2<=x;++e,t1=t2,t2*=pri[i]) {
            ans = ( ans + S(x/t1,i+1)*(pri[i]^e)%mod + (pri[i]^(e+1))%mod ) %mod;
        }
    }
    return ans;
}
main() {
    inv2 = ksm(2,mod-2);
    scanf("%lld",&n); sq = sqrt(n); oula(sq);
    for(int i=1,j;i<=n;i=j+1) {
        j = n/(n/i); w[++tot]=n/i;
        h[tot] = (w[tot]-1)%mod;
        g[tot] = ( w[tot]%mod*( (w[tot]+1)%mod )%mod*inv2%mod - 1)%mod;
        if(w[tot]<=sq) id1[w[tot]]=tot;
        else id2[j]=tot;
    }
    for(int j=1;j<=cnt;j++) { // up 
        for(int i=1;i<=tot&&pri[j]*pri[j]<=w[i];i++) { // down
            int k = (w[i]/pri[j]<=sq)?id1[w[i]/pri[j]]:id2[n/(w[i]/pri[j])];
            g[i] = ( g[i] - ( g[k] - sm[j-1] )%mod*pri[j]%mod )%mod;
            h[i] = ( h[i] - ( h[k] - (j-1) )%mod )%mod;
        }
    }
    printf("%lld", ( (S(n,1)+1)%mod + mod )%mod );
}

loj 572 Misaka Network 与求和 code:

#include<stdio.h>
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<map>
#include<cmath>
#define int unsigned int

using  namespace std;
const int maxn = 2e5;
int n,k,FF;
bool mk[maxn]; int pri[maxn],cnt;
int id1[maxn],id2[maxn],mu[maxn];
int gid(int x) {
    return (x<=FF)?id1[x]:id2[n/x];
}
int pk[maxn];
int ksm(int a,int b) {
    int ans=1;
    for(;b;b>>=1,a=a*a)
        if(b&1) ans=ans*a;
    return ans;
}
void oula() {
    mk[1]=mk[0]=1; mu[1]=1;
    for(int i=2;i<=FF;i++) {
        if(!mk[i]) { pri[++cnt]=i; mu[i]=-1; pk[cnt]=ksm(i,k); }
        for(int j=1;j<=cnt&&1ll*pri[j]*i<=FF;j++) {
            int k = pri[j]*i;
            mk[k]=1;
            if(i%pri[j]==0) {
                mu[k]=0; break;
            }
            mu[k]=-mu[i];
        }
    }
    for(int i=2;i<=FF;i++) mu[i]+=mu[i-1];
}
int w[maxn],tot;
int g[maxn],f[maxn];//pri cnt
void pre() {
    for(int i=1,j;i<=n;i=j+1) {
        w[++tot]=n/i; j=n/(w[tot]);
        if(w[tot]<=FF) id1[w[tot]]=tot;
        else id2[j]=tot;
        g[tot]=w[tot]-1;
    }
    for(int j=1;j<=cnt;j++) {
        int p=pri[j];
        for(int i=1;i<=tot&&1ll*p*p<=w[i];i++) {
            g[i]=g[i]-(g[gid(w[i]/p)]-j+1);
        }
    }
    for(int i=1;i<=tot;i++) f[i]=g[i];
    for(int j=cnt;j>=1;j--) {
        int p = pri[j];
        for(int i=1;p*p<=w[i]&&i<=tot;i++) {
            for(int k=1,now=p;1ll*now*p<=w[i];k++,now*=p) {
                int z = gid(w[i]/now);
                f[i] += (f[z]-g[z])+(g[z]-j)*pk[j]+pk[j];
            }
        }
    }
}
map<int,int>mm;
int gmu(int x) {
    if(x<=FF) return mu[x];
    if(mm.count(x)) return mm[x];
    int ans = 0;
    for(int i=2,j;i<=x;i=j+1) {
        j = (x/(x/i));
        ans += (j-i+1)*gmu(x/i);
    }
    return mm[x]=1-ans;
}
int erc(int x) {
    int ans = 0;
    int las = 0;
    for(int i=1,j;i<=x;i=j+1) {
        j = (x/(x/i));
        ans = ans + (gmu(j)-gmu(i-1))*(f[gid(x/i)]);
    }
    return ans;
}
main() {
    ios::sync_with_stdio(false);
    cin>>n>>k;
    FF = sqrt(n);
    oula(); pre();
    int ans = 0;
    int las = 0;
    for(int i=1,j;i<=n;i=j+1) {
        j = (n/(n/i));
        int now = erc(j);
        ans = ans + (n/i)*(n/i)*(now-las);
        las = now;
    }
    cout<<ans;
}