【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 <script type="math/tex" id="MathJax-Element-45">I</script>，都是積性函數，那麼右邊也是積性函數

然後自己推一推就可以發現了

這樣一來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=,mo=<<;
ll h[N+],h1[N+],n,ans;
int k,phi[N+],p[N+],f[N+],lemuda[N+],sum[N+],mx[N+],cnt[N+],mu[N+];
bool bz[N+],vis[N+],vis1[N+];
ll mi(ll x,int y) {
    ll z=;
    for(;y;y/=,x=x*x)
        if (y&) z=z*x;
    return z;
}
void prepare() {
    phi[]=;lemuda[]=;mu[]=;
    fo(i,,N) {
        if (!bz[i]) p[++p[]]=i,lemuda[i]=mu[i]=-,phi[i]=i-,mx[i]=cnt[i]=,f[i]=i;
        fo(j,,p[]) {
            int k=i*p[j];
            if (k>N) break;
            bz[k]=;lemuda[k]=-lemuda[i];
            if (p[j]==f[i]) cnt[k]=cnt[i]+;
            else cnt[k]=;
            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,,N) {
        sum[i]=lemuda[i]+mo;
        (sum[i]+=sum[i-])&=(mo-);
        phi[i]+=phi[i-];
        phi[i]=phi[i]&(mo-);
    }
    f[]=k;
    fo(i,,N) {
        if (mx[i]<=k) f[i]=lemuda[i]*(k-mx[i]+);
        else f[i]=;
        (f[i]+=f[i-]+mo)&=(mo-);
    }
}
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&) ans=(x%mo)*((x+)/%mo);
    else ans=(x/%mo)*((x+)%mo);
    ans=ans&(mo-);
    for(ll l=,r;l<=x;l=r+) {
        r=x/(x/l);
        ans+=mo-(r-l+)%mo*get_phi(x/l)%mo;
        ans=ans&(mo-);
    }
    vis[t]=;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,/);
    for(ll l=,r;l<=x;l=r+) {
        r=x/(x/l);
        ans+=mo-(r-l+)%mo*get_lemuda(x/l)%mo;
        ans=ans&(mo-);
    }
    vis1[t]=;h1[t]=ans;
    return ans;
}
ll solve(int d,ll n,ll y) {
    ll m=pow(n,/(d+));
    ll ans=;
    fo(i,,m) { 
        int z=mu[i];
        if (!z) continue;
        if (lemuda[i]==-&&!(d&)) z=-z;
        ll x=mi(i,d+);
        if (z==) {
            if (n/x<=N) (ans+=sum[n/x])&=(mo-);
            else (ans+=h1[x*y])&=(mo-);
        } else {
            if (n/x<=N) (ans+=mo-sum[n/x])&=(mo-);
            else (ans+=mo-h1[x*y])&=(mo-);
        }
    }
    return ans;
}
int main() {
    scanf("%lld%d",&n,&k);
    prepare();get_phi(n);get_lemuda(n);ll la=;
    for(ll l=,r;l<=n;l=r+) {
        r=n/(n/l);ll sum=,now=;
        if (n/l<=N) sum=*phi[n/l]-;
        else sum=*h[l]-;
        sum&=(mo-);
        if (r<=N) now=f[r];
        else fo(j,,k) (now+=solve(j,r,n/l))&=(mo-);
        ll tmp=now-la;
        if (tmp<) tmp+=mo;
        tmp=tmp&(mo-);
        tmp=tmp*sum&(mo-);
        (ans+=tmp)&=(mo-);
        la=now;
    }
    printf("%lld\n",ans);
}