51nod 1847 奇怪的数学题

求解∑i=1n∑j=1nsgcd(i,j),sgcd\sum_{i=1}^n\sum_{j=1}^nsgcd(i,j),sgcd∑i=1n∑j=1nsgcd(i,j),sgcd表示次大公约数,n≤1010n\le{10^{10}}n≤1010

那么首先有sgcd(i,j)=gcd(i,j)/mn(gcd(i,j))sgcd(i,j)=gcd(i,j)/mn(gcd(i,j))sgcd(i,j)=gcd(i,j)/mn(gcd(i,j))就是gcd除以最小质因子。现在可以考虑交换枚举顺序先枚举gcd，然后最终得到的式子是这个形式∑d=1ndkmn(d)(∑i=1⌊nd⌋2φ(i)−1)\sum_{d=1}^n\frac{d^k}{mn(d)}(\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor }2\varphi(i)-1)∑d=1nmn(d)dk(∑i=1⌊dn⌋2φ(i)−1)

然后看到这个形式我们就可以整除分块了，后面这一部分可以杜教筛，关键在于前面，这个东西和最小质因子相关，我们可以利用min25筛的第一部分，因为我们每次转移就相当于是枚举了最小质因子，那么剩下的就是答案，然后初始化的时候需要利用第二类斯特林数求自然数幂之和，因为模数不是质数。

