题解 SP34096 【DIVCNTK - Counting Divisors (general)】

· · 题解

DIVCNT2是杜教筛,DIVCNT3是洲阁筛……

好麻烦啊，不想写qaq……

怎么办啊？

写DIVCNTK然后把k改成2和3就行了啊

使用下面的算法只需要把代码中的k改成2和3就能通过divcnt2,3了，所以这是个3倍经验……

好了让我们看一下这题让我们求什么?

\sum_{i=1}^{n}\sigma _{0}(i^k)

其中\sigma_{0}是除数函数

那么我们不难发现函数f(x)=\sigma_{0}(x^k)满足这么几条性质

f(1)=1 f(p)=k+1 f(p^c)=kc+1 f(ab)=f(a)(b),gcd(a,b)=1

解释一下就是我们发现如果p是一个质数那么p^c只有p的0,1,2,3...c次方是它的约数，所以\sigma_{0}(p^c)=c+1然后就可以十分轻松的推出上面的式子了

看起来好像没啥性质啊，怎么做这题呢？

答案是一个被称为Min_25筛的黑科技，这个东西的复杂度在今年的集训队论文当中刚刚被证明，令人惊讶的是，它的复杂度居然是

O(\frac{n^{0.75}}{logn}+n^{1-\epsilon})

其实这个复杂度记号并不标准

那个O(n^{1-\epsilon})被称为“近线性复杂度”换句话说当n趋于无穷的时候这个东西的复杂度将会无限逼近于O(n)并且大于任意的O(n^a) (当然a是一个小于1的数字)

然后这并没有什么卵用，这东西的复杂度趋于O(n)是真的，但是也没说趋近速度如何，现实就是这东西在洲阁筛跑的动的数据范围内爆踩洲阁,如果对常数要求不严格甚至还能水一些杜教筛的题……,所以在10^{13}之内使用这个算法是完全安全的

好了让我们来看一下这个算法能筛什么积性函数f(x)的前缀和吧

首先这个函数必须是一个积性函数，

其次我们要求f(p)=f^{'}(p)其中f^{'}(p)是一个可以快速求前缀和的完全积性函数，不然做不了

最后一点就是f(p^c)可以快速计算

那么Min_25筛分为两部分，第一部分是求质数点值，而第二部分就是大力搜索

让我们先来介绍第一部分

Part1:求所有质数点值的和

众所周知给定一个n，\lfloor \frac{n}{i} \rfloor的取值恰好有2\lfloor \sqrt{n} \rfloor种，那么我们在第一部分要做的事情就是

求出一个函数Sp(n),它表示所有小于等于n且是质数的i的f(i)之和，并且对于\lfloor \frac{n}{i} \rfloor的每一个取值x都求出Sp(x)的值

那么我们怎么做呢,容易发现当i是质数的时候f(i)=f^{'}(i)所以我们只需要求出所有是质数的f^{'}(x)的和就行了

那么我们可以设计一个dp,不妨令P_{j}表示第j大的质数，那么我们定义函数S(n,j)表示满足以下两种限制之一的f^{'}(i)的和

1.i是一个小于等于n的质数

2.i是一个小于等于n的非质数，并且i的最小质因子大于P_{j}

那么我们发现一个比较显然的事实是

S(n,j)=S(n,j-1),P_{j}^2>n

为什么呢？因为我们发现对于第二类i来讲，最小的i应该是P_{j+1}^{2}如果连这个数字都大于n的话这部分就为0了，因此我们就可以在S(n,j),S(n,j-1)之间画等号了

如果很不巧，我们的P_{j}^2 \leq n那么我们有这样的转移式子

为了式子的美观我直接省略了下取整符号

S(n,j)=S(n,j-1)-f^{'}(P_{j})(S(\frac{n}{P_{j}},j-1)-\sum_{i=1}^{j-1}f^{'}(P_{i}))

为什么这个式子是对的呢？

我们发现S(n,j)和S(n,j)差的数是这样的数

i是一个小于等于n的非质数并且i的最小质因子恰好等于P_{j}

