【GDOI2018Day1模拟4.17】呼吸决定

Description：

1<=n<=10^9,1<=m<=200000

空间30720KB，时间1s,开O2

题解：

就是杜教筛+拉格朗日插值的结合题。

考场时没有打完，可惜了。

顺便回忆一下这两种算法。

设$s(n)=\sum_{i=1}^n i^m*\mu(i)$

$\sum_{i=1}^n i^m\sum_{j|i} \mu(j)$
$=\sum_{i=1}^n i^m\sum_{j=1}^{{\lfloor {n\over i} \rfloor}}\mu(j)*j^m$
$=\sum_{i=1}^n i^m s({\lfloor {n\over i} \rfloor})$
提出i=1的情况，移项，得：
$s(1)=\sum_{i=2}^n i^m s({\lfloor {n\over i} \rfloor})$

预处理$n^{2/3}$内的结果，之后的分块+hash。

现在问题在于快速求$\sum_{i=1}^n i^m$。

差分法可得这是一个m+1次多项式，于是我们可以通过插m+2个值来求解。

一般插0-m+1会好一点，注意一定是0-m+1。

拉格朗日插值：

$f(X)=\sum_{i=0}^{n} y_i\prod_{j≠i}{X - x_j\over x_i - x_j}$

显然这个东西展开后是n次的，并且把$x_i$代进去是对的。

对于拉格朗日插值求自然数幂和，要处理一下${X - x_j}$的前后缀积，和阶乘的逆元。

这样就可以$O(m)$算一次了。

这个预处理500000个，剩下200个，200*200000=400000，追求空间与时间的均衡。

Code：

#include<cstdio>
#include<bitset>
#define ll long long
#define fo(i, x, y) for(int i = x; i <= y; i ++)
#define fd(i, x, y) for(int i = x; i >= y; i --)
using namespace std;
const int mo = 998244353;
const int N = 1e6, C = 2e5 + 1, E = 5e6 + 5;
int n, m;
ll ksm(ll x, ll y) {
    ll s = 1;
    for(; y; y /= 2, x = x * x % mo)
        if(y & 1) s = s * x % mo;
    return s;
}
bitset<N + 1> bz;
int p0[350000], mu[N + 1];
int nf[C + 1], p[C + 1], q[C + 1], cm[E + 1];
void BB() {
    mu[1] = 1; cm[1] = 1;
    fo(i, 2, N)  {
        if(!bz[i]) p0[++ p0[0]] = i, mu[i] = -ksm(i, m);
        fo(j, 1, p0[0]) {
            int k = i * p0[j]; if(k > N) break;
            bz[k] = 1;
            if(i % p0[j] == 0) {mu[k] = 0; break;}
            mu[k] = (ll) mu[i] * mu[p0[j]] % mo;
        }
    }
    fo(i, 1, N) mu[i] = (mu[i - 1] + mu[i]) % mo;
    ll fac = 1;
    fo(i, 1, C) fac = fac * i % mo, nf[i] = ksm(fac, mo - 2);
    nf[0] = 1;
    fo(i, 1, E) cm[i] = (cm[i - 1] + ksm(i, m)) % mo;
    m ++;
}
const int M = 131321;
int h2[M][2];
int ha2(int x) {
    int y = x % M;
    while(h2[y][0] && h2[y][0] != x) y = (y + 1) % M;
    return y;
}
ll sum(ll n) {
    if(n <= E) return cm[n];
    int y = ha2(n);
    if(h2[y][0] == n) return h2[y][1];
    h2[y][0] = n;
    ll ans = 0;
    p[0] = n; fo(i, 1, m) p[i] = (ll) p[i - 1] * (n - i) % mo;
    q[m + 1] = 1; fd(i, m, 0) q[i] = (ll) q[i + 1] * (n - i) % mo;
    fo(i, 1, m) ans += (ll) cm[i] * p[i - 1] % mo * q[i + 1] % mo * nf[i] % mo * nf[m - i] * ((m - i) % 2 ? -1 : 1) % mo;
    ans = (ans % mo + mo) % mo;
    h2[y][1] = ans;
    return ans;
}
int h[M][2];
int ha(int x) {
    int y = x % M;
    while(h[y][0] && h[y][0] != x) y = (y + 1) % M;
    return y;
}
ll dg(int n) {
    if(n <= N) return (mu[n] + mo) % mo;
    int y = ha(n);
    if(h[y][0] == n) return h[y][1];
    h[y][0] = n;
    ll ans = 1;
    fo(i, 2, n) {
        int j = n / (n / i);
        ans -= dg(n / i) * (sum(j) - sum(i - 1) + mo) % mo;
        i = j;
    }
    ans = (ans % mo + mo) % mo; h[y][1] = ans;
    return ans;
}
int main() {
    freopen("calc.in", "r", stdin);
    freopen("calc.out", "w", stdout);
    scanf("%d %d", &n, &m);
    BB();
    printf("%lld\n", dg(n));
}