• [学习笔记] min_25 筛


    一些约定

    1. \(p_k\) 表示第 \(k\) 个质数,特别地,\(p_0=1\)
    2. \(n/k\)\(\left \lfloor \frac{n}{k}\right \rfloor\)
    3. \(\text{lpf}(x)\) 为数 \(x\) 的最小质因数。

    要干啥?

    积性函数 \(F(i)\) 的前缀和。

    推柿子

    \(h(i,j)=\sum_{x=2}^i[\text{lpf}(x)> p_j]F(x)\),那么答案就是 \(h(n,0)+F(1)=h(n,0)+1\)(积性函数的第一项一定为 \(1\)\(1\) 和任何数互质嘛)。求解的普遍思路似乎是递推,不妨先按照定义展开(注意,本次展开只包含可被分解成多个质数的合数):

    \[h(i,j)=\sum_{k=j+1}\sum_{e=1} F(p_k^e)\cdot h(i/p_{k}^e,k) \]

    其实就是枚举最小质因子及其幂次,然后进行递归(易得它们是互质的)。

    另外就还剩下质数、单个质数的幂的情况,后者也是比较容易写的:

    \[h(i,j)=\sum_{k=j+1}\sum_{e=1} F(p_k^e)\cdot h(i/p_{k}^e,k)+F(p_k^{e+1})\ \ \ \ (p_{k}^{e+1}\le i) \]

    需要注意后面的限制,可以发现,当 \(p_k^{e+1}>i\) 时,\(F(p_k^e)\cdot h(i/p_{k}^e,k)\) 一定也为零,所以这个限制是紧凑且合理的。

    \(P(n)\) 为小于等于 \(n\) 的质数的 \(F\) 值的前缀和,那么最终就有(注意到质数只需要枚举到 \(\sqrt n\) 级别):

    \[h(i,j)=P(i)-P(p_j)+\sum_{k=j+1}\sum_{e=1} F(p_k^e)\cdot h(i/p_{k}^e,k)+F(p_k^{e+1})\ \ \ \ (p_{k}^{e+1}\le i) \]

    于是我们只用计算 \(P(n)\).

    首先,我们需要 \(F(n)\) 是一个 次数较小 的多项式,这样拆开的每一项就可以视作一个 完全积性函数 —— 更一般地,不妨设 \(f_k(n)=n^k\),只需要计算这个函数的 \(P_k(n)\),乘上系数再相加即可得到 \(P(n)\).

    \(g(i,j)\) 为在 \([2,i]\) 之间且 \(\text{lpf}>p_j\) 的数字的 \(f\) 值之和,特别地,这一坨恒包含小于等于 \(i\) 的质数的 \(f\) 值。

    考虑递推地求取 \(g(i,j)\) —— 我们可以将 \(g(i,j)\) 理解为在 \(g(i,j-1)\) 的基础上,去掉 最小质因子\(p_j\)合数\(f\) 值:

    \[g(i,j)=g(i,j-1)-f_k(p_j)\cdot\left (g(i/p_j,j-1)-\sum_{q=1}^{j-1}f_k(p_q)\right)\ \ \ \ \ (i\ge p_j^2) \]

    由于这是完全积性函数,那么有 \(f(r\cdot p_j)=f(r)\cdot f(p_j)\),后面的一大坨相当于计算合法的 \(r\). 可以发现完全积性函数的性质省去了枚举质数幂次的操作。

    注意到后面 \(i\ge p_j^2\) 的条件。实际上,若 \(i<p_j^2\),那么就不可能找到合法的 \(r\),所以直接 \(g(i,j)\leftarrow g(i,j-1)\) 即可。

    需要注意的是,我们需要初始化 \(g(i,0)\),所以拆分出来的 \(f\) 的前缀和应当比较好求。

    容易发现,\(P_k(n)=g(n,cnt)\)\(cnt\) 是筛出质数的个数)。


    关于复杂度:\(\text{It is said that}\) 计算 \(h(i,j)\) 的复杂度为 \(\mathcal O(n^{1-\epsilon})\),计算 \(g(i,j)\) 的复杂度为 \(\mathcal O\left(\frac{n^{\frac{3}{4}}}{\log n}\right)\)(需要注意,拆成 \(k\) 个函数就要乘上常数 \(k\))。证明不会,可能会一直咕下去。

    代码实现

    \(\text{Update on 2022.1.29}\)

    救命之前忘了写一个很重要的东西,在 \(h(i,j)\) 函数中的第一层循环要加上 p[i]*p[i]<=x 的特判!

    然后还要注意我们只证明个数是 \(\sqrt n\) 级别,所以数组开大几倍!


    参照某谷 模板题

    在实现过程中,函数的 \(i,j\) 参数的 个数 都是小于 \(\sqrt n\) 的。\(j\) 的原因就不阐述,对于 \(i\),可以观察到它是一个类似 \(\left \lfloor \frac{\left \lfloor \frac{n}{a}\right \rfloor}{b}\right \rfloor\) 的向下嵌套的结构,可以证明它的取值属于 \(\left \lfloor \frac{n}{a}\right \rfloor\) 的取值集合,这样不仅可以证明它的个数在 \(\sqrt n\) 之内,而且还可以进行整除分块预处理。以下给出感性证明:不妨将 \(\left \lfloor \frac{n}{a}\right \rfloor\) 中的分母替换成 \(a'\) 使得 \(a'\mid n\)\(\frac{n}{a'}=\left \lfloor \frac{n}{a}\right \rfloor\),这样嵌套就可以换成 \(\left \lfloor \frac{n}{a'b}\right \rfloor\),结论得证。

    于是我们将合理的 \(i\) 存起来,需要注意的是,\(i\) 的取值可能很大,所以对于 \(i\le \sqrt n\) 的部分用 \(i\) 来查找,另一部分用 \(n/i\) 来查找(这个玩意实际上就是整除分块时有 \(i=n/l\),对应的 \(r\) 就是 \(n/i\))。

    回到此题,可知 \(F(n)=n^2-n\),可以分成 \(f_1,f_2\) 来计算,然后合并即可。

    #include <cstdio>
    #define print(x,y) write(x),putchar(y)
    
    template <class T>
    inline T read(const T sample) {
    	T x=0; char s; bool f=0;
    	while((s=getchar())>'9' || s<'0')
    		f |= (s=='-');
    	while(s>='0' && s<='9')
    		x = (x<<1)+(x<<3)+(s^48),
    		s = getchar();
    	return f?-x:x;
    }
    
    template <class T>
    inline void write(T x) {
    	static int writ[50],tp=0;
    	if(x<0) putchar('-'),x=-x;
    	do writ[++tp] = x-x/10*10, x/=10; while(x);
    	while(tp) putchar(writ[tp--]^48);
    }
    
    #include <cmath>
    #include <iostream>
    using namespace std;
    typedef long long ll;
    
    const int maxn = 1e6+5;
    const int mod = 1e9+7;
    const int inv = 166666668;
    
    bool is[maxn];
    ll n,id[maxn];
    int S,f[2][maxn],pc,p[maxn];
    int g[2][maxn],tot,ind[2][maxn];
    
    inline int inc(int x,int y) {
    	return x+y>=mod?x+y-mod:x+y;
    }
    
    inline int dec(int x,int y) {
    	return x-y<0?x-y+mod:x-y;
    }
    
    void sieve(int n) {
    	for(int i=2;i<=n;++i) {
    		if(!is[i]) {
    			p[++pc] = i;
    			f[0][pc] = inc(f[0][pc-1],i);
    			f[1][pc] = inc(f[1][pc-1],1ll*i*i%mod);
    		}
    		for(int j=1; j<=pc && i*p[j]<=n; ++j) {
    			is[i*p[j]] = 1;
    			if(i%p[j]==0) break;
    		}
    	}
    }
    
    inline int getSum(ll n,bool d) {
    	int ret; n %= mod; // qwq
    	if(!d) ret = (n*(n+1)>>1)%mod;
    	else ret = n*(n+1)%mod*(n*2+1)%mod*inv%mod;
    	return (ret-1<0?ret-1+mod:ret-1);
    }
    
    void getf() {
    	ll l,r,x;
    	for(l=1;l<=n;l=r+1) {
    		r = min(n,n/(n/l)); 
    		id[++tot] = (x=n/l);
    		if(x<=S) ind[0][x]=tot;
    		else ind[1][n/x]=tot; // for better restoration
    		g[0][tot] = getSum(x,0);
    		g[1][tot] = getSum(x,1);
    	}
        // 用滚动数组递推 g(i,j)
    	for(int i=1;i<=pc;++i) {
    		for(int j=1; j<=tot && 1ll*p[i]*p[i]<=id[j]; ++j) {
    			int pre = (id[j]/p[i]<=S)?
    					ind[0][id[j]/p[i]]:ind[1][n/(id[j]/p[i])];
    			g[0][j] = dec(g[0][j],1ll*p[i]*dec(g[0][pre],f[0][i-1])%mod);
    			g[1][j] = dec(g[1][j],1ll*p[i]*p[i]%mod*dec(g[1][pre],f[1][i-1])%mod);
    		}
    	}
    }
    
    int h(ll x,int y) {
    	if(x<=p[y]) return 0;
    	int ID = (x<=S)?ind[0][x]:ind[1][n/x];
    	int ans = dec(dec(g[1][ID],g[0][ID]),dec(f[1][y],f[0][y]));
    	for(int i=y+1; i<=pc && 1ll*p[i]*p[i]<=x; ++i) {
    		for(ll tmp=p[i];;tmp=tmp*p[i]) {
    			if(tmp*p[i]>x) break;
    			ans = inc(ans,inc(
    				tmp*(tmp-1)%mod*h(x/tmp,i)%mod,
    				tmp*p[i]%mod*((tmp*p[i]-1)%mod)%mod
    			));
    		}
    	}
    	return ans;
    }
    
    int main() {
    	n=read(9ll);
    	sieve(S=sqrt(n)); getf(); 
    	print(inc(h(n,0),1),'\n');
    	return 0;
    }
    

    题目小结

    例 1. \(\text{LOJ - 6053 }\)简单的函数

    这题的 \(F(n)\) 只在 \(p^k\) 上有定义,而且看不出来是个啥函数,难道不行?

    事实上,我们只需要保证 质数 的部分算对即可,因为其它数最后会被筛掉。也就是说,我们其实只需要研究 \(F(p)=p\text{ xor }1\) 这个函数。很明显,当 \(p\) 为奇质数时,\(F(p)=p-1\),反之(也就是 \(p=2\) 时)有 \(F(p)=p+1\).

    由于只有 \(2\) 一个特例,我们不妨计算出 \(F(n)=n-1\) 这个函数的答案再进行修改。事实上,我们只算错了 \(P(n)\),且第二维等于零(会计算 \(2\) )的 \(h(i,j)\) 只有一个即 \(h(n,0)\),所以只需要将答案加上 \(2\) 即可。别忘了最后算上 \(1\).

  • 相关阅读:
    vmware Unable to open kernel device "\.Globalvmx86": The system cannot find the file 的解决方法
    nc和telnet配合使用
    linux下批量替换文件内容
    Linux动态库的导出控制
    goang Receiver & interface
    Go与C语言的互操作 cgo
    Go fsm
    Git多账号登陆
    mysql 安装与配置、使用
    Reverse Integer
  • 原文地址:https://www.cnblogs.com/AWhiteWall/p/15841383.html
Copyright © 2020-2023  润新知