当前位置:网站首页>min25筛
min25筛
2022-06-12 02:32:00 【doctorZ_】
定义积性函数 f ( x ) f(x) f(x), f ( p c ) f(p^c) f(pc)可以表示为多个完全积性函数 h i ( p c ) h_i(p^c) hi(pc)相加,求
∑ i = 1 n f ( i ) \sum_{i=1}^nf(i) i=1∑nf(i)
∑ i = 1 n f ( i ) = ∑ p ≤ n f ( p ) + ∑ i i s n o t a p r i m e f ( i ) \sum_{i=1}^nf(i)=\sum_{p\le n}f(p)+\sum_{i\ is\ not\ a\ prime}f(i) i=1∑nf(i)=p≤n∑f(p)+i is not a prime∑f(i)
设 g ( i , j , k ) g(i,j,k) g(i,j,k)表示 ∑ w = 1 n [ w ∈ p r i m e ∣ m i n p > p j ] h k ( w ) \sum_{w=1}^n[w\in prime|minp>p_j]h_k(w) ∑w=1n[w∈prime∣minp>pj]hk(w)
g ( i , j , k ) = g ( i , j − 1 , k ) − [ g ( ⌊ i p j ⌋ , j − 1 , k ) − g ( p j − 1 , j − 1 , k ) ] h k ( p j ) g(i,j,k)=g(i,j-1,k)-[g(\lfloor\frac{i}{p_j}\rfloor,j-1,k)-g(p_j-1,j-1,k)]h_k(p_j) g(i,j,k)=g(i,j−1,k)−[g(⌊pji⌋,j−1,k)−g(pj−1,j−1,k)]hk(pj)
g ( i , m a x p r i m e , k ) g(i,maxprime,k) g(i,maxprime,k)则表示 ∑ p ≤ i h k ( p ) \sum_{p\le i}h_k(p) ∑p≤ihk(p)
设 S ( n , x ) S(n,x) S(n,x)表示 ∑ w = 1 n [ m i n p > p j ] f ( w ) \sum_{w=1}^n[minp>p_j]f(w) ∑w=1n[minp>pj]f(w)
S ( n , x ) = H ( n ) − H ( p x ) + ∑ p k e ≤ n , k > x [ S ( ⌊ n p e ⌋ , k ) + ( e ≠ 1 ) ] f ( p e ) S(n,x)=H(n)-H(p_x)+\sum_{p_k^e\le n,k>x}[S(\lfloor\frac{n}{p^e}\rfloor,k)+(e\not =1)]f(p^e) S(n,x)=H(n)−H(px)+pke≤n,k>x∑[S(⌊pen⌋,k)+(e=1)]f(pe)
∑ i = 1 n f ( i ) = S ( n , 0 ) \sum_{i=1}^nf(i)=S(n,0) ∑i=1nf(i)=S(n,0)
精细实现,可以做到 O ( n 3 4 log n ) O(\frac{n^\frac{3}{4}}{\log n}) O(lognn43)
P5325 【模板】Min_25筛
#include<cstdio>
#include<cmath>
#define ll long long
using namespace std;
const ll mod=1e9+7,inv2=500000004,inv3=333333336;
ll ksm(ll a,ll b)
{
ll res=1;
while(b)
{
if(b&1) res=res*a%mod;
a=a*a%mod,b>>=1;
}
return res;
}
const int N=2e6+1000;
ll n,w[N+10],sp1[N+10],sp2[N+10],g1[N+10],g2[N+10];int m,cnt,id1[N+10],id2[N+10],prime[N+10];
bool vis[N+10];
int id(ll x)
{
if(x<=m) return id1[x];
else return id2[n/x];
}
void init()
{
m=sqrt(n);
for(int i=2;i<=m;i++)
{
if(!vis[i]) prime[++prime[0]]=i,sp1[prime[0]]=(sp1[prime[0]-1]+i)%mod,sp2[prime[0]]=(sp2[prime[0]-1]+1ll*i*i%mod)%mod;
for(int j=1;j<=prime[0]&&i*prime[j]<=m;j++)
{
vis[i*prime[j]]=true;
if(i%prime[j]==0) break;
}
}
}
ll S(ll n,int x)
{
if(prime[x]>=n) return 0;
int l=id(n);
ll res=(g2[l]+mod-g1[l]+mod-(sp2[x]-sp1[x]))%mod;
for(int i=x+1;i<=prime[0]&&1ll*prime[i]*prime[i]<=n;i++)
{
ll pe=prime[i],tmp;
for(int e=1;pe<=n;e++,pe=pe*prime[i]) tmp=pe%mod,res=(res+tmp*(tmp-1)%mod*(S(n/pe,i)+(e!=1))%mod)%mod;
}
return res;
}
int main()
{
scanf("%lld",&n),init();
for(ll l=1,r;l<=n;l=r+1)
{
r=n/(n/l),w[++cnt]=n/l;
g1[cnt]=(w[cnt]%mod+2)*(w[cnt]%mod+mod-1)%mod*inv2%mod,g2[cnt]=(w[cnt]%mod*(w[cnt]%mod+1)%mod*(w[cnt]*2%mod+1)%mod*inv2%mod*inv3%mod+mod-1)%mod;
if(w[cnt]<=m) id1[w[cnt]]=cnt;
else id2[n/w[cnt]]=cnt;
}
for(int i=1;i<=prime[0];i++)
for(int j=1;j<=cnt&&1ll*prime[i]*prime[i]<=w[j];j++)
g1[j]=(g1[j]+mod-(g1[id(w[j]/prime[i])]+mod-sp1[i-1])%mod*prime[i]%mod)%mod,
g2[j]=(g2[j]+mod-(g2[id(w[j]/prime[i])]+mod-sp2[i-1])%mod*prime[i]%mod*prime[i]%mod)%mod;
printf("%lld",(S(n,0)+1)%mod);
return 0;
}
边栏推荐
- virsh创建/关闭/停止虚拟机常用的几条指令
- Swiftyjson analyse les fichiers json locaux
- adb命令 利用jks文件给apk签名
- [no title] 2022 coal mine safety inspection test questions and online simulation test
- 力扣解法汇总462-最少移动次数使数组元素相等 II
- 力扣解法汇总824-山羊拉丁文
- 力扣解法汇总436-寻找右区间
- ZABBIX notes: 6.0 lts source code installation
- ACL 2022 | 预训练语言模型和图文模型的强强联合
- The unique path of leetcode topic resolution II
猜你喜欢

