zoukankan      html  css  js  c++  java
  • Min25 筛与 Powerful Numbers

    Min25 筛与 Powerful Numbers

    Min25 筛

    大喊一声 Min25 NB!!!

    这是一个非常神奇的东西,用于求更加普遍的积性函数的前缀和。

    比如我们要求 (sum_{i=1}^{n}f(i)),其中 (f(1)=1)。我们考虑将质数与合数分开考虑。(由于 (1) 这俩都不是,我们先不考虑)所以答案就等于质数的答案加上合数的答案再加上 (f(1))

    但是这样还是不够优美。因为合数的范围太广泛了,我们考虑对于每个质数再分组。

    在这里,我们设 (sigma(x))(x) 的最小质因子。

    定义 (S(n,k)) 表示第 (2) 到第 (n) 个数中,(sigma (x) ge pri_k) 的数的 (f(x)) 之和。其中 (p_k) 表示第 (k) 大的质数。显然有 (max{p_k}le sqrt n)。这也决定了用到的质数不会很多。

    我们同样将质数和合数分开考虑。由于 (f) 是积性函数,所以有 (f(p_1^{k_1})f(p_2^{k^2})=f(p_1^{k_1}p_2^{k^2}))。这启发我们可以枚举质数幂进行转移。

    所以有

    [S(n,k)=sum_{ige k,p_ile n}f(p_i)+sum_{jge k}sum_{ege 1,p_j^{e+1}le n}(f(p_j^e)cdot S(lfloorfrac{n}{p_j^e} floor,j+1)+f(p_j^{e+1})) ]

    稍微解释一下就是前半部分是满足条件的质数的和。 (S(lfloorfrac{n}{p_j^e} floor,j+1)) 中的数的最小质因子 (ge p_{j+1}>p_j),所以这里我们枚举一个没有出现过的质数幂乘上去。同时由于 (f) 是积性函数,我们要对这一部分值全部乘上 (f(p_j^e))。同时注意答案中没有包含 (f(p_j^{e+1})) 这一项(因为我们是从 (2) 开始计数的,在前面做乘法的时候由于 (1) 的缺失导致 (p_j^e) 的缺失),要补上。

    最后的答案就是 (S(n,1)+f(1))

    这个东西可以递归进行计算,但是显然需要一些条件:

    • (f(p)) 是一个关于 (p) 的项数较小的多项式或可以快速求值。
    • (f(p^c)) 能够较方便的计算,即质数幂处的 (f) 值能够较方便的计算。
    • (sum_{ige k,p_ile n}f(p_i)) 能够快速计算。然后这东西可以拆成前缀和然后作差。(ile k) 的部分由于 (k) 很小,且 (f(p)) 是低次多项式,可以在线性筛的时候顺便处理掉,问题转化为快速求出当 (n) 较大时质数处的 (f) 的前缀和。

    第一个和第二个条件我们似乎没有什么办法克服,全看出题人是否良心。

    我们来考虑能否快速求出质数处的 (f) 的前缀和,设 (sum_k=sum_{ple k,pin mathrm{prime}}f(p))

    首先对上面的递推式进行观察,会发现我们需要的 (sum_k) 的个数其实很少,都是形如 (sum_{lfloorfrac n x floor}) 这样的值。

    显然可以证明这样的下标最多只有 (O(sqrt n )) 个,实在不放心可以打个表试试。

    然后我们发现,此时我们对于合数的值并不关心,即我们只需要构造一个 (g) 使得 (g(p)=f(p)),求出 (g) 的前缀和后再将合数位置处的 (g) 的值剪掉,同样能够得到我们想要得到答案。

    我们可以通过构造使得 (g) 是一个完全积性函数,即对于任意 (i,j) 均有 (g(i)g(j)=g(ij))

    (G(n,k)) 表示 (forall xin [2,n],xin mathrm{prime} or sigma(x)> p_k) 这样的 (x)(g(x)) 之和。其实这个状态非常类似与我们在进行埃氏筛的时候,枚举到第 (k) 个质数且该轮已筛完的状态。

    于是我们枚举当前应该筛掉哪些数,有

    [G(n,k)=G(n,k-1)-g(p_k)(G(lfloor frac{n}{p_k} floor,k-1)-sum_{p_{k-1}}) ]

    我们当前要筛掉 (sigma (x)=p_k) 的合数。

    (p_k^2ge n) 时,显然我们无数可筛。此时 (G(n,k)=G(n,k-1))

    (p_k^2le n) 时,由于我们构造了 (g) 使其为完全积性函数,所以减掉 (g(p_k)G(lfloor frac{n}{p_k} floor,k-1)),就相当于减掉了这一部分合数的 (g(x)) 之和。但是根据定义,(G(lfloor frac{n}{p_k} floor,k-1)) 包含了一部分 (le p_k) 的质数,这样会使得 (sigma (x) eq p_k),所以要加上 (g(p_k)sum_{p_{k-1}})

    显然答案为 (G(n,mathrm{cnt})),其中 (mathrm{cnt})(le sqrt n) 的质数个数。

    然后根据前文所提到的,我们最后需要的只有 (O(sqrt n))(G(N,mathrm{cnt})),直接做就完事了。

    初始状态也非常显然:(G(n,0)=sum_{i=1}^{n}g(i))

    所以我们构造的这个函数应该需要比较好算。

    事实上若 (f(x)) 是一个 (k) 次多项式,我们有(f(x)=sum_{i=0}^ka_ix^i)

    我们把 (f(x)) 拆开,令 (g_k(x)=x^k),则有 (f(x)=sum_{i=0}^ka_ig_i(x))

    显然每一个 (g_k) 都是一个完全积性函数,而且它的前缀和也非常好算(拉格朗日插值即可),所以我们对每一个 (g_k) 筛一遍再合并即可。

    Q:如果 k 很大咋办呢?

    A:k 很大... k 很大...... 自求多福。——zzq

    然后我们就做完了。

    我也不太会证明复杂度,反正就能跑 (10^{10})(10^{11}) 可以抢救一下。

    实际上这个算法比杜教筛快,限制还比杜教筛小。

    杜教筛,永远的垃圾。——zzq

    「LOJ 6053」 简单的函数

    模板题。对于所有奇质数 (p),有 (f(p)=p-1)。我们先令 (f(p)=p) 然后再减掉素数个数即可。

    同时注意考虑唯一的偶质数 (2) 的贡献。

    /*---Author:HenryHuang---*/
    /*---Never Settle---*/
    #include<bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    const ll maxn=5e5+5;
    const ll p=1e9+7;
    const ll inv2=5e8+4;
    ll n,P;
    ll pri[maxn],pp[maxn],cnt;
    ll id[maxn],id2[maxn],dex;
    ll w[maxn],g[maxn],h[maxn],f[maxn];
    ll pre[maxn];
    ll F(ll x,ll y){
    	return x^y;
    } 
    void init(){
    	for(ll i=2;i<=P;++i){
    		if(!pp[i]) pri[++cnt]=i,pre[cnt]=(pre[cnt-1]+i)%p;
    		for(ll j=1;j<=cnt&&pri[j]*i<=P;++j){
    			pp[pri[j]*i]=1;
    			if(i%pri[j]==0) break;
    		}
    	}
    }
    
    ll getid(ll x){
    	if(x<=P) return id[x];
    	else return id2[n/x];
    }
    void prework(){
    	for(ll l=1,r,v;l<=n;l=r+1){
    		v=(n/l),r=(n/v);
    		if(v<=P) id[v]=++dex;
    		else id2[r]=++dex;
    		w[dex]=v;v%=p;
    		g[dex]=1ll*(2+v)*(v-1+p)%p*inv2%p;
    		h[dex]=(v-1);
    	}
    	for(ll i=1;i<=cnt;++i){
    		ll z=pri[i];
    		for(ll j=1;j<=dex&&z*z<=w[j];++j){
    			ll op=getid(w[j]/z);
    			g[j]=(g[j]-z*(g[op]-pre[i-1]+p)%p+p)%p;
    			h[j]=(h[j]-h[op]+i-1+p)%p;
    		}
    	}
    	for(ll j=1;j<=dex;++j) f[j]=(g[j]-h[j]+p)%p;
    }
    ll solve(ll x,ll i){
    	if(x<=1||pri[i]>x) return 0;
    	ll ans=f[getid(x)]-(pre[i-1]-(i-1));
    	if(i==1) ans+=2;
    	for(ll j=i;j<=cnt&&pri[j]*pri[j]<=x;++j){
    		ll r=pri[j];
    		for(ll e=1;r*pri[j]<=x;++e,r*=pri[j]){
    			ans=(ans+F(pri[j],e)*solve(x/r,j+1)+F(pri[j],e+1))%p;
    		}
    	}
    	return ans;
    }
    int main(){
    	ios::sync_with_stdio(0);
    	cin.tie(0),cout.tie(0);
    	cin>>n;P=sqrt(n+0.5);
    	init();
    	prework();
    	cout<<(solve(n,1)+1)%p<<'
    ';
    	return 0;
    }
    


    Powerful Numbers

    这也是一个非常神奇的东西,用于求极为特殊的积性函数的前缀和。

    比如我们要求 (sum_{i=1}^{n}f(i)),但是这个 (f) 非常神奇,其满足 (f(1)=1,f(p)=1),其中 (p) 为质数。

    然后我们发现这和我们一个非常喜欢的积性函数 (g(x)=1) 非常像,大家都会求它的前缀和。它们在质数处的值是一样的,即 (f(p)=g(p))

    我们尝试构造一个积性函数 (h(x)) 使得 (sum_{ij=x}h(i)g(j)=f(x))

    首先由 (f(1)=1) 可推出 (h(1)=1)

    再由 (f(p)=1) 可以推出 (f(p)=h(1)g(p)+h(p)g(1)=1+h(p))

    从这里我们发现 (h(p)=0),也就是说,(h) 在质数位置上的值均为 (0)

    再进一步,由于 (h) 是一个积性函数,由 (h(p_ip_j)=h(p_i)h(p_j))(h) 在所有存在一个质数因子的次数 (=1) 的位置均为 (0)

    换句话说,只有当一个数,其包含的所有质因子的次数均 (ge 2)(h) 才会有值。而这样的数,就被称为 ( exttt{Powerful Numbers})

    可以证明其个数为 (O(sqrt n)) 的。你实在不放心可以打个表跑一跑。

    如何求?枚举质因子爆搜就完事了,但需要注意必要的剪枝。

    由于 ( exttt{Powerful Numbers}) 有这样美好的性质,我们可以把原来的式子转换一下:

    [egin{aligned} &sum_{x=1}^nf(x)\ =&sum_{x=1}^nsum_{ij=x}h(i)g(j)\ =&sum_{i=1}^nh(i)sum_{j=1}^{lfloor frac n i floor}g(j) end{aligned} ]

    本质上我们需要求出所有 (ijle n)(h(i)g(j)),所以换个方式枚举即可。

    现在我们的问题就变成了:求 ( exttt{Powerful Numbers})(h(x)),以及求 (g) 的前缀和。

    在这个例子中,(g) 的前缀和大家都会求。即 (sum_{i=1}^ng(i)=n)

    问题变为如何快速求出 (h(x))

    事实上你经常可以通过质数幂处的 (f(p^c)) 来猜出 (h(p^c)) 的值,又因为其是一个积性函数,你就能在爆搜的同时把答案求出来了。

    事实上如果你猜不出来可能需要把质数幂处的值全部整出来然后暴力多项式除法?

    不管能猜就猜猜不出来再说

    总结一下,能够用这种筛法的条件

    • 你能构造出一个积性函数 (g) 使得 (f(p)=g(p))
    • 这个 (g) 的前缀和比较好算。

    Example

    (f(x)=p_1^{lfloor frac{c_1}{2} floor}p_2^{lfloor frac{c_2}{2} floor}cdots p_k^{lfloor frac{c_k}{2} floor}),其中 (p) 为对 (x) 唯一分解后的质因子,(c) 为每个质因子的次数。求 (sum_{i=1}^nf(i))


    我们发现这东西显然是满足 (f(p)=1) 的,所以非常愉快地令 (g(x)=1)

    现在问题转变为求 (h(x))

    我们多列几项:

    [egin{aligned} f(1)=h(1)g(1)& ightarrow h(1)=1\ f(p)=h(1)g(p)+h(p)g(1)& ightarrow h(p)=0\ f(p^2)=h(1)g(p^2)+h(p)g(p)+h(p^2)g(1)& ightarrow h(p^2)=p-1\ f(p^3)=h(1)g(p^3)+h(p)g(p^2)+h(p^2)g(p)+h(p^3)g(1)& ightarrow h(p^3)=0\ f(p^4)=h(1)g(p^4)+h(p)g(p^3)+h(p^2)g(p^2)+h(p^3)g(p)+h(p^4)g(1)& ightarrow h(p^4)=p^2-p=p(p-1)\ end{aligned} ]

    我觉得列到这里已经差不多能猜出来了。

    (k) 为奇数时,(h(p^k)=0)

    (k)(2) 时,(h(p^k)=p-1)

    (k) 为其他偶数时,(h(p^k)=(p-1)p^{frac{k-2}{2}})

    然后算就完事了。

    时间复杂度为 (O(sqrt n ))。事实上由于当 (k) 是奇数时 (h(p^k)) 均为 (0),所以可以进一步剪枝。

    /*---Author:HenryHuang---*/
    /*---Never Settle---*/
    #include<bits/stdc++.h>
    using namespace std;
    const int maxn=3.5e6+5;
    typedef unsigned long long ll;
    ll pri[maxn],pp[maxn],p,cnt;
    void init(){
    	for(ll i=2;i<=p;++i){
    		if(!pp[i]) pri[++cnt]=i;
    		for(ll j=1;j<=cnt&&pri[j]*i<=p;++j){
    			pp[pri[j]*i]=1;
    			if(i%pri[j]==0) break;
    		}
    	}
    }
    ll ans,n,owo;
    void dfs(ll num,ll now,ll phi){
    	if(now==cnt+1||num*pri[now]>n||num*pri[now]*pri[now]>n){
    		ans+=phi*(n/num);
    		++owo;
    		return ;
    	}
    	dfs(num,now+1,phi);
    	for(int i=2;;i+=2){
    		num*=pri[now];
    		if(num>n) break;
    		num*=pri[now];
    		if(num>n) break;
    		if(i==2) phi*=(pri[now]-1);
    		else phi*=pri[now];
    		dfs(num,now+1,phi);
    	}
    }
    int main(){
    	ios::sync_with_stdio(0);
    	cin.tie(0),cout.tie(0);
    	int T;cin>>T;
    	p=sqrt(1e13);
    	init();
    	while(T--){
    		cin>>n;ans=0;
    		dfs(1,1,1);
    		cout<<ans<<'
    ';
    	}
    	return 0;
    }
    

    Submit Link

    在繁华中沉淀自我,在乱世中静静伫立,一笔一划,雕刻时光。
  • 相关阅读:
    Python生成验证码
    Django设置
    OpenStack安装后检查流程总结
    利用src.rpm包修改源码后重新制作rpm包
    Python知识点:distutils常用子模块
    libvirt, libvirt-python, libvirtd 关系浅析
    Python知识点: os.popen
    Python知识点: __import__
    修改initrd.img里ko文件的一个小tips
    关于openstack自动化安装的一点思考
  • 原文地址:https://www.cnblogs.com/HenryHuang-Never-Settle/p/14422983.html
Copyright © 2011-2022 走看看