[学习笔记] SPOJ DIVCNTK - Min_25筛

这些各种乱七八糟的筛法真难懂……
首先Min_25筛的基本思想就是在不停的枚举最小质因子。
zzt的论文根本看不懂。

下文所有除法为了方便都是下取整。

过程是这样的，我们先求这样一个东西：
$$g(n)=\sum _{i=1}^n[\mathrm{i}\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}]F(i)$$
通常题目会保证$F(p)$是一个系数固定的低次多项式，假设是k次，那么我们就是要求：
$$g_k(n)=\sum _{i=1}^n[\mathrm{i}\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}]i^k$$
这个怎么算，首先将$g(n)$初始化为${\sum }_{i=1}^n{i}^k$，然后减去不合法的部分（即$i$不是质数的情况）。
即我们从小到大枚举质数$p$，并且定义：
$$g_k(n)=\sum _{i=1}^n[i是质数或者i的最小质因子大于p]i^k$$
那么随着$p$的增大，$g_k(n)$会越来越小。每次会减小哪些呢？显然是那些以$p$做最小质因子的情况，也就是$x=tp$，其中$t$的最小质因子不小于$p$。这些$x$的和就是$p^k{g}_k(n/p)$，但是这其中$x$还有可能是小于$p$的质数，需要减去，因此答案准确的说应当是$p^k(g_k(n/p)-g_k(p-1))$。这就是要从$g_k(n)$减去这个东西。

下文我们使用的时候，设$s=\lfloor \sqrt{N}\rfloor ,x\le s$，要么我们要问$g_k(x)$，要么问$g_k(N/x)$。（$N$和$n$是不一样的，$N$是输入的，全局固定的）。

我们可以注意到，若$a_1,a_2\le s$并且$n/a_1=n/a_2$，那么可以说明$a_1=a_2$。
因此我们用$A(x)=g_k(x),B(x)=g_k(N/x),1\le x\le s$。（代码里面是s和l）。

显然上文过程$t$是可能包含$p$作为最小质因子的，因此$n$要从大到小枚举，并且对于$n<p^2$的$n$是没有讨论必要的。
转移（注意转移顺序）：
$$\begin{gathered} B(n)=g_k(N/n),-=p^k(g_k(n/ip)-g(p-1)),n\le \min (N/p^2,s) \\ \therefore \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}n\le s/p,B(n)-=p^k(B(np)-A(p-1)) \\ \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}B(n)-=p^k(A(n/(ip))-A(p-1)) \\ A(n)=g_k(n),-=p^k(g_k(n/p)-g_k(p-1))=p^k(A(n/p)-A(p-1)) \end{gathered}$$
（上面虽然看上去好长其实很好推，之不过我把所有细节写下来了而已）

然后考虑维护答案。
首先记$S(x,y)$表示$1～x$中，最小质因子不小于$prime[y]$的数字的$F$的和。那么答案显然就是$S(n,1)+F(1)$。

考虑去计算$S(x,y)$，答案分成两部分，

第一部分是$i$是质数的情况，这个直接用刚刚求出的$g$函数做一下就可以了；

后半部分考虑枚举$i$的最小质因子$prime[z]>=prime[y]$，那么这个$i=prime[z{]}^e\times t$，其中$t$的最小质因子$>prime[z]$（注意是大于）。

这样枚举$z$和$e\ge 1$使得$prime[z{]}^{e+1}\le x$，然后答案加上$S(x/(prime[z{]}^e),z+1)\times F(prime[z{]}^e)+F(prime[z{]}^{e+1})$（后面那个是去计算$t=1$的情况，由于$F(prime[z])$已经在之前的$g$函数中处理过了，所以这里就不用算了）。

注意到这里至少是$prime[z{]}^2\le x$，因此只要枚举到根号$x$的质数即可。
这里一个小细节是为啥没有加上$S(x/(prime[z{]}^{e+1}),z+1)\times F(prime[z{]}^{e+1})$（$e$取到最大），这是因为$prime[z]$的$e+2$次是超过了$n$的，那么$prime[z]$的$e+1$次再乘以后面某个大于$prime[z]$的质数也就爆炸了）。
……终于说完了。这东西复杂度为啥是对的我也不知道，反正跑的还挺快就是了。
代码是divcntk的代码，满足$F(p)=k+1$（是个零次多项式）。

#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);
}