Graduation design of fire hydrant monitoring system --- thesis (add the most comprehensive hardware circuit design - > driver design - > Alibaba cloud Internet of things construction - > Android App D

Intel case

Ue4\ue5 touch screen touch event: single finger and double finger

(9) Serial port interrupt

Intel case

Start ticwatch2

Bochuang smart sprint technology innovation board: annual revenue of RMB 1.1 billion, book value of accounts receivable of RMB 300million

How to make div 100% page (not screen) height- How to make a div 100% of page (not screen) height?

maya前臺渲染插件mel脚本工具

DbNull if statement - DbNull if statement
随机推荐
ACL2022 | DCSR:一种面向开放域段落检索的句子感知的对比学习方法
Alertmanager alarm configuration
Abaqus中批量对节点施加集中力荷载
Force deduction solution summary 1022- sum of binary numbers from root to leaf
UE4\UE5触摸屏touch事件:单指、双指
Research Report on Chinese psoriasis drug market evaluation and investment direction (2022 Edition)
力扣解法汇总905-按奇偶排序数组
Xcall cluster script (view JPS command)
力扣解法汇总-剑指 Offer II 114. 外星文字典
Is there a female Bluetooth headset suitable for girls? 38 Bluetooth headsets worth getting started
adb命令 利用jks文件给apk签名
el-upload上传文件
2022福建省安全员C证(专职安全员)考试模拟100题及答案
力扣解法汇总713- 乘积小于 K 的子数组
Force deduction solution summary 433- minimum gene change
ADB command uses JKS file to sign apk
没有文笔,大家多多包涵
Force deduction solution summary interview question 01.05 Edit once
String number with special style
Grafana data visualization