ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

Powerful Number 筛学习笔记

2021-03-31 20:02:21  阅读:224  来源: 互联网

标签:val int Powerful Number 笔记 number rg nans


Powerful Number 筛学习笔记

用途

\(Powerful\ number\) 筛可以用来求出一类积性函数的前缀和,最快可以达到根号复杂度。

实现

\(Powerful\ number\) 的定义是每个质因子次数都 \(\ge 2\) 的数。

有如下的性质:

\(1\)、一个 \(Powerful\ number\) 一定可以表示为 \(a^2b^3\) 的形式。

\(2\)、\(n\) 以内的 \(Powerful\ number\) 个数是 \(O(\sqrt n)\) 级别的。

所以找 \(Powerful\ number\) 可以直接暴力 \(dfs\)。

如果要求的函数是 \(f\),那么我们需要找到一个积性函数 \(g\),使得 \(f\) 和 \(g\) 在质数处的取值相同。

同时还要找到一个积性函数 \(h\),使得 \(f=g*h\)。

根据狄利克雷卷积的定义

\(f(p)=g(p)h(1)+g(1)h(p)=g(p)+h(p)=f(p)+h(p)\)。

所以 \(h(p)=0\),因为 \(h\) 是一个积性函数,所以所有非 \(Powerful\ number\) 在 \(h\) 函数中的取值都是 \(0\)。

\(\sum_{i=1}^nf(i)=\sum_{i=1}^n\sum_{d|i}h(d)g(\frac{i}{d})=\sum_{d=1}^nh(d)\sum_{j=1}^{\frac{n}{d}}g(j)\)。

因为 \(h\) 只在 \(Powerful\ number\) 处有值,所以我们只需要求出 \(\sqrt{n}\) 个 \(g\) 函数的前缀和即可。

例题

题目描述

给定一个积性函数 \(f\),满足 \(f(1)=1\),并且对于任意质数 \(p\) 和正整数 \(e\),都有 \(f(p^e)=p^k\),\(k\) 为给定的数,\(n \leq 10^{13},k \leq 20\)。

分析

构造积性函数 \(g\),满足对于任意 \(x\),都有 \(g(x)=x^k\)。

\(f(p^2)=g(p^2)h(1)+g(p)h(p)+g(1)h(p^2)\)。

那么 \(p^{k}=p^{2k}+h(p^2)\),\(h(p^2)=p^{k}-p^{2k}\)。

\(f(p^3)=g(p^3)h(1)+g(p^2)h(p)+g(p)h(p^2)+g(1)h(p^3)\)。

\(p^{k}=p^{3k}+p^k(p^{k}-p^{2k})+h(p^3)\),\(h(p^3)=p^{k}-p^{2k}\)。

多算几项就会发现 \(h(p^e)=p^{k}-p^{2k},e \ge 2\),

直接做就行了。

代码

#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cmath>
#include<map>
#define rg register
template<typename T>void read(rg T& x){
	x=0;rg int fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	x*=fh;
}
const int maxn=1e7+5,maxm=35,mod=1e9+7;
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
inline int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
int pri[maxn],k,sqr,b[maxm],c[maxm][maxm],ny[maxm],tot,val[maxn],mi[maxn],ans;
bool not_pri[maxn];
long long n,w[maxn];
void pre(){
	for(rg int i=0;i<maxm;i++) c[i][0]=1;
	for(rg int i=1;i<maxm;i++){
		for(rg int j=1;j<=i;j++){
			c[i][j]=addmod(c[i-1][j-1],c[i-1][j]);
		}
	}
	ny[1]=1;
	for(rg int i=2;i<maxm;i++) ny[i]=mulmod(mod-mod/i,ny[mod%i]);
	b[0]=1;
	for(rg int i=1;i<=20;i++){
		for(rg int j=0;j<=i-1;j++){
			b[i]=addmod(b[i],mulmod(c[i+1][j],b[j]));
		}
		b[i]=delmod(0,mulmod(b[i],ny[i+1]));
	}
}
std::map<int,int> mp;
int getsum(rg int val){
	if(mp.find(val)!=mp.end()) return mp[val];
	val++;
	rg int nans=0,tmp=val;
	for(rg int i=k;i>=0;i--){
		nans=addmod(nans,mulmod(c[k+1][i],mulmod(b[i],tmp)));
		tmp=mulmod(tmp,val);
	}
	nans=mulmod(nans,ny[k+1]);
	return mp[val-1]=nans;
}
void xxs(rg int mmax){
	not_pri[0]=not_pri[1]=1;
	for(rg int i=2;i<=mmax;i++){
		if(!not_pri[i]){
			pri[++pri[0]]=i;
			mi[pri[0]]=delmod(ksm(i,k),ksm(i,k<<1));
		}
		for(rg int j=1;j<=pri[0] && i*pri[j]<=mmax;j++){
			not_pri[i*pri[j]]=1;
			if(i%pri[j]==0) break;
		}
	}
}
void dfs(rg int now,rg long long nw,rg int nval){
	w[++tot]=nw,val[tot]=nval;
	for(rg int i=now;i<=pri[0] && nw<=n/pri[i]/pri[i];i++){
		rg long long tmp=1LL*nw*pri[i];
		for(;tmp<=n/pri[i];tmp*=pri[i]) dfs(i+1,tmp*pri[i],mulmod(nval,mi[i]));
	}
}
int main(){
	pre();
	read(n),read(k);
	sqr=sqrt(n);
	xxs(sqr+5);
	dfs(1,1,1);
	for(rg int i=1;i<=tot;i++){
		ans=addmod(ans,mulmod(val[i],getsum((n/w[i])%mod)));
	}
	printf("%d\n",ans);
	return 0;
}

标签:val,int,Powerful,Number,笔记,number,rg,nans
来源: https://www.cnblogs.com/liuchanglc/p/14603461.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有