【GDOI2018模拟8.12】求和

Description

求∑i=1n∑j=1n∑d=1kfd(gcd(i,j))

其中当n=∏ipaii时fd(n)=∏i(−1)ai[ai<=d]

答案对2^30取模

n<=10^10,k<=40

Solution

首先是个懂数论的人都能把式子化成这样：

Ans=∑i=1n∑d=1kfd(i)∑j=1⌊ni⌋φ(j)

后面这东西直接用杜教筛做

关键我们要怎么求前面这东西呢？

考虑当n较小的时候，我们可以预处理出g(n)=∑ni=1∑kd=1fd(i)

你只需要线筛出每个数的最大指数mx

如果mx>k那么g(i)显然是0

否则我们引入λ函数，lambda(n)=∏i(−1)ai

那么g(i)=λk−mx+1(i)

但是当n较大的时候呢？

我们设Fd=∑ni=1fd(i)

我们可以发现Fd(n)=∑d+1(√n)i=1μ(i)∑⌊nid+1⌋j=1λ(id+1j)

这个根据容斥原理退一下就好了

然后我们设W(n)=∑ni=1λ(i)（实在不知道大写λ怎么打）

那么根据λ是完全积性函数，式子可以写成Fd(n)=∑d+1(√n)i=1λd+1(i)μ(i)W(⌊nid+1⌋)

如果我们求出了W，那么这个东西我们就可以直接暴力求了。

现在问题是W怎么求？

我们可以发现∑d|nλ(d)=[n是完全平方数]

因为狄利克雷卷积左边是λ和I，都是积性函数，那么右边也是积性函数

然后自己推一推就可以发现了

这样一来W我们就可以通过杜教筛来解决，总复杂度O(n^2/3)

如果常数写打了被卡，发现模数是2的幂可以直接位运算加速

其实也可以unsigned long long自然溢出_ (:з」∠) _

Code

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
typedef long long ll;
const int N=1e7,mo=1<<30;
ll h[N+5],h1[N+5],n,ans;
int k,phi[N+5],p[N+5],f[N+5],lemuda[N+5],sum[N+5],mx[N+5],cnt[N+5],mu[N+5];
bool bz[N+5],vis[N+5],vis1[N+5];
ll mi(ll x,int y) {
ll z=1;
for(;y;y/=2,x=x*x)
if (y&1) z=z*x;
return z;
}
void prepare() {
phi[1]=1;lemuda[1]=1;mu[1]=1;
fo(i,2,N) {
if (!bz[i]) p[++p[0]]=i,lemuda[i]=mu[i]=-1,phi[i]=i-1,mx[i]=cnt[i]=1,f[i]=i;
fo(j,1,p[0]) {
int k=i*p[j];
if (k>N) break;
bz[k]=1;lemuda[k]=-lemuda[i];
if (p[j]==f[i]) cnt[k]=cnt[i]+1;
else cnt[k]=1;
mx[k]=max(mx[i],cnt[k]);
if (!(i%p[j])) {
phi[k]=phi[i]*p[j];
f[k]=p[j];
break;
}
phi[k]=phi[i]*phi[p[j]];
if (!f[k]) f[k]=p[j];
mu[k]=-mu[i];
}
}
fo(i,1,N) {
sum[i]=lemuda[i]+mo;
(sum[i]+=sum[i-1])&=(mo-1);
phi[i]+=phi[i-1];
phi[i]=phi[i]&(mo-1);
}
f[1]=k;
fo(i,2,N) {
if (mx[i]<=k) f[i]=lemuda[i]*(k-mx[i]+1);
else f[i]=0;
(f[i]+=f[i-1]+mo)&=(mo-1);
}
}
ll get_phi(ll x) {
if (x<=N) return phi[x];int t=n/x;
if (vis[t]) return h[t];ll ans;
if (x&1) ans=(x%mo)*((x+1)/2%mo);
else ans=(x/2%mo)*((x+1)%mo);
ans=ans&(mo-1);
for(ll l=2,r;l<=x;l=r+1) {
r=x/(x/l);
ans+=mo-(r-l+1)%mo*get_phi(x/l)%mo;
ans=ans&(mo-1);
}
vis[t]=1;h[t]=ans;
return ans;
}
ll get_lemuda(ll x) {
if (x<=N) return sum[x];int t=n/x;
if (vis1[t]) return h1[t];
ll ans=pow(x,1.0/2);
for(ll l=2,r;l<=x;l=r+1) {
r=x/(x/l);
ans+=mo-(r-l+1)%mo*get_lemuda(x/l)%mo;
ans=ans&(mo-1);
}
vis1[t]=1;h1[t]=ans;
return ans;
}
ll solve(int d,ll n,ll y) {
ll m=pow(n,1.0/(d+1));
ll ans=0;
fo(i,1,m) {
int z=mu[i];
if (!z) continue;
if (lemuda[i]==-1&&!(d&1)) z=-z;
ll x=mi(i,d+1);
if (z==1) {
if (n/x<=N) (ans+=sum[n/x])&=(mo-1);
else (ans+=h1[x*y])&=(mo-1);
} else {
if (n/x<=N) (ans+=mo-sum[n/x])&=(mo-1);
else (ans+=mo-h1[x*y])&=(mo-1);
}
}
return ans;
}
int main() {
scanf("%lld%d",&n,&k);
prepare();get_phi(n);get_lemuda(n);ll la=0;
for(ll l=1,r;l<=n;l=r+1) {
r=n/(n/l);ll sum=0,now=0;
if (n/l<=N) sum=2*phi[n/l]-1;
else sum=2*h[l]-1;
sum&=(mo-1);
if (r<=N) now=f[r];
else fo(j,1,k) (now+=solve(j,r,n/l))&=(mo-1);
ll tmp=now-la;
if (tmp<0) tmp+=mo;
tmp=tmp&(mo-1);
tmp=tmp*sum&(mo-1);
(ans+=tmp)&=(mo-1);
la=now;
}
printf("%lld\n",ans);
}