题解 SP20174 【DIVCNT3 - Counting Divisors (cube)】

· · 题解

首先这题是一道洲阁筛板子题，就当做是一篇洲阁筛讲解吧不要跟我提啥min25，另外部分内容参考自__debug大神的博客在此%%%OrzstO%%%一下。

洲阁筛和杜教筛一样，是用来在低于O(n)的时间内求一类积性函数的前缀和的。然后这些函数F(x)必须在当x是质数的时候是容易求的，emm……或者应该这样说，当x是质数的时候，F(x)是关于x的一个常数次的多项式。洲阁筛的关键就是利用了在这些质数位置的函数值，然后再利用积性函数的性质钦定这些质数是mindiv或者maxdiv从而把答案递推出来。

好了，这个题让我们求这个：

\sum_{i=1}^nσ_0(i^3)

然后我们不妨以\sqrt n为界限把这n个F(x)=σ_0(x^3)的求和给他分开，就是这n个数里的每个数都最多只会有一个大于\sqrt n的质因子吧，如果真的存在的话把这个大质数除去剩下的部分显然是不超过\sqrt n的，我们不妨在[1,\sqrt n]内枚举这个部分，那么由于F(x)=σ_0(x^3)是积性函数，当x\in[1,\sqrt n]或者x存在大于\sqrt n的质因子，这些F(x)的和可以表示为：

\sum_{i=1}^nF(i)(1+\sum_{j=\sqrt n+1}^{\lfloor\frac{n}{i}\rfloor}[j\text{为质数}]F(j))……(1)

然后剩下的数既大于\sqrt n，又没有超过\sqrt n的质因子，我们把他们的F(x)单独求个和：

\sum_{i=\sqrt n+1}^n[maxdiv(i)<=\sqrt n]F(i)……(2)

对于①式因为i<=\sqrt n我们可以把他们的F预处理出来，然后对于每个i\in[1,\sqrt n]我们需要求出[\sqrt n+1,\lfloor\frac n i \rfloor]内所有质数的F(x)之和，在最开头我们说了F(x)需要在当x是质数的时候是一个关于x的常数次多项式，实际上我们发现F(x)=σ_0(x^3)当x是质数的时候恒等于4……于是我们的问题其实就是求这段区间内的质数个数。

首先这个区间内的质数它必须是大于\sqrt n的，这是一个很厉害的条件，那么我们就可以把这个问题等价转化为在[1,\lfloor\frac n i \rfloor]内有多少数和[1,\sqrt n]内的质数互质，显然其他合数是达不到这个条件的。

于是就可以有这样一个递推的思路（我发现这种把问题递推出来的思路我自己肯定想不出来……），就是设g(i,j)为在[1,j]内与P_{1…i}（P是质数数组）的质数都互质的数的个数。我们肯定要考虑从g(i-1,j)->g(i,j)的转移，即新加了一个质数P_i，显然随着i的不断增大满足条件的数会不断变少，这之前g(i-1,j)是仅与前i-1个质数互质的数，这其中势必会有P_i的倍数，显然P_i是这些数的mindiv，那么我们要把这些数减掉（也就是“筛”），已经确定了他们的mindiv是P_i，那么剩下的部分显然就在[1,\lfloor\frac j {P_i} \rfloor]内，也不能出现小于P_i的质因子（等于P_i也没事），那么就是g(i-1,\lfloor\frac j {P_i} \rfloor)，于是g的递推式就是：

g(i,j)=g(i-1,j)-g(i-1,\lfloor\frac j {P_i}\rfloor)

注意g(0,j)就是没有互质限制的数，显然就是j了，这是边界问题。

然后我们来考虑他的优化。先明确一下我们的目标是对每个i\in [1,\sqrt n]都求一个g(m,\lfloor\frac n i\rfloor)（我们令P_m为<=\sqrt n的最大的质数，注意是<=而不是<我被这个坑了好久……），然后我们用递归的那种方式去想，发现j一开始就是n除一个整数，然后在递归的过程中j将会不断的除一些整数，根据等式\lfloor\frac{\lfloor\frac n {i}\rfloor}{j}\rfloor=\lfloor\frac n {ij}\rfloor可以发现j实际上是n不断除一个整数除出来的，我们知道n除[1,n]内的整数只有O(\sqrt n)级别个不同的结果，再多了也就是0，所以我们断言j只有O(\sqrt n)个，我们可以先把这些j求出来用哈希表离散化一下再去递推。

然而这样做复杂度仍然是暴力，但是我们可以接着用一开头就用过的根号思想，就是在[1,j]内的数最多只会有一个超过\sqrt j的质因子，那么当P_i^2>j的时候我们不是要把P_i的倍数筛掉吗？显然最多只有一个数会被筛掉，而且i越大P_i也是越大的，那么当我们在递推的时候发现P_i^2>j的时候就在以后不再枚举这个j了，而是记录一下最后访问他的位置i_{last}与g(i_{last},j)，当我们以后要访问一个g(i,j)的时候就用g(i_{last},j)减去P[i_{last}+1,i]内有多少个<=j的质数（这个比较坑因为这些质数不一定都在[1,j]内，可能就没有这个质数不需要筛掉），这个乱搞一下还是很好求的。

