[51nod1847]奇怪的数学题

description

51nod
求$$\sum_{i=1}^{{n}\sum_{j=1}}sgcd(i,j)^k$$其中\(sgcd(i,j)\)表示\(i,j\)的次大公约数,如果\(gcd(i,j)=1\)那么\(sgcd(i,j)=0\)。

solution

记答案为\(Ans\)。
首先考虑直接枚举\(sgcd(i,j)\)。

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

其中当\(n\not=1\)时\(\xi^k(n)=\frac{n}{p_{min}(n)}\)。
这个时候如果你像菜鸡fdf一样瞎反演出
\(\sum_{i=1}^{n}\sum_{j=1}^n[gcd(i,j)==d]=\sum_{t|d}\lfloor\frac{n}{t}\rfloor^2\mu(\frac{t}{d})\)就会变成这样

\[\begin{aligned} Ans=&\sum_{d=1}^{n}\xi^k(d)\sum_{t|d}\lfloor\frac{n}{t}\rfloor^2\mu(\frac{t}{d}) \\ =&\sum_{t=1}^{n}\lfloor\frac{n}{t}\rfloor^2\sum_{t|d}\xi^k(d)\mu(\frac{t}{d}) \\ \end{aligned} \]

假设我们对于前面数论分块,那么现在要求

\[\begin{aligned} \sum_{i=1}^{n}\sum_{n|d}\xi^k(d)\mu(\frac{d}{n})=&\sum_{i=1}^{n}\xi^k(i)\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\mu(j)\\ \end{aligned} \]

你会发现这玩意又要数论分块一次。
然后菜鸡fdf就华丽地得到了一个\(O(n)\)的做法。O(n)可能也可以做

我们注意\(\sum_{i=1}^{n}\sum_{j=1}^n[gcd(i,j)==d]\)中\(i,j\)的上界都是\(n\)。
转化一下变成求

\[\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[gcd(i,j)==1]=2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\phi(i)-1 \]

这玩意不就是枚举其中一个数+去重么。
这应该是一个经典等式,然而我因为太菜没有做过。
于是原式变为

\[Ans=\sum_{d=1}^{n}\xi^k(d)(2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\phi(i)-1) \]

这样只要一次数论分块即可。

Code

我不管我不管我不会杜教筛只会暴力Min_25

#include<bits/stdc++.h>
#define FL "a"
using namespace std;
typedef unsigned int ll;
typedef long double dd;
const int N=1e5+10;
const int mod=1e9+7;
inline ll read(){
  ll data=0,w=1;char ch=getchar();
  while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
  if(ch=='-')w=-1,ch=getchar();
  while(ch<='9'&&ch>='0')data=data*10+ch-48,ch=getchar();
  return data*w;
}
inline void file(){
  freopen(FL".in","r",stdin);
  freopen(FL".out","w",stdout);
}

inline ll poww(ll a,ll b){
  ll res=1;
  for(;b;b>>=1,a*=a)
    if(b&1)res*=a;
  return res;
}

int cnt;ll n,k,pri[N],pw[N][52],sum[N],sumk[N],s[52][52];bool vis[N];
inline void sieve(){
  register int i,j;
  for(vis[1]=1,i=2;i<N;i++){
    if(!vis[i])pri[++cnt]=i;
    for(j=1;j<=cnt&&i*pri[j]<N;j++){
      vis[i*pri[j]]=1;if(i%pri[j]==0)break;
    }
  }
  for(i=1;i<=cnt;i++)
    for(j=pw[i][0]=1;j<=50;j++)pw[i][j]=pw[i][j-1]*pri[i];
  for(i=1;i<=cnt;i++){
    sum[i]=sum[i-1]+pri[i];
    sumk[i]=sumk[i-1]+poww(pri[i],k);
  }
  s[0][0]=1;
  for(i=0;i<=50;i++)
    for(j=1;j<=i;j++)
      s[i][j]=s[i-1][j-1]+j*s[i-1][j];
}
inline ll getsum(ll sn,ll sk){
  register ll res=0,tmp,i,j;
  for(i=0;i<=sk;i++){
    tmp=1;
    for(j=0;j<=i;j++)
      if((sn-j+1)%(i+1)==0)tmp*=(sn-j+1)/(i+1);
      else tmp*=(sn-j+1);
    res+=s[sk][i]*tmp;
  }
  return res;
}

ll w[N],g[N],xi[N],ps[N],phi[N];int sqr,tot,id1[N],id2[N];
#define ID(x) (x<=sqr?id1[x]:id2[n/(x)])
inline void Min25(){
  register ll i,j,a,b,e,ans;
  for(sqr=sqrt(n),tot=0,i=1;i<=n;i=j+1){
    j=n/(n/i);w[++tot]=n/i;
    w[tot]<=sqr?id1[w[tot]]=tot:id2[n/w[tot]]=tot;
  }
  
  for(i=1;i<=tot;i++)g[i]=getsum(w[i],k)-1;
  for(i=1;i<=cnt&&1ll*pri[i]*pri[i]<=n;i++)
    for(j=1;1ll*pri[i]*pri[i]<=w[j];j++){
      g[j]-=pw[i][k]*(g[ID(w[j]/pri[i])]-sumk[i-1]);
      xi[j]+=g[ID(w[j]/pri[i])]-sumk[i-1];
    }
  
  for(i=1;i<=tot;i++)g[i]=w[i]-1;
  for(i=1;i<=cnt&&1ll*pri[i]*pri[i]<=n;i++)
    for(j=1;1ll*pri[i]*pri[i]<=w[j];j++)
      g[j]-=g[ID(w[j]/pri[i])]-i+1;
  for(i=1;i<=tot;i++)ps[i]-=g[i],xi[i]+=g[i];
  
  for(i=1;i<=tot;i++)g[i]=1ll*w[i]*(w[i]+1)/2-1;
  for(i=1;i<=cnt&&1ll*pri[i]*pri[i]<=n;i++)
    for(j=1;1ll*pri[i]*pri[i]<=w[j];j++)
      g[j]-=pri[i]*(g[ID(w[j]/pri[i])]-sum[i-1]);
  for(i=1;i<=tot;i++)ps[i]+=g[i];

  for(i=cnt;i;i--)
    if(1ll*pri[i]*pri[i]<=n)
      for(j=1;1ll*pri[i]*pri[i]<=w[j];j++)
	for(e=1,a=pri[i];1ll*a*pri[i]<=w[j];e++,a*=pri[i])
	  phi[j]+=(pw[i][e]-pw[i][e-1])*(phi[ID(w[j]/a)]+ps[ID(w[j]/a)]-(sum[i]-i))+(pw[i][e+1]-pw[i][e]);
  for(i=1;i<=tot;i++)phi[i]+=ps[i]+1;  
  ans=0;
  for(i=1;i<=n;i=j+1)
    j=n/(n/i),ans+=(xi[ID(j)]-xi[ID(i-1)])*(2*phi[ID(n/i)]-1);
  printf("%u\n",ans);
}

int main()
{
  n=read();k=read();sieve();Min25();
  return 0;
}