Wannafly Camp 2020 Day 3D 求和

杜教筛求 (phi(n))，

[S(n)=n(n+1)/2-sum_{d=2}^n S(frac{n}{d}) ]

答案为

[sum_{d=1}^n phi(d) h(frac{n}{d}) ]

其中 (h(n)=sum_{i=1}^n i^2)

顺便学习了一波 unordered_map

#include <bits/stdc++.h>
#include <unordered_map>
using namespace std;

#define int long long
const int N = 2e6+5;

int i2,i4,i6,mod,n;

bool isNotPrime[N + 5];
int mu[N + 5], phi[N + 5], primes[N + 5], cnt, sum[N+5];
unordered_map<int,int> mp;
inline void euler() {
    isNotPrime[0] = isNotPrime[1] = true;
    mu[1] = 1;
    phi[1] = 1;
    for (int i = 2; i <= N; i++) {
        if (!isNotPrime[i]) {
            primes[++cnt] = i;
            mu[i] = -1;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= cnt; j++) {
            int t = i * primes[j];
            if (t > N) break;
            isNotPrime[t] = true;
            if (i % primes[j] == 0) {
                mu[t] = 0;
                phi[t] = phi[i] * primes[j];
                break;
            } else {
                mu[t] = -mu[i];
                phi[t] = phi[i] * (primes[j] - 1);
            }
        }
    }
}

inline void exgcd(int a,int b,int &x,int &y) {
    if(!b) {
        x=1,y=0;
        return;
    }
    exgcd(b,a%b,x,y);
    int t=x;
    x=y,y=t-(a/b)*y;
}

inline int inv(int a,int b) {
    int x,y;
    return exgcd(a,b,x,y),(x%b+b)%b;
}

int S(int n) {
    if(n<=N) return sum[n];
    if(mp[n]) return mp[n];
    int ans=n*(n+1)%mod*i2%mod;
    int l=2,r;
    while(l<=n) {
        r=n/(n/l);
        ans-=(r-l+1)*S(n/l)%mod;
        ans%=mod;
        ans+=mod;
        ans%=mod;
        l=r+1;
    }
    mp[n]=ans;
    return ans;
}

int h(int n) {
    return (n*(n+1)%mod*(2*n+1)%mod*i6%mod);
}

signed main() {
    cin>>n>>mod;
    i2=inv(2,mod); i4=inv(4,mod); i6=inv(6,mod);
    euler();
    for(int i=1;i<=N;i++) sum[i]=sum[i-1]+phi[i], sum[i]%=mod;
    int l=1,r,ans=0;
    while(l<=n) {
        r=n/(n/l);
        ans+=(S(r)-S(l-1))*h(n/l)%mod;
        ans%=mod;
        ans+=mod;
        ans%=mod;
        l=r+1;
    }
    cout<<ans;
}