这些各种乱七八糟的筛法真难懂……
首先Min_25筛的基本思想就是在不停的枚举最小质因子。
zzt的论文根本看不懂。
过程是这样的,既然F是个低阶多项式,那么先要求:
这个怎么算,首先将g_k(n)初始化为1~n的k次方和,然后减去不合法的部分。
即我们从小到大枚举质数p,定义g_k(n)为
。那么随着p的增大,g会越来越小。每次会减小哪些呢?显然是那些以p做最小质因子的情况,也就是x=tp,t的最小质因子不小于p。这些x的和就是
,但是x还有可能是小于p的质数,因此答案准确的说应当是
。也就是要从g_k(n)减去这个东西。显然n要i从大到小枚举,并且当
的时候就要停止了(因为不可能减掉任何东西了)。
实现起来是这样的,维护两个数组s[x]=g_k(x),l[x]=g_k(n/x),这样x的范围不超过根号n,然后转移的时候去分类讨论一下即可。
(这里有个问题,对于n/x相同的x,我们都存到一起去了,他们的答案不会有差别么?不是很清楚)
然后考虑维护答案。首先记S(x,y)表示1~x中,最小质因子不小于prime[y]的数字的F的和。那么答案显然就是S(n,1)+F(1)。
考虑去计算S,答案分成两部分,第一部分是i是质数的情况,这个直接用刚刚求出的g函数做一下就可以了;后半部分考虑还是枚举i的最小质因子prime[z]>=prime[y],那么这个
,t的最小质因子>prime[y](注意是大于),这样枚举z和e>=1使得
,然后答案加上
(后面那个是去计算t=1的情况,由于F(prime[z])已经在之前的g函数中处理过了,所以这里就不用算了)。注意到这里至少是
,因此只要枚举到根号x的质数即可。这里一个小细节是为啥没有加上
(e取到最大),这是因为prime[z]的e+2次是超过了n的,那么prime[z]的e+1次再乘以后面某个大于prime[z]的质数也就爆炸了)。
……终于说完了。这东西复杂度为啥是对的我也不知道,反正跑的还挺快就是了。
代码是divcntk的代码。
#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<cmath>
#define ull unsigned long long
#define MXS 100010
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
int notp[MXS],pri[MXS];
ull s[MXS],l[MXS];
inline int my_sqrt(ull n)
{
int x=(int)sqrt(n);
while(1ull*x*x<n) x++;
while(1ull*x*x>n) x--;
return x;
}
inline int prelude(int n)
{
notp[1]=1;int pc=0;
for(int i=2;i<=n;i++) if(!notp[i])
{
pri[++pc]=i;
for(int j=i;(ull)j*i<=(ull)n;j++)
notp[i*j]=1;
}
return pri[++pc]=n+1;
}
inline int get_g(ull n,int sn,ull k)
{
for(int i=1;i<=sn;i++) s[i]=i,l[i]=n/i;
for(int p=2;p<=sn;p++)
{
if(notp[p]) continue;ull r=(ull)p*p,v=s[p-1];
int ed1=sn/p,ed2=(int)min(n/r,(ull)sn);
for(int i=1;i<=ed1;i++) l[i]-=l[i*p]-v;
for(int i=ed1+1;i<=ed2;i++) l[i]-=s[n/i/p]-v;
for(int i=sn;(ull)i>=r;i--) s[i]-=s[i/p]-v;
}
for(int i=1;i<=sn;i++) s[i]*=(k+1),l[i]*=(k+1);
return 0;
}
inline ull F(ull n,ull x,int y,int sn,ull k)
{
if((ull)pri[y]>x) return 0;ull res=-s[pri[y]-1];
if(x<=(ull)sn) res+=s[x];else res+=l[n/x];
for(int i=y;(ull)pri[i]*pri[i]<=x;i++)
{
ull v=pri[i];
for(int j=1;v*pri[i]<=x;j++,v*=pri[i])
res+=F(n,x/v,i+1,sn,k)*(j*k+1)+(j+1)*k+1;
}
return res;
}
int main()
{
ull n,k;int sn;scanf("%llu%llu",&n,&k);
sn=my_sqrt(n),prelude(sn),get_g(n,sn,k);
return !printf("%llu\n",F(n,n,1,sn,k)+1);
}