筛法瞎写
看见APJifengc在写min25筛发现我这个东西都不会了。赶紧复习了一把去写点题。写一个更一个。
loj6682 梦中的数论
答案显然 \(\frac 12\sum_{i=1}^n\binom{d(i)}{2}=\frac 12\sum_{i=1}^nd^2(i)-d(i)\) 。推导略,以前好像写过。
对了根号的 \(\sum_{i=1}^nd(i)\) 是个傻逼问题,我拿这个诈骗joke3579然后他被骗了。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
const int mod=998244353;
long long n,m,sq,cnt,sum,g[200010],w[200010];
int p[200010],pos1[200010],pos2[200010];
bool v[200010];
long long find(long long x){
return x<sq?pos1[x]:pos2[n/x];
}
void get(int n){
for(int i=2;i<=n;i++){
if(!v[i])p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<=n;j++){
v[i*p[j]]=true;
if(i%p[j]==0)break;
}
}
}
long long s(long long x,long long y){
if(p[y]>=x)return 0;
long long ret=(g[find(x)]-g[find(p[y])]+mod)%mod;
for(int i=y+1;i<=cnt&&1ll*p[i]*p[i]<=x;i++){
for(long long j=p[i],e=1;j<=x;j=j*p[i],e++){
ret=(ret+1ll*(e+1)*(e+1)*(s(x/j,i)+(e>1)))%mod;
}
}
return ret;
}
int main(){
scanf("%lld",&n);
sq=sqrt(n);
get(sq);
for(long long l=1,r;l<=n;l=r+1){
r=n/(n/l);
w[++m]=n/l;
sum=(sum+1ll*w[m]*(r-l+1))%mod;
if(w[m]<sq)pos1[w[m]]=m;
else pos2[n/w[m]]=m;
g[m]=4ll*(w[m]-1)%mod;
}
for(int j=1;j<=cnt;j++){
for(int i=1;i<=m&&1ll*p[j]*p[j]<=w[i];i++){
g[i]=(g[i]-1ll*(g[find(w[i]/p[j])]-g[find(p[j-1])])%mod+mod)%mod;
}
}
long long ans=1ll*(s(n,0)+1-sum+mod)%mod*((mod+1)>>1)%mod;
printf("%lld\n",ans);
return 0;
}
DIVCNTK
我本机开氧气极限数据1.4s。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define int unsigned long long
int n,m,sq,cnt,sum,k,g[2000010],w[2000010];
int p[2000010],pos1[2000010],pos2[2000010];
bool v[2000010];
int find(int x){
return x<sq?pos1[x]:pos2[n/x];
}
void get(int n){
for(int i=2;i<=n;i++){
if(!v[i])p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<=n;j++){
v[i*p[j]]=true;
if(i%p[j]==0)break;
}
}
}
int s(int x,int y){
if(p[y]>=x)return 0;
int ret=1ll*(k+1)*(g[find(x)]-g[find(p[y])]);
for(int i=y+1;i<=cnt&&1ll*p[i]*p[i]<=x;i++){
for(int j=p[i],e=1;j<=x;j=j*p[i],e++){
ret=(ret+(1ll*k*e+1)*(s(x/j,i)+(e>1)));
}
}
return ret;
}
signed main(){
int tim;scanf("%llu",&tim);
get(100000);
while(tim--){
scanf("%llu%llu",&n,&k);
sq=sqrt(n);m=0;
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
w[++m]=n/l;
sum=(sum+1ll*w[m]*(r-l+1));
if(w[m]<sq)pos1[w[m]]=m;
else pos2[n/w[m]]=m;
g[m]=(w[m]-1);
}
for(int j=1;j<=cnt;j++){
for(int i=1;i<=m&&1ll*p[j]*p[j]<=w[i];i++){
g[i]=(g[i]-(g[find(w[i]/p[j])]-g[find(p[j-1])]));
}
}
int ans=(s(n,0)+1);
printf("%llu\n",ans);
}
return 0;
}
DIVCNT3
这题卡你小点,线筛就行了。
#pragma GCC optimize(3)
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define int long long
int n,m,sq,cnt,sum,k,g[1000010],w[1000010];
int p[1000010],pos1[1000010],pos2[1000010];
int low[1000010],pw[1000010],f[1000010];
bool v[1000010];
int find(int x){
return x<sq?pos1[x]:pos2[n/x];
}
void get(int n){
f[1]=1;
for(int i=2;i<=n;i++){
if(!v[i])p[++cnt]=i,low[i]=i,pw[i]=1,f[i]=4;
for(int j=1;j<=cnt&&i*p[j]<=n;j++){
v[i*p[j]]=true;
if(i%p[j]==0){
low[i*p[j]]=low[i]*p[j];
pw[i*p[j]]=pw[i]+1;
f[i*p[j]]=f[i/low[i]]*(3ll*pw[i*p[j]]+1);
break;
}
low[i*p[j]]=p[j];pw[i*p[j]]=1;
f[i*p[j]]=f[i]*f[p[j]];
}
}
for(int i=1;i<=n;i++)f[i]+=f[i-1];
}
int s(int x,int y){
if(p[y]>=x)return 0;
int ret=4ll*(g[find(x)]-g[find(p[y])]);
for(int i=y+1;i<=cnt&&1ll*p[i]*p[i]<=x;i++){
for(int j=p[i],e=1;j<=x;j=j*p[i],e++){
ret=(ret+(3ll*e+1)*(s(x/j,i)+(e!=1)));
}
}
return ret;
}
signed main(){
int tim;scanf("%lld",&tim);
get(500000);
while(tim--){
scanf("%lld",&n);
if(n<=500000){
printf("%lld\n",f[n]);
continue;
}
sq=sqrt(n);m=0;
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
w[++m]=n/l;
sum=(sum+1ll*w[m]*(r-l+1));
if(w[m]<sq)pos1[w[m]]=m;
else pos2[n/w[m]]=m;
g[m]=(w[m]-1);
}
for(int j=1;j<=cnt;j++){
for(int i=1;i<=m&&1ll*p[j]*p[j]<=w[i];i++){
g[i]=(g[i]-(g[find(w[i]/p[j])]-g[find(p[j-1])]));
}
}
int ans=(s(n,0)+1);
printf("%lld\n",ans);
}
return 0;
}
DIVCNT2
今天的成就:爆切了三个水黑题。
卡常。
#pragma GCC optimize(3)
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define int unsigned long long
int n,m,sq,cnt,sum,k,g[2000010],w[2000010];
int p[2000010],pos1[2000010],pos2[2000010];
bool v[2000010];
int find(int x){
return x<sq?pos1[x]:pos2[n/x];
}
void get(int n){
for(int i=2;i<=n;i++){
if(!v[i])p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<=n;j++){
v[i*p[j]]=true;
if(i%p[j]==0)break;
}
}
}
int s(int x,int y){
if(p[y]>=x)return 0;
int ret=3*(g[find(x)]-y);
for(int i=y+1;i<=cnt&&p[i]*p[i]<=x;i++){
for(int j=p[i],e=1;j<=x;j=j*p[i],e++){
ret+=(2*e+1)*(s(x/j,i)+(e!=1));
}
}
return ret;
}
signed main(){
int tim;scanf("%llu",&tim);
get(1000000);
while(tim--){
scanf("%llu",&n);
sq=sqrt(n);m=0;
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
w[++m]=n/l;
sum=(sum+1ll*w[m]*(r-l+1));
if(w[m]<sq)pos1[w[m]]=m;
else pos2[n/w[m]]=m;
g[m]=(w[m]-1);
}
for(int j=1;j<=cnt;j++){
for(int i=1;i<=m&&1ll*p[j]*p[j]<=w[i];i++){
g[i]-=(g[find(w[i]/p[j])]-j+1);
}
}
int ans=(s(n,0)+1);
printf("%llu\n",ans);
}
return 0;
}
DIVCNT1
真的还算筛法吗。
首先你发现这是个傻逼式子:
\[\sum_{i=1}^n\lfloor\frac ni\rfloor
\]
然后你发现这题数据范围 longlong 以内所以过不了。进行一步没有什么卵用的转化:
\[-\lfloor \sqrt n\rfloor^2+\sum_{d=1}^{\lfloor \sqrt n\rfloor}\lfloor \frac nd \rfloor
\]
它还是根号的。因为这玩意的图象关于 \(y=x\) 对称所以成立。
我们用一堆横纵坐标都是整数的向量来拟合这个图象。为了完美拟合,需要 Stern-Brocot Tree (SB树)。
开个单调栈来存一堆向量拟合这个东西。显然单调栈里的元素斜率单调递减。假如说我们当前拟合到 \((x,y)\) 。那么一直执行如下过程:
- 取出栈顶向量,不断加上它,每步累计贡献,直到下一步会走进区域(就是 \(y=\dfrac 1x\) 和坐标轴围成的区域)内。
- 不断弹栈,直到加上栈顶刚好走进区域内。那么我们现在就有了两个向量,一个 \(l\) 加上刚好能进去,一个 \(r\) 刚好进不去。
- 用 SB 树不断二分一个中间向量 \(mid\) 。如果 \(mid\) 能走进去,那么进栈,\(r=mid\) ,接着二分。反之如果 \(r\) 的斜率小于图象在 \(x+x_mid\) 的斜率,那么之后二分的所有向量都不能走进去,直接停。反之 \(l=mid\) 继续。
边界是 \(y\le \sqrt[3]{n}\) ,直接暴力。总复杂度 \(O(\sqrt[3]n\log n)\) 。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define int __int128
int read(){
int x=0;char ch=getchar();
while(!isdigit(ch))ch=getchar();
while(isdigit(ch))x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
return x;
}
void print(int x){
if(x>9)print(x/10);
putchar(x%10+'0');
}
struct node{
int x,y;
node operator+(const node &s)const{
return {x+s.x,y+s.y};
}
}stk[1000010],l,r,p;
int ans,n,top,sq,cb;
signed main(){
int tim=read();
while(tim--){
n=read();ans=top=0;
sq=sqrt(n);cb=cbrt(n);
p={n/sq,sq+1};
stk[++top]={1,0};stk[++top]={1,1};
while(1){
l=stk[top];top--;
while(n<(p.x+l.x)*(p.y-l.y)){
ans+=p.x*l.y+(l.x-1)*(l.y+1)/2;
p.x+=l.x;p.y-=l.y;
}
if(p.y<=cb)break;
r=stk[top];
while(1){
node mid=l+r;
if(n<(p.x+mid.x)*(p.y-mid.y)){
r=mid;stk[++top]=r;
}
else if(n*r.x<=1ll*(p.x+mid.x)*(p.x+mid.x)*r.y)break;
else l=mid;
}
}
for(int i=1;i<p.y;i++)ans+=n/i;
ans=ans*2-sq*sq;
print(ans);putchar('\n');
}
return 0;
}
快踩