然后递推出这个g数组的复杂度据说是O(\frac{n^{\frac 3 4}}{logn})的，下面除一个logn大概是因为i是枚举的质数密度会除个logn……算了不口胡了反正我不会证……感性理解一下就好

然后我们就成功地解决了①式。这个②式我们可以采用和上面一样的递推的做法，显然我们可以把这个式子前缀和相减一下，那么就是求一段区间内maxdiv不超过\sqrt n的所有数的F(x)之和，于是我们照葫芦画瓢的设f(i,j)为[1,j]中maxdiv不超过P_i的数的\sum F(x)，边界就是f(0,j)=1显然一个质因子都不能有只有1能够满足。这个的转移方式跟上面的g正好相反，我们考虑f(i-1,j)->f(i,j)的转移时是考虑在当前质数集合里添加一个P_i，我们把maxdiv=P_i的所有数都添加进来，但是我们在这里就不能让去掉P_i之后剩下的那部分带有P_i，因为那个也要从f(i-1,?)转移过来，所以我们还要枚举一下P_i在这个数里的次数，这一点跟g还是有所不同的。于是再根据积性函数的性质（所以求积性函数的和还是相对较容易的……）f的递推式就长成这样：

f(i,j)=f(i-1,j)+\sum_{q=1}^{P_i^q<=j}F(P_i^q)f(i-1,\lfloor\frac{j}{P_i^q}\rfloor)

然后这个式子应该也可以用和g一样的优化？当时我在学洲阁筛的时候说f需要把i的意义改为后i个质数然后就可以用上面那个优化了，现在突然发现好像直接套用g的方法好像也是可以的……不过并没有试过qwq……把i倒过来做就是让f(i,j)变成不超过\sqrt n的后i个质数，我们在递推的时候只枚举P_i^2<=j的j，当他第一次被枚举到就把没把他枚举到的质数添加进去（在[1,j]中只出现一次），然后复杂度也是O(\frac{n^{\frac 3 4}}{logn})的。

然后这题就做完了……

然后你在理论上分析完这题可做，然后一写代码发现越写越恶心……因为全都是一堆细节问题……建议在写代码之前在一张纸上把整个过程一点不漏的写出来。然后下面是一些实现上的问题。

1.我们对每次询问需要的f和g的i都是m就是不超过\sqrt n的最大质数的下标，那么我们可以先把表打出来，并且i显然是可以用滚动数组滚掉的。

2.我们要把j用哈希表离散化……不要引入带log的二分或map……对了不要忘记清空哈希表。

3.在求一段区间内的质数有多少个小于x的时候，可以不用二分，因为这都是\sqrt n级别的，所以我们在最开始线性筛预处理质数的时候可以对每个数记一下不超过他的最大的质数的下标是啥就可以直接O(1)算了。

4.n比较小的时候可以直接用前缀和，这样会快一点。

另外不要高估了这题的数据……虽然有1e4的询问但是这个复杂度还是可以过的……数据里不会全是极端数据。