那么由于f^{'}(x)是一个完全积性函数，所以我们可以将所有最小质因子是P_{j}的数字除一个P_{j}这样就递归到了S(\frac{n}{P_{j}},j-1)的子问题上，不过有一个问题是S(\frac{n}{P_{j}},j-1)还包含小于P_{j}的质数，因此我们要把这些东西减去

那么我们预处理一下\sum_{i=1}^{n}f^{'}(P_{i})就可以O(1)的实现递推了

直接暴力的实现这个dp是O(\frac{n}{logn})的会非常非常的慢

但是我们发现一件事情是如果我们用滚动数组实现这个dp的话，当P_{j}^2>n的时候我们可以直接跳过这个状态而不去更新它而且更妙的是如果我们将所有需要dp的x排好序的话,我们还可以拿两根双指针维护一下还需要更新的x值，如此这般dp的话根据积分估计我们的复杂度就被优化到了

O(\frac{n^{0.75}}{logn})

好了这样dp到最后的话我们就求出了Sp(n)了

边界条件S(n,0)=\sum_{i=2}^{n}f^{'}(x)由于这个函数能快速求前缀和因此可以算出来

注意一点不要拿map或者hash表来存，而是对于小于\sqrt{n}的值直接拿数组存而大于\sqrt{n}的值将他的dp值存在\frac{n}{x}这个下标上，可以证明这样存储是没有冲突的

对于这题来讲我们发现f'(p)并不完全积性,不过我们可以用上面的技巧dp出来每个x以内有几个质数，这样就可以把dp值乘上一个(K+1)就能算Sp(n)了

Part2 筛积性函数前缀和

怎么筛呢？

暴力dfs，对你没看错，就是暴 力 搜

我们设G(n,j)表示所有小于等于n并且最小质因子大于等于P_{j}的i的f(i)之和,注意我们并不考虑1

那么怎么计算G(n,j)呢？答案是暴力，将贡献拆成两部分，一部分是质数的贡献而另一部分是合数的贡献

质数的贡献显然相当好算，就是Sp(n)-\sum_{i=1}^{j-1}f^{'}(P_{i}),因为需要保证最小质因子大于等于P_{j}所以要减去一些东西

而合数部分的贡献就是暴力递归,大力枚举最小质因子的值以及次幂,我们可以得到这样一个式子

G(n,j)=Sp(n)-\sum_{i=1}^{j-1}f^{'}(P_{i})+\sum_{j \leq k}\sum_{e}f(P_{k}^{e})(G(\frac{n}{P_{k}^{e}},k+1)+[e>1])

解释一下就是把最小质因子全部除掉然后直接递归下去算，由于G的定义不包含1因此需要加上[e>1]这个修正因子

好了接下来是更加玄学的地方，实现这个东西，不 需 要 记 忆 化

直接暴力搜索即可，然后也不难理解为什么这个东西会变成近线性复杂度了，但是这东西就是快，你没办法……

然后对于刚才的式子我们直接把f^{'}(P)和f(P_{k}^{e})带进具体的表达式算就行了

备注:把这份代码的K改成2和3可以直接通过DIVCNT2和3

上代码~

#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;const int N=1e5+1000;typedef unsigned long long ll;
int zhi[N];int ct;bool book[N];int lim;ll n;ll K;ll q[N];int hd;int T;
ll dp1[N];ll dp2[N];ll pre[N];
inline ll solve(ll p,int q)//灵魂搜索代码 
{
    if(p<zhi[q])return 0;ll ret=((p<=lim)?dp1[p]:dp2[n/p])-pre[q-1];
    if(q==ct+1)return (p<=(ll)zhi[ct])?0:ret;ll pw;int t;
    for(int k=q;k<=ct&&(ll)zhi[k]*zhi[k]<=p;k++)
        for(pw=zhi[k],t=1;pw<=p;pw*=zhi[k],t++) 
            ret+=(K*t+1)*(solve(p/pw,k+1)+(t>1));return ret;    
}
inline void calc1()//dp预处理 
{
    for(int i=1;i<=lim;i++)dp1[i]=i-1;for(int i=1;i<=hd;i++)dp2[i]=q[i]-1;
    for(int j=1,t1=hd;j<=ct;j++)
    {
        ll lt=(ll)zhi[j]*zhi[j];while(lt>q[t1]&&t1>=1)t1--;
        for(int i=1;i<=t1;i++)
        {ll fv=q[i]/zhi[j];fv=(fv<=lim)?dp1[fv]:dp2[n/fv];dp2[i]-=fv-j+1;}
        for(int i=lim;i>=lt;i--)dp1[i]-=dp1[i/zhi[j]]-j+1;
    }for(int i=1;i<=lim;i++)dp1[i]=(K+1)*dp1[i];
    for(int i=1;i<=hd;i++)dp2[i]=(K+1)*dp2[i];
}
inline void solve()
{
    scanf("%llu%llu",&n,&K);lim=sqrt(n);for(ct=1;(ll)zhi[ct]<=lim;ct++);ct--;hd=0;
    for(ll i=1,r;i<=n;i=r+1){r=n/i;if(r>lim)q[++hd]=r,r=n/r;else break;}
    for(int i=1;i<=ct;i++)pre[i]=(K+1)*i;calc1();
    printf("%llu\n",solve(n,1)+1);
}
int main()
{
    for(int i=2;i<=N-10;i++)//线性筛素数 
    {
        if(!book[i])zhi[++ct]=i;
        for(int j=1;j<=ct&&(ll)i*zhi[j]<=N-10;j++)
        {book[i*zhi[j]]=true;if(i%zhi[j]==0)break;}
    }scanf("%d",&T);for(int z=1;z<=T;z++)solve();return 0;//拜拜程序~ 
}