luogu P5325 Min_25筛
LINK：Min_25筛
新版感觉有点鬼畜 而且旧版的也够用了至少.
这个并不算很简单也不算很困难的知识点 学起来还是很麻烦的。
(误入了很多dalao的blog 说的云里雾里的 甚是懵逼 这里推荐几个blog一起看 能看出很多门道
网上资源辣么多 我自然也不会去写一个非常正常的学习笔记辣.. 只会写几个容易疑惑的地方。
注意 学会 和会写代码是两码事 因为代码中有一些细节需要细细揣摩。
关于g数组的求出 其转移静下心来理解还是可以看懂的这里不再赘述。
注意 为了方便\(f(1)\)最后考虑.
设 \(maxx=\sqrt n\)
\(g_{n,j}=\sum_{i=1}^{n}f(i)[i\in P|min_p(i)>p_j]\)
关于g数组 第二维的范围显然只有\(maxx\)的大小。
因为大于maxx的转移都是承接上一个的状态 所以不必要求 或者可以理解欧拉筛的思想 筛到maxx大小的质数之后 所有数字都被筛完了。
考虑第一维 这是根据 我们后面的S组数来定的 但是注意到后面的S每次只会用到\(\frac{n}{p_i}\)
虽然我们不知道有多大但是数量级还是 是固定的 因为\(\frac{n}{i}\)最多只有\(maxx\)数量级 注意 真正的大小应该是\(2\cdot maxx\)
从状态转移方程 来看 第二维可以滚动。看似状态数量很多 实际上关于 \(p_j\cdot p_j>n\)的那些状态都是没有必要的 所以状态数量在一个可控范围。
一个笔者还未搞懂的问题 以后可能会完成回答：实际上在求g数组的时候利用到了完全积性函数的性质 当f不是完全积性函数的时候
不过很多blog上说可以这样做 且 只有\(g_{n,|P|}\)的位置上的值有效 这个地方还不太清楚为什么。
下面关于求答案的部分。
设\(S(n,j)=\sum_{i=1}^{n}f(i)[min_p(i)>=p_j]\)
容易得到\(S(n,1)+f(1)\) 且有\(S(i,j)=g_{i,|P|}-sum_{j-1}+\sum_{k>=j}\sum_{e}f(p^e)(S(\frac{i}{p^e},k+1)+[e\neq 1])\)
然后就没了 复杂度的大概就是1s 1e10,3s 1e11的样子。
code
//#include<bits/stdc++.h>
#include<iostream>
#include<iomanip>
#include<cstdio>
#include<ctime>
#include<cstdlib>
#include<cctype>
#include<cstring>
#include<cmath>
#include<string>
#include<utility>
#include<queue>
#include<vector>
#include<algorithm>
#include<deque>
#include<stack>
#include<list>
#include<bitset>
#include<set>
#include<map>
#define INF 1000000000000000000ll
#define rep(p,n,i) for(int i=p;i<=n;++i)
#define fep(n,p,i) for(int i=n;i>=p;--i)
#define vep(p,n,i) for(int i=p;i<n;++i)
#define db double
#define get(x) x=read()
#define put(x) printf("%d\n",x)
#define pb push_back
#define ll long long
#define db double
#define mod 1000000007
#define putl(x) printf("%lld\n",x)
using namespace std;
char *fs,*ft,buf[1<<15];
inline char getc()
{
return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;
}
inline ll read()
{
ll x=0,f=1;char ch=getc();
while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getc();}
while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getc();}
return x*f;
}
const int MAXN=200010,N=5,inv3=333333336;
ll n;
int maxx,cnt,top;
ll g1[MAXN],g2[MAXN],w[MAXN];
int sum1[MAXN],sum2[MAXN],id1[MAXN],id2[MAXN];
int v[MAXN],p[MAXN];
inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline void prepare()
{
rep(2,maxx,i)
{
if(!v[i])
{
v[i]=p[++top]=i;
sum1[top]=add(sum1[top-1],i);
sum2[top]=add(sum2[top-1],(ll)i*i%mod);
}
rep(1,top,j)
{
if(p[j]>maxx/i)break;
v[p[j]*i]=p[j];
if(v[i]==p[j])break;
}
}
}
inline void calc()
{
for(ll i=1,j;i<=n;i=j+1)
{
w[++cnt]=n/i;
j=n/w[cnt];
g1[cnt]=w[cnt]%mod;
g2[cnt]=g1[cnt]*(g1[cnt]+1)/2%mod*(2*g1[cnt]+1)%mod*inv3%mod-1;//平方和.
g1[cnt]=(1+g1[cnt])*g1[cnt]/2%mod-1;
if(w[cnt]<=maxx)id1[w[cnt]]=cnt;
else id2[j]=cnt;
}
rep(1,top,i)
{
rep(1,cnt,j)
{
if((ll)p[i]*p[i]>w[j])break;
ll cc=w[j]/p[i],k;
if(cc<=maxx)k=id1[cc];else k=id2[n/cc];
g1[j]=(g1[j]-p[i]*(g1[k]-sum1[i-1]))%mod;
g2[j]=(g2[j]-p[i]*(g2[k]-sum2[i-1])%mod*p[i])%mod;
}
}
}
inline int S(ll x,int y)
{
if(x<=1||p[y]>x)return 0;
int k=x<=maxx?id1[x]:id2[n/x];
ll ans=(g2[k]-g1[k]-sum2[y-1]+sum1[y-1])%mod;
for(int i=y;i<=top;++i)
{
if((ll)p[i]*p[i]>x)break;
ll pe=p[i],pp=pe*pe;
for(int e=1;pp<=x;++e,pe=pp,pp*=p[i])
{
ans=(ans+pe%mod*((pe-1)%mod)%mod*S(x/pe,i+1)+pp%mod*((pp-1)%mod)%mod)%mod;
}
}
return ans;
}
signed main()
{
freopen("1.in","r",stdin);
get(n);maxx=(int)sqrt(n*1.0)+1;
prepare();calc();put((S(n,1)+1+mod)%mod);
return 0;
}