上代码~

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define int unsigned long long
#define p 19260817
using namespace std;
typedef unsigned long long ll;
namespace ywy{
    typedef struct _n{
        ll data;int dp;int nxt;
    }
    node;node memchi[2000001];
    int gn=1;int heads[19260818];
    inline void insert(ll data,int dp){
        memchi[gn].data=data;
        memchi[gn].dp=dp;
        memchi[gn].nxt=heads[data%p];
        heads[data%p]=gn;gn++;
    }
    inline int get(ll data){
        for(register int i=heads[data%p];i;i=memchi[i].nxt){
            if(memchi[i].data==data)return(memchi[i].dp);
        }return(0);
    }
    inline ll get(){
        ll n=0;char c;while((c=getchar())||23333){
            if(c>='0'&&c<='9')break;if(c=='-')goto s;
        }n=c-'0';while((c=getchar())||23333){
            if(c>='0'&&c<='9')n=n*10+c-'0';else return(n);
        }s:while((c=getchar())||23333){
            if(c>='0'&&c<='9')n=n*10-c+'0';else return(n);
        }return(0);
    }
    int nearby[400001];
    void iprint(ll num){
        if(num>=10)iprint(num/10);
        putchar(num%10+'0');
    }
    int prime[400001];
    unsigned char bv[400001];
    int ptr=1;
    int mindiv[400001];
    ll sig[400001],sums[400001];
    inline void prim(int n){
        for(register int i=2;i<=n;i++){
            if(!bv[i]){
                nearby[i]=ptr;
                prime[ptr]=i;ptr++;
                mindiv[i]=i;
            }for(register int j=1;j<ptr;j++){
                int cjr=i*prime[j];
                if(cjr>n)break;
                mindiv[cjr]=prime[j];
                bv[cjr]=1;
                if(i%prime[j]==0)break;
            }
        }for(register int i=1;i<=n;i++){
            if(!nearby[i])nearby[i]=nearby[i-1];
            ll ans=1;int tmp=i;
            while(tmp!=1){
                int d=mindiv[tmp];int k=1;
                while(mindiv[tmp]==d)k+=3,tmp/=d;ans*=k;
            }
            sig[i]=ans;
            sums[i]=sums[i-1]+ans;
        }
    }
    ll f[2][4000001],g[2][4000001];
    inline int erfen(ll n){
        int ans=0,l=1,r=ptr-1;
        while(l<=r){
            int mid=(l+r)>>1;
            if((ll)prime[mid]*(ll)prime[mid]<=n)ans=mid,l=mid+1;
            else r=mid-1;
        }return(ans);
    }
    ll ints[4000001];
    ll n;ll sqr;int sqrp;int pos=0;
    int i0[1000001];ll rg[1000001];
    inline int how(int l,int r,ll x){
        if(l>r)return(0);
        if(prime[r]<=x)return(r-l+1);
        if(prime[l]>x)return(0);
        return(nearby[x]-l+1);
    }
    inline void getg(){
        for(register int i=0;i<pos;i++)g[0][i]=ints[i],i0[i]=0,rg[i]=0;
        for(register int i=1;i<=sqrp;i++){
            for(register int j=pos-1;~j&&!i0[j];j--){
                ll ans=g[(i-1)&1][j];
                if(prime[i]>ints[j]/prime[i]){
                    i0[j]=i;ans-=(ints[j]>=prime[i]);rg[j]=ans;
                }else{
                    int dst=get(ints[j]/prime[i]);
                    if(i0[dst])ans-=(rg[dst]-how(i0[dst]+1,i-1,ints[j]/prime[i]));
                    else ans-=g[(i-1)&1][dst];
                }g[i&1][j]=ans;
            }
        }
    }
    int fi0[4000001];
    inline void getf(){
        for(register int i=0;i<pos;i++)fi0[i]=0,f[0][i]=f[1][i]=0;
        for(register int i=1;i<=sqrp;i++){
            for(register int j=pos-1;~j&&prime[sqrp-i+1]<=ints[j]/prime[sqrp-i+1];j--){
                ll ans=f[(i-1)&1][j];
                if(!fi0[j])fi0[j]=i,ans=1ll+4ll*(ll)how(sqrp-i+2,sqrp,ints[j]);
                ll ji=prime[sqrp-i+1],F=4;
                while(ji<=ints[j]){
                    int dst=get(ints[j]/ji);
                    if(!fi0[dst]){
                        ans+=F*(1ll+4ll*(ll)how(sqrp-i+2,sqrp,ints[j]/ji));
                    }
                    else{
                        ans+=F*f[(i-1)&1][dst];
                    }
                    ji*=prime[sqrp-i+1];F+=3;
                }
                f[i&1][j]=ans;
            }
        }
    }
    void ywymain(){
        prim(400000);int t=get();
        while(t){
            t--;n=get();
            if(n<=400000){
                iprint(sums[n]);putchar('\n')；continue;
            }
            sqr=sqrt(n),sqrp=erfen(n);pos=0;
            for(register ll i=1;i<=n;){
                ints[pos]=n/i;pos++;i=n/(n/i)+1;
            }
            for(register int i=1;i<=sqr;){
                ints[pos]=sqr/i;pos++;i=sqr/(sqr/i)+1;
            }
            sort(ints,ints+pos);
            pos=unique(ints,ints+pos)-ints;
            for(register int i=0;i<pos;i++){
                insert(ints[i],i);
            }
            getg();getf();
            ll a=f[sqrp&1][get(n)],b=f[sqrp&1][get(sqr)];
            if(!fi0[get(n)])a=1ll+4ll*how(1,sqrp,n);
            if(!fi0[get(sqr)])b=1ll+4ll*how(1,sqrp,sqr);
            ll ans=a-b;
            for(register int i=1;i<=sqr;i++){
                int ywy=get(n/i);
                ll cjr=g[sqrp&1][ywy];
                if(i0[ywy])cjr=rg[ywy]-how(i0[ywy]+1,sqrp,n/i);
                ans+=((cjr-1ll)*4ll+1ll)*sig[i];
            }
            iprint(ans);putchar('\n');
            for(register int i=0;i<pos;i++){
                heads[ints[i]%p]=0;
            }gn=1;
        }
    }
}
signed main(){
    ywy::ywymain();return(0);
}