【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);}