【51nod1847】奇怪的数学题（Min_25筛+杜教筛）

最新推荐文章于 2021-08-21 13:11:40 发布

转载最新推荐文章于 2021-08-21 13:11:40 发布 · 184 阅读

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/bztMinamoto/p/10420143.html

本文深入探讨了一种解决特定数学问题的方法，利用质因数分解和杜教筛算法来高效计算一系列数学表达式的和。通过引入次大质因子的概念，文章详细解析了如何将复杂求和问题转化为更易于处理的形式，并使用杜教筛算法进行优化。此外，还讨论了如何在非质数模意义下计算多项式和，采用第二类斯特林数进行优化。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题面

传送门

题解

这题有毒……不知为啥的错误调了半天……

令\(f(i)={sgcd(i)}\)，那么容易看出\(f(i)\)就是\(i\)的次大质因子，用\(i\)除以它的最小质因子即可计算

于是题目所求即为

\[\sum_{i=1}^n\sum_{j=1}^n{f(\gcd(i,j))}^k\]

\[\sum_{d=1}^n {f(d)}^k\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,j)=1]\]

\[\sum_{d=1}^n {f(d)}^k(2\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\varphi(i)-1)\]

后面那个根据\(\varphi\)的定义就可以推出来

因为后面括号中的式子只与\({\left\lfloor\frac{n}{d}\right\rfloor}\)的值有关，我们可以考虑整除分块，那么括号里的东西就是一个杜教筛

那么我们现在的问题就是对于前半部分怎么快速求和了。我们考虑\(Min\_25\)筛的过程，其中\(g(n,j)\)表示所有\(1\)到\(n\)中的质数或最小质因子大于\(P_j\)的所有数的\(k\)次方之和，即

\[g(n,j)=\sum_{i=1}^ni^k[i\in P\ or\ min_p(i)>P_j]\]

中间的转移中有这么一步

\[g(n,i)=g(n,j-1)-{P_j}^k(g({\left\lfloor\frac{n}{P_j}\right\rfloor},j-1)-\sum_{i=1}^{j-1}{P_i}^k)\]

考虑后面减去的东西，是所有最小质因子等于\(P_j\)且自身为合数的数的\(k\)次方之和，即

\[\sum_xx^k[x\notin P\ and\ min_p(x)=P_j]\]

然后我们惊喜地发现，后面减去的东西中，括号里的东西就是上面所有满足条件的\(x\)的\({f(x)}^k\)之和！

那么我们就可以直接\(Min\_25\)筛来计算\({f(x)}^k\)的前缀和了，那么一段区间只要差分一下就可以了

最后有个尴尬的问题，就是这个模数不是质数，所以我们在初始化的时候需要快速计算\(\sum_{i=1}^n i^k\)就不能拉格朗日插值了

那么只好用第二类斯特林数优化了

\[\begin{aligned} \sum\limits_{i=1}^{n}i^K&=\sum_{i=1}^n\sum_{j=1}^K\begin{Bmatrix}K\\j\end{Bmatrix}i^\underline{j}\\ &=\sum_{j=1}^K\begin{Bmatrix}K\\j\end{Bmatrix}\sum_{i=1}^ni^\underline{j}\\ &=\sum_{j=1}^K\begin{Bmatrix}K\\j\end{Bmatrix}\frac{{(n+1)}^\underline{j+1}}{j+1} \end{aligned}\]

然后没有然后了

//minamoto
#include<bits/stdc++.h>
#define R register
#define int unsigned int
#define IT map<int,int>::iterator
#define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
const int N=1e5+5;
int w[N<<1],g[N<<1],h[N<<1],G[N<<1],id1[N],id2[N],p[N],pw[N],sp[N],phi[N],S[55][55];
map<int,int>mp;bitset<N>vis;int n,m,k,tot,sqr,id,ans,now,las;
int ksm(R int x,R int y){
    R int res=1;
    for(;y;y>>=1,x*=x)if(y&1)res*=x;
    return res;
}
void init(int n){
    phi[1]=1;
    fp(i,2,n){
        if(!vis[i])p[++tot]=i,pw[tot]=ksm(i,k),sp[tot]=sp[tot-1]+pw[tot],phi[i]=i-1;
        for(R int j=1;j<=tot&&1ll*i*p[j]<=n;++j){
            vis[i*p[j]]=1;
            if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}
            phi[i*p[j]]=phi[i]*(p[j]-1);
        }
    }
    fp(i,1,n)phi[i]+=phi[i-1];
    S[0][0]=1;
    fp(i,1,k){
        S[i][1]=1;
        fp(j,2,i)S[i][j]=S[i-1][j]*j+S[i-1][j-1];
    }
}
int Phi(int n){
    if(n<=sqr)return phi[n];
    IT it=mp.find(n);
    if(it!=mp.end())return it->second;
    int res=n*(n+1)>>1;
    for(R int i=2,j;i<=n;i=j+1)
        j=n/(n/i),res-=Phi(n/i)*(j-i+1);
    return mp[n]=res;
}
int calc(int n){
    int res=0,tmp;
    fp(i,1,min(k,n)){
        tmp=S[k][i];
        fp(j,n-i+1,n+1)(j%(i+1)==0)?tmp*=j/(i+1):tmp*=j;
        res+=tmp;
    }return res;
}
signed main(){
//  freopen("testdata.in","r",stdin);
    scanf("%u%u",&n,&k),init(sqr=sqrt(n));
    for(R int i=1,j;i<=n;i=j+1){
        j=n/(n/i),w[++m]=n/i;
        w[m]<=sqr?id1[w[m]]=m:id2[n/w[m]]=m;
        g[m]=calc(w[m])-1,h[m]=w[m]-1;
    }
    fp(j,1,tot)for(R int i=1;1ll*p[j]*p[j]<=w[i];++i){
        id=(w[i]/p[j]<=sqr)?id1[w[i]/p[j]]:id2[n/(w[i]/p[j])];
        g[i]-=pw[j]*(g[id]-sp[j-1]);
        h[i]-=h[id]-j+1;
        G[i]+=g[id]-sp[j-1];
    }
    for(R int i=2,j;i<=n;i=j+1){
        j=n/(n/i),id=(j<=sqr)?id1[j]:id2[n/j];
        now=G[id]+h[id],ans+=((Phi(n/i)<<1)-1)*(now-las),las=now;
    }
    printf("%u\n",ans);
    return 0;
}

转载于:https://www.cnblogs.com/bztMinamoto/p/10420143.html