#include<bits/stdc++.h>
#define LL long long
#define uint unsigned int
#define V inline void
#define I inline uint
#define FOR(i,a,b) for(register int i=a,end##i=b;i<=end##i;++i)
#define REP(i,a,b) for(register int i=a,end##i=b;i>=end##i;--i)
#define go(i,x) for(int i=hed[x];i;i=e[i].pre)
using namespace std;
inline int read()
{char x='\0';int fh=1,sum=0;for(x=getchar();x<'0'||x>'9';x=getchar())if(x=='-')fh=-1;for(;x>='0'&&x<='9';x=getchar())sum=sum*10+x-'0';return fh*sum;
}
const int N=1e6+9;
uint n,k;
I ksm(uint a,uint b)
{uint sum=1;while(b){if(b&1)sum*=a;b>>=1;a*=a;}return sum;
}
uint s[59][59];
I getsum(uint x)
{uint ret=0;FOR(i,1,k){uint sum=1;FOR(j,0,i){if((x+1-j)%(i+1)==0)sum*=(x+1-j)/(i+1);else sum*=x+1-j;}ret+=s[k][i]*sum;}return ret;
}
uint phi[N],pri[N],pcnt;
bool vis[N];
uint sp0[N],spk[N];
V oula(uint mx)
{vis[0]=vis[1]=true;phi[1]=1;for(uint i=2;i<=mx;i++){if(!vis[i])pri[++pcnt]=i,phi[i]=i-1,sp0[pcnt]=sp0[pcnt-1]+1,spk[pcnt]=spk[pcnt-1]+ksm(i,k);for(uint j=1;j<=pcnt&&i*pri[j]<=mx;j++){vis[i*pri[j]]=true;phi[i*pri[j]]=phi[i]*(pri[j]-1);if(i%pri[j]==0){phi[i*pri[j]]=pri[j]*phi[i];break;}}}FOR(i,1,mx)phi[i]+=phi[i-1];
}
map<uint,uint>sphi;
I getphi(uint x)
{if(x<=N-9)return phi[x];if(sphi.find(x)!=sphi.end())return sphi[x];uint ans=x*(x+1)/2;for(uint lp=2,rp=0;lp<=x;lp=rp+1){uint now=x/lp;rp=x/now;ans-=(rp-lp+1)*getphi(now);}sphi[x]=ans;return ans;
}
uint sq;
uint w[N],tot,ind1[N],ind2[N];
uint g0[N],gk[N],ans[N];
int main()
{n=read(),k=read();sq=ceil(sqrt(n));oula(1e6);s[0][0]=1;FOR(i,1,k)FOR(j,1,i)s[i][j]=j*s[i-1][j]+s[i-1][j-1];for(uint lp=1,rp=0;lp<=n;lp=rp+1){uint now=n/lp;rp=n/now;w[++tot]=now;g0[tot]=now-1,gk[tot]=getsum(now)-1;
//      cout<<now<<' '<<g0[tot]<<' '<<gk[tot]<<endl;//if(now<=sq)ind1[now]=tot;else ind2[n/now]=tot;}for(uint i=1;i<=pcnt&&pri[i]<=sq;i++){for(uint j=1;j<=tot&&1LL*pri[i]*pri[i]<=w[j];j++){uint now=w[j]/pri[i];uint h=(now<=sq)?ind1[now]:ind2[n/now];gk[j]-=ksm(pri[i],k)*(gk[h]-spk[i-1]);ans[j]+=gk[h]-spk[i-1];g0[j]-=g0[h]-sp0[i-1];
//          cout<<pri[i]<<" "<<w[j]<<' '<<g0[j]<<' '<<gk[j]<<' '<<ans[j]<<endl;}}uint sum=0;for(uint lp=1,rp=0;lp<=n;lp=rp+1){uint now=n/lp;rp=n/now;uint x=(rp<=sq)?ind1[rp]:ind2[n/rp];uint y=(lp-1<=sq)?ind1[lp-1]:ind2[n/(lp-1)];
//      cout<<"now "<<now<<' '<<getphi(now)<<endl;
//      cout<<"ans "<<ans[x]<<' '<<ans[y]<<endl;//
//      cout<<"g "<<g0[x]<<' '<<g0[y]<<endl;//sum+=(2*getphi(now)-1)*(ans[x]+g0[x]-ans[y]-g0[y]);
//      cout<<"sum "<<sum<<endl;//}printf("%u",sum);return 0;
}

51nod 1847 奇怪的数学题（数论/min25筛/杜教筛/斯特林数）相关推荐

【bzoj4176】Lucas的数论莫比乌斯反演+杜教筛
题目描述去年的Lucas非常喜欢数论题,但是一年以后的Lucas却不那么喜欢了. 在整理以前的试题时,发现了这样一道题目"求Sigma(f(i)),其中1<=i<=N" ...
P3768-简单的数学题【莫比乌斯反演,杜教筛】
正题题目链接:https://www.luogu.com.cn/problem/P3768 题目大意给出n,pn,pn,p求∑i=1n∑j=1ngcd(i,j)∗i∗j\sum_{i=1}^n\s ...
【数论】【杜教筛】选数（P3172）
正题 P3172 题目大意在 [L,R] 选n个数,问gcd=k的方案数解题思路因为gcd=k,那么所选的数都是k的倍数,那么可以让L,R整除k,那么有 ∑a1=LR∑a2=LR...∑an=L ...
[LOJ]#572. 「LibreOJ Round #11」Misaka Network 与求和 min_25筛+杜教筛
Solution 推一下式子,容易得到一个线性做法:∑d=1nfk(d)((2∑i=1⌊ni⌋φ(i))−1)\sum_{d=1}^nf^k(d)((2\sum_{i=1}^{\lfloor{n\ov ...
【LOJ #572】【LibreOJ Round #11】—Misaka Network 与求和（min_25筛+杜教筛）
传送门首先可以随便莫反一波,可以得到 ans=∑d=1nf(d)k∑i=1nd∑j=1ndgcd(i,j)=1ans=\sum_{d=1}^{n}f(d)^k\sum_{i=1}^{\frac n ...
【Luogu3768】简单的数学题（莫比乌斯反演/杜教筛/欧拉函数）
[Luogu3768]简单的数学题 https://www.cnblogs.com/cjyyb/p/8298339.html
P3768 简单的数学题（杜教筛）
P3768 简单的数学题推式子 ∑i=1n∑j=1mijgcd(i,j)=∑d=1nd∑i=1n∑j=1mij(gcd(i,j)=d)=∑d=1nd3∑i=1nd∑j=1ndij∑k∣gcd(i,j ...
matlab狄利克雷函数,数论入门1——莫比乌斯函数,欧拉函数,狄利克雷卷积,线性筛,莫比乌斯反演,杜教筛...
数论入门1 一个菜鸡对数论的一点点理解... 莫比乌斯函数定义函数$\mu(n)$为: 当n有平方因子时,$\mu(n)=0$. 当n没有平方因子时,$\mu(n)=(-1)^{\omega(n)} ...
杜教筛及其时间复杂度分析
文章目录杜教筛方法举例莫比乌斯函数欧拉函数时间复杂度杜教筛杜教筛用于求一类积性函数的前缀和,时间复杂度可以做到 O(n23)O(n^{\frac{2}{3}})O(n32). 方法 ...

51nod 1847 奇怪的数学题（数论/min25筛/杜教筛/斯特林数）

51nod 1847 奇怪的数学题

51nod 1847 奇怪的数学题（数论/min25筛/杜教筛/斯特林数）相关推荐

最新文章

热门文章