EOJ Monthly 2019.11 E-数学题 莫比乌斯反演+杜教筛+拉格朗日插值

传送门： https://acm.ecnu.edu.cn/contest/231/problem/E/#report16

题意

$$求解\sum _{i=1}^n\sum _{a_1=1}^i...\sum _{a_k=1}^i[gcd(a_1,a_2..a_k,i)=1]$$

思路

$$\sum _{i=1}^n\sum _{a_1=1}^i...\sum _{a_k=1}^i[gcd(a_1,a_2..a_k,i)=1]$$

$$\sum _{i=1}^n\sum _{a_1=1}^i...\sum _{a_k=1}^i\sum _{d|a_1..d|a_{k}d|i}\mu (i)$$

$$\sum _{i=1}^n\sum _{d|i}\mu (d)(\frac{i}{d}{)}^k$$

$$\sum _{d=1}^n\mu (d)\sum _{i=1}^{\frac{n}{d}}{i}^k$$

$因为n有1e9，所以杜教筛处理，{\sum }_{i=1}^n{i}^k用拉格朗日插值/伯努利数处理即可。$

Code

#include "bits/stdc++.h"
using namespace std;

typedef long long ll;

const ll mod = 998244353;

const int N = 1e6 + 10;

namespace polysum {
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
    const int D=1010000;///可能需要用到的最高次
    ll a[D],f[D],g[D],p[D],p1[D],p2[D],b[D],h[D][2],C[D], num[D];
    ll powmod(ll a,ll b){ll res=1;a%=mod;assert(b>=0);for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}

    ///函数用途：给出数列的（d+1）项，其中d为最高次方项
    ///求出数列的第n项，数组下标从0开始
    ll calcn(int d,ll *a,ll n) { /// a[0].. a[d]  a[n]
        if (n<=d) return a[n];
        p1[0]=p2[0]=1;
        rep(i,0,d+1) {
            ll t=(n-i+mod)%mod;
            p1[i+1]=p1[i]*t%mod;
        }
        rep(i,0,d+1) {
            ll t=(n-d+i+mod)%mod;
            p2[i+1]=p2[i]*t%mod;
        }
        ll ans=0;
        rep(i,0,d+1) {
            ll t=g[i]*g[d-i]%mod*p1[i]%mod*p2[d-i]%mod*a[i]%mod;
            if ((d-i)&1) ans=(ans-t+mod)%mod;
            else ans=(ans+t)%mod;
        }
        return ans;
    }
    void init(int M) {///用到的最高次
        f[0]=f[1]=g[0]=g[1]=1;
        rep(i,2,M+5) f[i]=f[i-1]*i%mod;
        g[M+4]=powmod(f[M+4],mod-2);
        per(i,1,M+4) g[i]=g[i+1]*(i+1)%mod;///费马小定理筛逆元
    }

    ///函数用途：给出数列的（m+1）项，其中m为最高次方
    ///求出数列的前（n-1）项的和（从第0项开始）
    ll polysum(ll m,ll *a,ll n) { /// a[0].. a[m] \sum_{i=0}^{n-1} a[i]
        for(int i=0;i<=m;i++) b[i]=a[i];

        ///前n项和，其最高次幂加1
        b[m+1]=calcn(m,b,m+1);
        rep(i,1,m+2) b[i]=(b[i-1]+b[i])%mod;
        return calcn(m+1,b,n-1);
    }
 
    ll solve(ll n, int k) {
        ll ans = polysum(k + 1, num, n) % mod;
        return ans;
    }
}


bool is_prime[N];
int prime[N], cnt, mu[N];
ll n, k;

void init() {
    mu[1] = 1;
    for(int i = 2;i < N; i++) {
        if(!is_prime[i]) prime[++cnt] = i, mu[i] = -1;
        for(int j = 1;j <= cnt && i * prime[j] < N; j++) {
            is_prime[i * prime[j]] = 1;
            if(i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            }
            else mu[i * prime[j]] = -mu[i];
        }
    }
    for(int i = 2;i < N; i++) mu[i] += mu[i - 1];

    polysum::init(k + 10);
    for(int i = 0;i <= k + 1; i++) polysum::num[i] = polysum::powmod((ll)i + 1, k);
}

map<ll, ll> mp;
ll S(ll x) {
    if(x < N) return mu[x];
    if(mp[x]) return mp[x];
    ll ans = 1;
    for(int l = 2, r;l <= x; l = r + 1) {
        r = min(x, x / (x / l));
        ans = ans - (r - l + 1) * S(x / l);
    }
    return mp[x] = ans;
}

void solve() {
    cin >> n >> k;
    init();
    ll ans = 0;
    for(ll l = 1, r;l <= n; l = r + 1) {
        r = min(n, n / (n / l));
        ans = (ans + (S(r) - S(l - 1)) * polysum::solve(n / l, k) % mod) % mod;
    }
    cout << (ans % mod + mod) % mod << endl;
}

signed main() {
    solve();
}