<学习笔记> 杜教筛
杜教筛
处理数论函数的前缀和问题,可以在低于线性的复杂度里求出
对于任意一个数论函数
那么可以得到:
时间复杂度最小为
一般情况下,可以设出方便求前缀的
例如根据
单位元
根据狄利克雷卷积得到几个性质:
对于一,我们直接展开就是
对于二,还是展开
证明
设
我们还知道
对于三,我们可以另
一些题#
DZY Loves Math IV #
发现
我们设
我们另
尝试枚举
当
code
#include<bits/stdc++.h>
const int mod=998244353;
#define int long long
using namespace std;
const int N=1e5;
int prime[N+5];
bool nprime[N+5];
int np[N+5],phi[N+5],sumphi[N+5];
void get_prime(){
phi[1]=1;
for(int i=2;i<=N;i++){
if(!nprime[i]){
prime[++prime[0]]=i;
phi[i]=i-1;
np[i]=i;
}
for(int j=1;j<=prime[0] && i*prime[j]<=N;j++){
int k=i*prime[j];
nprime[k]=1;
if(i%prime[j]==0){
np[k]=np[i];
phi[k]=phi[i]*prime[j];
break;
}
phi[k]=phi[i]*phi[prime[j]];
np[k]=np[i]*np[prime[j]];
}
}
for(int i=1;i<=N;i++) sumphi[i]=(sumphi[i-1]+phi[i])%mod;
}
unordered_map<int,int> Sphi,Ans[N+5];
int get_phi(int n){
if(n<=N) return sumphi[n];
if(Sphi[n]) return Sphi[n];
int l=2,r;
int ans=(n+1)*n/2;
while(l<=n){
int t=n/l;
int r=n/t;
ans=(ans-get_phi(n/l)*(r-l+1)%mod+mod)%mod;
l=r+1;
}
return Sphi[n]=ans;
}
int S(int n,int m){
if(m==0) return 0;
if(Ans[n][m]) return Ans[n][m];
if(n==1){
Ans[n][m]=get_phi(m);
return Ans[n][m];
}
int ans=n/np[n];
int sum=0;
int t=sqrt(np[n]);
for(int d=1;d<=t;d++){
if(np[n]%d==0){
sum=(sum+phi[np[n]/d]*S(d,m/d)%mod)%mod;
if(d!=np[n]/d){
int tmp=np[n]/d;
sum=(sum+phi[np[n]/tmp]*S(tmp,m/tmp)%mod)%mod;
}
}
}
return Ans[n][m]=ans*sum%mod;
}
signed main(){
get_prime();
int T;
scanf("%lld",&T);
while(T--){
int n,m;
scanf("%lld%lld",&n,&m);
int ans=0;
for(int i=1;i<=n;i++){
ans=(ans+S(i,m))%mod;
}
printf("%lld\n",ans);
}
}
Pyh 的求/毒瘤之神的考验 #
比上面那道题多一个多测,并且
发现上面做法行不通,答案是没啥关系,唯一的关系是我都不会。
有一个结论:
证明
如果在乘上一个
那么我们那它推式子
变为枚举
另
我们可以套路的另
另
那么答案就变成了
但是后面的那个东西好像无法直接整除分块,所以我们考虑根号分治,直接预处理一部分答案,我们另
我们定义阈值为
复杂度:预处理的话是
code
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
const int N=1e5;
const int B=50;
int prime[N+2];
bool nprime[N+2];
int phi[N+2],mul[N+2],F[N+2],inv[N+2],sumF[N+2];
vector<int> G[N+2];
vector<int> Ans[B+2][B+2];
inline int qpow(int x,int p){
int ans=1;
while(p){
if(p&1) ans=(1ll*ans*x)%mod;
x=(1ll*x*x)%mod;
p>>=1;
}
return ans;
}
inline void get_prime(){
phi[1]=mul[1]=1;
for(int i=2;i<=N;i++){
if(!nprime[i]){
phi[i]=i-1;
mul[i]=mod-1;
prime[++prime[0]]=i;
}
for(int j=1;j<=prime[0] && i*prime[j]<=N;j++){
int k=i*prime[j];
nprime[k]=1;
if(i%prime[j]==0){
mul[k]=0;
phi[k]=1ll*phi[i]*prime[j]%mod;
break;
}
phi[k]=1ll*phi[i]*phi[prime[j]]%mod;
mul[k]=1ll*mul[i]*mul[prime[j]]%mod;
}
}
}
int ask_G(int x,int y){
if(y<=B) return G[x][y];
int ans=0;
for(int i=1;i<=x;i++) ans=(1ll*ans+phi[i*y])%mod;
return ans;
}
inline int read(){
int x(0);bool f(0);char ch=getchar();
for(;ch<'0'||ch>'9';ch=getchar()) f^=ch=='-';
for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+(ch^48);
return f?x=-x:x;
}
inline void write(int x){
x<0?x=-x,putchar('-'):0;static short Sta[50],top(0);
do{Sta[++top]=x%10;x/=10;}while(x);
while(top) putchar(Sta[top--]|48);
putchar('\n');
}
signed main(){
get_prime();
for(int i=1;i<=N;i++) inv[i]=qpow(i,mod-2);
for(int d=1;d<=N;d++){
for(int k=1;d*k<=N;k++){
int T=d*k;
F[T]=(1ll*F[T]+1ll*d*inv[phi[d]]%mod*mul[k]%mod)%mod;
}
}
for(int i=1;i<=N;i++) sumF[i]=(1ll*sumF[i-1]+F[i])%mod;
G[0].resize(N+5);
for(int x=1;x<=N;x++){
G[x].resize(N/x+5);
for(int y=1;x*y<=N;y++){
G[x][y]=(1ll*G[x-1][y]+phi[x*y])%mod;
}
}
for(int y=1;y<=B;++y){
for(int z=1;z<=B;++z){
Ans[y][z].resize(min(N/y,N/z)+5);
for(int x=1;x*y<=N && x*z<=N;++x){
Ans[y][z][x]=(1ll*Ans[y][z][x-1]+1ll*F[x]*ask_G(y,x)%mod*ask_G(z,x)%mod)%mod;
}
}
}
int T;
T=read();
while(T--){
int n,m;
n=read(),m=read();
if(n>m) swap(n,m);
long long ans=0;
for(int t=1;t<=m/B;t++){
ans=(ans+1ll*F[t]*G[n/t][t]%mod*G[m/t][t]%mod)%mod;
}
int l=m/B+1,r=0;
while(l<=n){
int ta=n/l,tb=m/l;
if(!ta || !tb) break;
r=min(n/ta,m/tb);
ans=(ans+1ll*(Ans[ta][tb][r]-Ans[ta][tb][l-1]+mod)%mod)%mod;
l=r+1;
}
write(ans);
}
}
简单的最大公约数#
首先莫反做法,设
那么答案就是
后面的式子相当于
直接用杜教筛就可以。
然后假如用欧拉反演的话,让
那么就有
code
//简单的最大公约数
#include<bits/stdc++.h>
#define int unsigned long long
using namespace std;
const int N=3*1e6;
int prime[N+5];
bool nprime[N+5];
int sumphi[N+5],phi[N+5];
unordered_map<int,int> Sumphi;
inline void get_prime(){
phi[1]=1;
for(register int i=2;i<=N;++i){
if(!nprime[i]){
prime[++prime[0]]=i;
phi[i]=i-1;
}
for(register int j=1;j<=prime[0] && i*prime[j]<=N;j++){
int k=prime[j]*i;
nprime[k]=1;
if(i%prime[j]==0){
phi[k]=phi[i]*prime[j];
break;
}
phi[k]=phi[i]*phi[prime[j]];
}
}
for(register int i=1;i<=N;++i){
sumphi[i]=sumphi[i-1]+phi[i];
}
}
inline int get_phi(int n){
if(n<=N) return sumphi[n];
if(Sumphi[n]) return Sumphi[n];
int ans=0;
if(n&1)
ans=(n+1)/2*n;
else
ans=n/2*(n+1);
int l=2;
while(l<=n){
int t=n/l;
int r=n/t;
ans-=(r-l+1)*get_phi(n/l);
l=r+1;
}
return Sumphi[n]=ans;
}
int qpow(int x,int p){
int ans=1;
while(p){
if(p&1) ans=(ans*x);
x=(x*x);
p>>=1;
}
return ans;
}
signed main(){
get_prime();
int n,m;
cin>>n>>m;
int ans=0;
int l=1;
while(l<=m){
int t=m/l;
int r=min(m/t,m);
ans=(ans+qpow(t,n)*(get_phi(r)-get_phi(l-1)));
l=r+1;
}
cout<<ans<<endl;
}
作者:bloss
出处:https://www.cnblogs.com/jinjiaqioi/p/17978192
版权:本作品采用「署名-非商业性使用-相同方式共享 4.0 国际」许可协议进行许可。
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 单线程的Redis速度为什么快?
· 展开说说关于C#中ORM框架的用法!
· Pantheons:用 TypeScript 打造主流大模型对话的一站式集成库
· SQL Server 2025 AI相关能力初探
· 为什么 退出登录 或 修改密码 无法使 token 失效