min25筛

这也是个用于求积性函数前缀和的筛法。

对于积性函数 f(n),满足这样几个要求时可以用min25来做。

  • n \in prime 时, f(n) 可以使用关于 n 的多项式表示。

  • n \in prime 时,f(n^k) 要能够快速计算答案。

预处理 Σf(p)

在min25筛的过程中,我们需要预处理多项式的 x^k 相关的这个东西

\sum_{i=2}^n i^k [i \in prime]

那么记所有质数的集合为 P,第 j 个素数为 P_jx 的最小质因子为 minp(x),我们需要预处理这样一个东西

g(n,j) = \sum_{i=2}^n i^k [i \in P ~ or ~ minp(i) > P_j]

举个例子吧

\begin{aligned} g(n,0) &= 2^k + 3^k + 4^k + 5^k + 6^k + 7^k + \cdots + n^k \\ g(2n,1) &= 2^k + 3^k + 5^k + 7^k + 9^k + 11^k + \cdots + (2n-1)^k \\ g(6n,2) &= 2^k + 3^k + 5^k + 7^k + 11^k + \cdots + (6n-5)^k + (6n-1)^k \end{aligned}

相当于在这个求和中,g(n,j) 将前 j 个素数的倍数全部去掉了,只留下了那些素数的 k 次方和。

\begin{aligned} g(2n,0) – g(2n,1) &= 4^k + 6^k + 8^k + 10^k + 12^k + \cdots + (2n)^k \\ &= 2^k \left( 2^k + 3^k + 4^k + 5^k + 6^k + \cdots + n^k \right) \\ \\ g(6n,1) – g(6n,2) &= 9^k + 15^k + 21^k + 27^k + 33^k + \cdots + (6n-3)^k \\ &= 3^k \left( 3^k + 5^k + 7^k + 9^k + 11^k + \cdots + (2n-1)^k \right) \end{aligned}

那么从 g(n,j-1) 转移到 g(n,j) 的过程中,我们考虑去掉 p_j^k 的倍数 \sum (p_j \prod p_{j+i}^{e_i})^k

差不多就是有这样一个转移

g(n,j) = g(n,j-1) – p_j^k \left( g\left(\left\lfloor\frac{n}{p_j}\right\rfloor, j-1\right) – \sum_{x=1}^{j-1} p_x^k \right)

为什么要减去前 j-1 个素数 k 次方和的转移呢?因为 (p_x p_j)^k 在这之前都被筛掉了啊。

如果 p_j^2 > n,那肯定没有东西能被删掉了。所以有这样一个转移

g(n,j) = \begin{cases} g(n,j-1), &p_j^2 \gt n \\ g(n,j-1) – p_j^k \left( g(\lfloor n / p_j\rfloor, j-1) – \sum_{x=1}^{j-1} p_x^k \right), &p_j^2 \le n \end{cases}

我们需要一些转换,才能把上面那个东西从递归变成迭代。

我们会发现,我们可能需要经常性的计算 g(\lfloor n/x\rfloor, j) 的值,那我们就一起维护了吧。

这里是预处理的代码。

// euler sieve, spk for sum p^k from p_1 ~ p_pcn
int pri[MAXN], pcn, spk[MAXN];
int m, Sqrt, id1[MAXN], id2[MAXN];
long long gk[MAXN], w[MAXN];

inline long long sum_ik(lld n); // for i=2..n, sum i^k

void init()
{
    m = 0, Sqrt = sqrt(n);

    for (long long L = 1, R; L <= n; L = R+1)
    {
        R = n / (n / L), w[++m] = n / L;
        (w[m] <= Sqrt ? id1[w[m]] : id2[R]) = m;
        gk[m] = sum_ik(w[m]);
    }

    for (int j = 1; j <= pcn; j++) // p_pcn^2 should be greater than n
    {
        for (int i = 1; i <= m && pri[j] <= w[i]/pri[j]; i++)
        {
            int d = (w[i]/pri[j]<=Sqrt) ? id1[w[i]/pri[j]] : id2[n/(w[i]/pri[j])];
            gk[i] -= (spk[j] - spk[j-1]) * (gk[d] - spk[j-1]);
        }
    }
}

其中,id1[i] 表示数论分块中整除值 \lfloor n/l \rfloor = i 的块编号,id2[n/i] 表示整除值 \lfloor n / l \rfloor = i 的块编号,w[s] 表示数论分块第 s 个块的整除值。这样比较节省空间。

然后考虑下面那个迭代。首先,w_i 就是 g 函数的第一个参数。

p_j^2 > w_i,也就是第 i 个块整除值的时候,没有必要继续处理了。内层的 i=1 \cdots m 相当于从大往小枚举块,也就不用利用滚动数组什么的记录上一轮 j-1 时的具体值了。

最后,init执行完时,gk[i] 就是第 i 个整除值 w_i 时,g(w_i, pcn) 的值了,也就是我们想要的东西了。

由于那个多项式可能会有好几个 p^x,所以都维护出来吧。

求 Σf(i)

我们刚才去掉了合数处的点值,那么我们还要根据积性函数的性质把那些点值加回来。

那么我们现在构造这样一个函数

S(n,j) = \sum_{x=2}^n f(x)~[ minp(x) \ge p_j ]

例如

\begin{aligned} S(20,3) &= f(5) + f(7) + f(11) + f(13) + f(17) + f(19) \\ S(20,2) – S(20,3) &= f(3) + f(9) + f(15) \\ S(20,1) – S(20,2) &= f(2) + f(4) + f(6) + \cdots + f(18) + f(20) \end{aligned}

我们需要从最大的质数开始,一个个往回补上点值。

S(n,j) = \sum_{k \ge j} \sum_{e} f(p_k^e)\left(1 + S\left(\left\lfloor n/p_k^e\right\rfloor, k+1\right)\right)

考虑 p_k^e 的倍数,显然当 e=0 的时候不能算上贡献;当 e>0 的时候,\sum p_k^e \prod p_{k+i}^{e_i} 的贡献由 f(p_k^e) 乘上 \sum f(\prod p_{k+i}^{e_i}) 得到。为什么还有个 1+ 呢?S(x,k) 里面可没有 p_j 的东西,不加就算不上 \sum f(p_k^e) 的贡献了呀。

那我们整理一下这个鬼一样的式子。

\begin{aligned} S(n, j) &= \sum_{k \ge j} f(p_k) + \sum_{k \ge j} \sum_{e} \left( f(p_k^e) S(\lfloor n/p_k^e\rfloor, k+1) + f(p_k^{e+1}) \right) \\ &= g(n, |P|) – \sum_{i=1}^{j-1} f(p_i) + \sum_{k \ge j} \sum_{e} \left( f(p_k^e) S(\lfloor n/p_k^e\rfloor) + f(p_k^{e+1}) \right) \end{aligned}

所以这一部分的代码可以写成

long long S(long long x, int y)
{
    if (x <= 1 || pri[y] > x) return 0;
    int d = x <= Sqrt ? id1[x] : id2[n/x];
    long long ans = g(x) - g(pri[y]-1);

    for (int i = y; i <= pcn && pri[i] <= x / pri[i]; i++)
    {
        long long p1 = pri[i], p2 = p1 * p1;
        for (int e = 1; p1 <= x / pri[i]; ++e, p1 = p2, p2 *= pri[i])
            ans += S(x/p1, i+1) * f(i,e) + f(i,e+1);
    }

    return ans;
}

其中 g(n)\sum_{i=1}^n f(i) [i \in prime],也就是上面式子中的 g(n,|P|) 可以从上面预处理的式子中得到。而 f(i,e) 就是 f(p_i^e)

最后你要求的 \sum_{i=1}^n f(i) = S(n,1) + 1

举例来说,我已知积性函数 f(p^e) = C_{e+A-1}^{A-1},那么 f(p)=A,我可以预处理上面预处理部分 sp0, g0 等内容,那么 g(x) - g(pri[y]-1) 应该写为 A * (g0[d]-(y-1)),而 f(i,e) 则写为 C(e+A-1, A-1)。对应例题为 HDU6537。

举例来说,我已知积性函数 f(p^e) = Ae+1,那么 f(p)=A+1,我可以预处理 sp0, g0 等内容,那么 g(x) - g(pri[y]-1) 应该写为 (A+1) * (g0[d]-(y-1)),而 f(i,e) 则写为 A*e+1。对应具体例题为 SPOJ DIVCNT3。

再举例来说,我已知积性函数 f(p^e) = p \oplus e,那么可以在 p=2f(2)=3,然后其他时候有 f(p)=p-1,预处理 sp0, sp1, g0, g1 等内容,那么 g(x) - g(pri[y]-1) 应该写为 g1[d]-g0[d]-sp1[y-1]+y-1,并特判 if (y == 1) ans += 2;,而 f(i,e) 则写为 pri[i]^e。对应具体例题为 LOJ6053。

后记

所以说这是不是很像埃式筛啊,所以它又被称为拓展埃式筛。

据说时间复杂度是

O\left( \frac{n^{\frac{3}{4}}}{\log n} \right)

咱也不会证,咱也不敢证。

2019杭电多校第三场 E. Easy Math Problem

题目想要对于给定 nk

\sum_{i=1}^n \sum_{j=1}^n \gcd(i,j)^k ~ lcm(i,j) ~[\gcd(i,j) \in prime]

那我们考虑一下枚举某个素数 d = \gcd(i,j) 的情况吧。

\sum_{d \in P} d^{k-1} \sum_{i=1}^n \sum_{j=1}^n ij ~ [\gcd(i,j) = d]

我们先来考虑一下

f(N) = \sum_{i=1}^N \sum_{j=1}^N ij [\gcd(i,j)=1]

先根据对称性,可以得到

\sum_{i=1}^N \sum_{j=1}^i ij [\gcd(i,j)=1] = \frac{f(N)+1}{2}

再根据非常显然的 \gcd(i,j)=1 \Leftrightarrow \gcd(i,j-i)=1,可以配对

\sum_{j=1}^i j [\gcd(i,j)=1] = \frac{i \varphi(i) + [i=1]}{2}

那么就有

\begin{aligned} \frac{f(N)+1}{2} &= \sum_{i=1}^N \sum_{j=1}^i ij [\gcd(i,j)=1] \\ &= \frac{\sum_{i=1}^N i^2 \varphi(i) +1}{2} \end{aligned}

所以说刚才求的 f(N) 就相当于对 i^2 \varphi(i) 求了个前缀和。根据杜教筛一般套路,给配上个 g(n) = n^2 就可以筛筛筛了。

f(n) = \frac{n(n+1)(2n+1)}{6} – \sum_{i=2}^n i^2 f\left(\left\lfloor\frac{n}{i}\right\rfloor\right)

再回到原式,其实相当于求

\sum_{d \in P} d^{k+1} \sum_{i=1}^{\lfloor n/d\rfloor} i^2 \varphi(i)

那么我们对 n 进行数论分块,就可以获得答案了。考虑块 [L,R] 内前面一个东西贡献的值,其实相当于

\sum_{i=L}^R i^{k+1} [i \in prime]

那么把它变成一个前缀和的差,构造一个积性函数

h(p) = p^{k+1}; h(p^e) = 0, e>1

然后考虑用min25筛,就可以做了。

由于min25筛需要快速求出 p^k 的前缀和,那么就得上拉格朗日插值一下。

推荐看 这篇文章 学习min25。我也要好好花时间学习一下啦 QAQ

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MAXN = 7e6+5;
const int MOD = 1e9+7;
lld N, K, ans;

lld fpow(lld a, lld k) {
    lld b = 1;
    while (k) {
        if (k & 1) b = b * a % MOD;
        a = a * a % MOD;
        k >>= 1;
    }
    return b;
}

namespace interpolation
{
    typedef vector<lld> Poly;

    Poly& operator*=(Poly &a, const Poly &b) {
        Poly c(a.size()+b.size()-1, 0);
        for (int i = 0; i < a.size(); i++)
            for (int j = 0; j < b.size(); j++)
                (c[i+j] += a[i] * b[j] % MOD) %= MOD;
        return a = c;
    }

    Poly& operator*=(Poly &a, lld qwq) {
        for (lld &ai : a)
            ai = ai * qwq % MOD;
        return a;
    }

    Poly& operator+=(Poly &a, const Poly &b) {
        a.resize(max(a.size(), b.size()));
        for (int i = 0; i < b.size(); i++)
            (a[i] += b[i]) %= MOD;
        return a;
    }

    Poly Lagrange(const Poly &X, const Poly &FX) {
        Poly c;
        int n = X.size();
        for (int i = 0; i < n; i++) {
            Poly x({FX[i]});
            for (int j = 0; j < n; j++)
                if (j != i)
                    x *= Poly({-X[j],1}),
                    x *= fpow(X[i]-X[j], MOD-2);
            c += x;
        }
        return c;
    }

    lld evaluate(const Poly &a, lld x) {
        lld xx = 1, ans = 0;
        for (int i = 0; i < a.size(); i++)
            (ans += a[i] * xx % MOD) %= MOD,
            xx = xx * x % MOD;
        return ans;
    }

    Poly init_lag() {
        Poly a, b;
        lld last = 0;
        for (int i = 0; i <= K+2; i++) {
            a.push_back(i);
            (last += fpow(i, K+1)) %= MOD;
            b.push_back(last);
        }
        return Lagrange(a, b);
    }
}

namespace euler_seive
{
    int pri[MAXN], _B1[MAXN], pcn;
    bool isnp[MAXN];
    unordered_map<lld, int> _B2;
    const lld inv6 = fpow(6, MOD-2);

    void init_prime()
    {
        _B1[1] = 1;

        for (int i = 2; i < MAXN; i++)
        {
            if (!isnp[i])
            {
                pri[++pcn] = i;
                _B1[i] = i-1;
            }

            for (int j = 1; j <= pcn; j++)
            {
                if (1ll * i * pri[j] >= MAXN) break;
                isnp[i * pri[j]] = true;
                if (i % pri[j] == 0) {
                    _B1[i * pri[j]] = _B1[i] * pri[j];
                    break;
                } else {
                    _B1[i * pri[j]] = _B1[i] * (pri[j]-1);
                }
            }
        }

        for (int i = 2; i < MAXN; i++)
            _B1[i] = (_B1[i-1] + 1ll * _B1[i] * i % MOD * i % MOD) % MOD;
    }

    inline lld sum2(lld a) { a %= MOD; return a * (a+1) % MOD * (2*a+1) % MOD * inv6 % MOD; }

    int sumId2Phi(lld n)
    {
        if (n < MAXN) return _B1[n];
        else if (_B2.count(n)) return _B2[n];
        lld ans = n % MOD * (n%MOD+1) % MOD * ((MOD+1)/2) % MOD;
        ans = ans * ans % MOD; // \sum n^3
        for (lld l = 2, r; l <= n; l = r+1)
            r = n / (n / l),
            ans -= (sum2(r) - sum2(l-1)) * sumId2Phi(n/l) % MOD;
        ans = (ans % MOD + MOD) % MOD;
        return _B2[n] = ans;
    }
}

namespace min25_seive
{
    using euler_seive::pri;
    using euler_seive::pcn;
    interpolation::Poly _sumik_f;
    int _sumpk[MAXN], m, Sqrt;
    int id1[MAXN], id2[MAXN], tot;
    lld g[MAXN], w[MAXN];

    inline lld sumik(lld n) {
        return interpolation::evaluate(_sumik_f, n % MOD) - 1;
    }

    inline void init()
    {
        for (int i = 1; i <= pcn; i++)
            _sumpk[i] = (_sumpk[i-1] + fpow(pri[i], K+1)) % MOD;
        _sumik_f = interpolation::init_lag();

        m = 0;
        Sqrt = sqrt(N);

        for (lld l = 1, r; l <= N; l = r+1)
        {
            r = N / (N / l), w[++m] = N / l;
            (w[m] <= Sqrt ? id1[w[m]] : id2[N/w[m]]) = m;
            g[m] = sumik(w[m]);
        }

        tot = lower_bound(pri+1, pri+1+pcn, Sqrt) - pri;
        for (int j = 1; j <= tot; j++)
        {
            for (int i = 1; i <= m && pri[j] <= w[i]/pri[j]; i++)
            {
                int k = (w[i]/pri[j]<=Sqrt) ? id1[w[i]/pri[j]] : id2[N/(w[i]/pri[j])];
                (g[i] -= (_sumpk[j] - _sumpk[j-1]) * (g[k] - _sumpk[j-1]) % MOD) %= MOD;
            }
        }
    }

    inline lld sumpk(lld n) {
        return g[n <= Sqrt ? id1[n] : id2[N/n]];
    }
}

void solve()
{
    scanf("%lld %lld", &N, &K);
    min25_seive::init();

    ans = 0;
    for (lld l = 1, r, s = 0; l <= N; l = r+1)
    {
        r = N / (N / l);
        lld ss = min25_seive::sumpk(r);
        ans += euler_seive::sumId2Phi(N/l) * (ss-s) % MOD;
        s = ss;
    }

    printf("%lld\n", (ans % MOD + MOD) % MOD);
}

int main()
{
    euler_seive::init_prime();
    int T;
    scanf("%d", &T);
    while (T--) solve();
    return 0;
}

我好菜啊,反演不会推,min25也不会写,一道题写了四个小时到处是bug,到处TLE和MLE。

std好快啊。qwq

2019杭电多校第一场 K. Function

题目想要求

\sum_{i=1}^n \gcd\left(\left\lfloor\sqrt[3]i\right\rfloor,i\right) mod~998244353

首先可以想到,利用 \lfloor\sqrt[3]{i}\rfloor 的不同值进行分块。我们先记 m=\lfloor\sqrt[3]{n}\rfloor, 那么原式变为

\sum_{c=1}^{m} \sum_{i=c^3}^{\min((c+1)^3-1, n)} \gcd(c,i)

看着那个min比较碍眼,不如考虑当 c 的情况,那么有

H(c)=\sum_{i=c^3}^{(c+1)^3-1} \gcd(c, i) = \sum_{d|c} d \sum_{i=c^3}^{(c+1)^3-1} \left[\gcd(c,i)=d\right]

看起来有个什么类似于前缀和的东西?再设

f(n,d) = \sum_{i=1}^{n} \left[ \gcd(d,i) = 1 \right]

根据莫比乌斯函数的容斥意义容易得到

f(n,d) = \sum_{g|d} \mu(g) \left\lfloor\frac{n}{g}\right\rfloor

回到上面那个式子

\begin{aligned} H(c) &= \sum_{d|c} d \times \left( f\left(\left\lfloor\frac{(c+1)^3-1}{d}\right\rfloor,\frac{c}{d}\right) – f\left(\left\lfloor\frac{c^3-1}{d}\right\rfloor,\frac{c}{d}\right) \right) \\ &= \sum_{d|c} d \sum_{g|\frac{c}{d}} \mu(g) \left( \left\lfloor\frac{(c+1)^3-1}{gd}\right\rfloor – \left\lfloor\frac{c^3-1}{gd}\right\rfloor \right) \\ &= \sum_{gkd=c} d ~\mu(g) \left( \left\lfloor\frac{(c+1)^3-1}{gd}\right\rfloor – \left\lfloor\frac{c^3-1}{gd}\right\rfloor \right) \\ &= \sum_{kt=c} \left( \sum_{gd=t} d ~\mu(g) \right) \left( \left\lfloor\frac{(c+1)^3-1}{gd}\right\rfloor – \left\lfloor\frac{c^3-1}{gd}\right\rfloor \right) \\ &= \sum_{kt=c} \varphi(t) \left( \left\lfloor\frac{(c+1)^3-1}{t}\right\rfloor – \left\lfloor\frac{c^3-1}{t}\right\rfloor \right) \end{aligned}

好像可以发现 c | (c+1)^3-1,那后面那个玩意又有点碍眼,看看好像只有 t=1 的时候不太一样?不如我们……

\begin{aligned} H(c) &= \sum_{kt=c} \varphi(t) \left( \frac{(c+1)^3-1}{t} – \left\lfloor\frac{c^3}{t} \right\rfloor + 1 \right) \\ &= \sum_{kt=c} \varphi(t)~k~(3c+3) + \sum_{kt=c} \varphi(t)k \\ &= c + 3(c+1) (\varphi * id)(c) \end{aligned}

那我们继续看后面那部分吧

\begin{aligned} H'(m) &= \sum_{d|m} d \sum_{i=m^3}^n [\gcd(m,i) = d] \\ &= \sum_{d|m} d \sum_{g|\frac{m}{d}} \mu(g) \left( \left\lfloor\frac{n}{gd}\right\rfloor – \left\lfloor\frac{m^3-1}{gd}\right\rfloor \right) \\ &= \sum_{t|m} \varphi(t) \left( \left\lfloor\frac{n}{t}\right\rfloor – \left\lfloor\frac{m^3-1}{t}\right\rfloor \right) \end{aligned}

再回去看看最开始要的东西

H'(m) + \sum_{i=1}^{m-1} H(i)

不过 H(i) 看起来可以线性筛再前缀和?

好像 O(n^{\frac{1}{6}}) 单组询问的复杂度可以接受?Let’s go~

#include <bits/stdc++.h>
typedef long long lld;
using namespace std;
const size_t MAXN = 1e7+5;
const int MOD = 998244353;
bool isnp[MAXN];
int pcn, pri[MAXN/10];
lld g[MAXN], s[MAXN];
int phi[MAXN], q[MAXN];
int mod;

void init_prime()
{
    g[1] = q[1] = phi[1] = 1;

    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i])
        {
            g[i] = s[i] = 2*i-1;
            q[i] = i-1;
            phi[i] = i-1;
            pri[pcn++] = i;
        }

        for (int j = 0; j < pcn; j++)
        {
            int k = i * pri[j];
            if (k >= MAXN) break;
            isnp[k] = true;

            if (i % pri[j] == 0)
            {
                s[k] = (s[i] + q[i] + 0ll) * pri[j];
                q[k] = q[i] * pri[j];
                g[k] = g[i] / s[i] * s[k];
                phi[k] = phi[i] * pri[j];
                break;
            }
            else
            {
                phi[k] = phi[i] * (pri[j] - 1);
                q[k] = q[pri[j]];
                s[k] = s[pri[j]];
                g[k] = g[i] * s[k];
            }
        }
    }

    for (int i = 1; i < MAXN; i++)
    {
        g[i] = (g[i-1] + i + (3 * i + 3) * (g[i] % MOD)) % MOD;
    }
}

int two(__int128 n)
{
    if (n < 8) return 1;
    int l = 1, r = 1e9+7, m;
    while (r - l > 1)
    {
        m = (l + r) >> 1;
        if (m > n / m / m)
            r = m;
        else
            l = m;
    }
    return r-1;
}

void solve()
{
    __int128 n = 0;
    char buf[30];
    scanf("%s", buf);
    for (int i = 0; buf[i]; i++)
        n = n * 10 + (buf[i] - '0');
    int m = two(n);
    lld ans = g[m-1];
    __int128 qwq = m; qwq *= m; qwq *= m; qwq -= 1;
    for (int i = 1; i * i <= m; i++)
    {
        if (m % i == 0)
        {
            ans += phi[i] * ((n / i % MOD) - (qwq / i % MOD)) % MOD;
            if (m != i * i)
                ans += phi[m / i] * (n / (m / i) % MOD - (qwq / (m / i) % MOD)) % MOD;
        }
    }
    ans %= MOD;
    if (ans < 0) ans += MOD;
    printf("%lld\n", ans);
}

int main()
{
    init_prime();
    int T;
    scanf("%d", &T);
    while (T--) solve();
    return 0;
}

2019CCPC湘潭邀请赛 F. Neko and function

题目链接:HDU-6537

题目有 f(n,k) 的定义,但是显然不够好直接计算。不妨先将限定条件改为 a_i \ge 1,考虑 g(n,k)

此时就可以发现,每个素数和其幂次对 g(n,k) 的贡献是独立的,也就是说,g(n,k)是个积性函数。所以,我们来考虑 g(p^e, k)。这个过程中,相当于将 e 个小球放进 k 个盒子,并且盒子有序可空,球无序。

g(n, k) = \prod_{p^e || n} C_{e+k-1}^{k-1}

用容斥原理去除重复项。考虑长度为 s\lbrace a_n \rbrace 存在 1 的方案对 g(n,x) 的贡献。那么,在长度为 k 的数组里选择 x 个放 1,剩下的将之前方案排列进去,也就是 C_k^x

f(n, k) = \sum_{x=0}^k (-1)^{k-x} g(n, x) C_k^x

所以

\sum_{i=1}^n f(n, k) = \sum_{x=0}^k (-1)^{k-x} C_k^x \sum_{i=1}^n g(n, x)

那么用 min_25 筛就可以解决 g(n,k) 的前缀和啦,再特判一些可能浪费时间的样例。

吐槽一下,这个时限卡的太狠了,用long long过不了。本身就是个人类智慧题,竟然还要卡常。

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 5e6+5;
const int MOD = 1e9+7;
bool isnp[MAXN];
int pri[MAXN], Pcn, id1[MAXN], id2[MAXN], m;
int w[MAXN], n, Sqrt, h[MAXN]; // h:sum 1
int inv[MAXN], invs[MAXN], fac[MAXN], pcn, K, kk;
inline int mul(int a, int b) { return int(1ll * a * b % MOD); }
inline int mul(int a, int b, int c) { return int(1ll * a * b % MOD * c % MOD); }
inline void modu(int &a) { while (a >= MOD) a -= MOD; while (a < 0) a += MOD; }
inline int add(int a, int b) { modu(a); modu(b); modu(a += b); return a; }
inline int add(int a, int b, int c) { return add(a, add(b, c)); }
inline int C(int N, int M) { return N < 0 || M > N || M < 0 ? 0 : mul(fac[N], invs[M], invs[N-M]); }

void sieve()
{
    inv[1] = invs[0] = fac[0] = invs[1] = fac[1] = 1;

    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i]) pri[++Pcn] = i;
        inv[i] = mul(MOD-MOD/i, inv[MOD%i]);
        fac[i] = mul(fac[i-1], i);
        invs[i] = mul(invs[i-1], inv[i]);
        for (int j = 1; j <= Pcn && i * pri[j] < MAXN; j++)
        {
            isnp[i * pri[j]] = true;
            if (i % pri[j] == 0) break;
        }
    }
}

int S(int x, int y)
{
    if (x <= 1 || pri[y] > x) return 0;
    int k = (x <= Sqrt) ? id1[x] : id2[n/x];
    int ans = mul(add(h[k], MOD - y + 1), K);

    for (int i = y; i <= pcn && pri[i] <= x / pri[i]; i++)
    {
        int p1 = pri[i], p2 = p1 * p1;
        for (int e = 1; p1 <= x / pri[i]; ++e, p1 = p2, p2 *= pri[i])
            ans = add(ans, mul(S(x/p1, i+1), C(e+K-1, K-1)), C(e+K, K-1));
    }

    return ans;
}

inline int g(int x)
{
    if (K == 0) return 1;
    else if (K == 1) return x % MOD;
    else return S(n,1) + 1;
}

inline int calc()
{
    if (kk == 1) return (n-1) % MOD;
    else if (n < (1<<kk)) return 0;
    else if (n < (1<<kk)+(1<<(kk-1))) return 1;

    Sqrt = sqrt(n);
    for (pcn = 1, m = 0; pri[pcn] <= Sqrt; pcn++);

    for (int l = 1, r; l <= n; l = r+1)
    {
        r = n / (n / l); w[++m] = n / l;
        if (w[m] <= Sqrt) id1[w[m]] = m; else id2[n/w[m]] = m;
        h[m] = add(w[m], MOD - 1);
    }

    for (int j = 1; j <= pcn; j++)
    {
        for (int i = 1; i <= m && pri[j] <= w[i] / pri[j]; i++)
        {
            int k = (w[i]/pri[j]<=Sqrt)?id1[w[i]/pri[j]]:id2[n/(w[i]/pri[j])];
            h[i] = add(h[i], MOD - h[k], j - 1);
        }
    }

    int ans = 0;
    for (K = 0; K <= kk; K++)
        ans = add(ans, mul(g(n), C(kk, K), ((kk-K)&1) ? MOD-1 : 1));
    return (ans % MOD + MOD) % MOD;
}

int main()
{
    sieve();
    while (scanf("%d %d", &n, &kk) == 2)
        printf("%d\n", calc());
    return 0;
}

vjudge 的 Leaderboard 上还有个神仙做法,看懂后更新。

2019CCPC东北地区赛 I. Temperature Survey

题目链接:Gym 102220I

这个题目的动态规划解法应该挺好想的。

dp_{i,j} = \sum_{k=1}^{j} dp_{i,k} = dp_{i,j-1} + dp_{i-1,j}, k \in [1,a_i]

然而这显然是 O(n^2) 的做法,过不了这题。

那么我们来看看quailty的Editorial啦~

这里仅仅补充一下关于矩形的转移的部分。

考虑宽度为 n,高度为 m 的一个矩形的二维前缀和

显然由于对称性,我们可以把a和b分开计算。

先考虑如何得到下底边的值。

\begin{aligned} H(x) & = \sum_{i=0}^{n-1} \left( \sum_{j=0}^i C_{m-1+j}^j a_{i-j} \right) x^i \\ & = \left( \sum_{i=0}^{n-1} a_i x^i \right) \left( \sum_{j=0}^{n-1}C_{m-1+j}^j x^j \right) \end{aligned}

再考虑右边的值

V(x) = \sum_{i=0}^{m-1} \left( \sum_{j=0}^{n-1} C_{n-1-j+i}^{i} a_j \right) x^i

先考虑

x^{n-1}V(x) = \sum_{i=n-1}^{n+m-2} \left( \sum_{j=0}^{n-1} C_{i-j}^{i-n+1} a_j \right) x^i

a_n, a_{n+1}, \cdots = 0,然后补充几项,可以得到

\begin{aligned} x^{n-1}V(x) & = \sum_{i=n-1}^{n+m-2} \left( \sum_{j=0}^i C_{i-j}^{i-n+1} a_j \right) x^i \\ & = \sum_{k=n-1}^{n+m-2} \left( \sum_{j=0}^k \frac{(k-j)!}{(n-1-j)!} a_j \right) \frac{x^k}{(k-n+1)!} \end{aligned}

然后让我们换个式子再看看这个看起来可以NTT的东西

\begin{aligned} W(x) & = \sum_{k=n-1}^{n+m-2} \left( \sum_{i+j=k} i! \frac{a_j}{(n-1-j)!} \right) x^k \\ & = \left( \sum_{j=0}^{n-1} \frac{a_j}{(n-1-j)!} x^j \right) \left( \sum_{i=0}^{n+m-2} i!x^i \right) \end{aligned}

就可以很快的计算这个方格的右腰和下底啦~

也是一种类似于分治NTT的思想呢。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MOD = 998244353, L=20, MAXL = 1<<L;
int rev[MAXL], wmk[MAXL], inv[MAXL], len, N;
int fac[MAXL], invs[MAXL];

void init()
{
    wmk[0] = inv[1] = fac[1] = 1;
    invs[1] = fac[0] = invs[0] = 1;
    long long t = 15311432; //G^119
    for (int i = 0; i < 23-L; i++)
        t = t * t % MOD;
    for (int i = 1; i < MAXL; i++)
        wmk[i] = wmk[i-1] * t % MOD;
    for (int i = 2; i < MAXL; i++)
        fac[i] = 1ll * fac[i-1] * i % MOD;
    for (int i = 2; i < MAXL; i++)
        inv[i] = 1ll * (MOD-MOD/i) * inv[MOD%i] % MOD;
    for (int i = 2; i < MAXL; i++)
        invs[i] = 1ll * invs[i-1] * inv[i] % MOD;
}

int combine(int n, int m)
{
    return 1ll * fac[n] * invs[m] % MOD * invs[n-m] % MOD;
}

void discreteFourierTransform(vector<int> &a)
{
    if (!wmk[1]) init();
    for (int i = 0; i < N; i++)
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int m = 2, m2 = 1; m <= N; m <<= 1, m2 <<= 1)
        for (int k = 0; k < N; k += m)
            for (int j = 0, t, u; j < m2; j++)
                t = 1ll * wmk[MAXL/m*j] * a[k+j+m2] % MOD,
                u = a[k+j], a[k+j] = (u+t)%MOD, a[k+j+m2] = (u-t+MOD)%MOD;
}

void multiply(vector<int> &a, vector<int> b)
{
    int need = int(a.size() + b.size() - 1);

    if (need <= 128)
    {
        vector<int> c = a;
        a.assign(need, 0);
        for (int i = 0; i < int(c.size()); i++)
            for (int j = 0; j < int(b.size()); j++)
                a[i+j] = (a[i+j] + 1ll * c[i] * b[j]) % MOD;
    }
    else
    {
        len = 0, N = 1;
        while (N < need) ++len, N <<= 1;
        for (int i = 0; i < N; i++)
            rev[i] = (rev[i>>1]>>1)|((i&1)<<(len-1));
        if (a.size() < N) a.resize(N);
        if (b.size() < N) b.resize(N);
        bool equals_ab = a == b;
        discreteFourierTransform(a);
        if (equals_ab) b = a;
        else discreteFourierTransform(b);
        for (int i = 0; i < N; i++)
            a[i] = int(1ll*a[i]*b[i]%MOD*inv[N]%MOD);
        reverse(a.begin()+1, a.begin()+N);
        discreteFourierTransform(a);
        a.resize(need);
    }
}

struct solver : map<pair<int,int>,int>
{
    int operator()(int a, int b) const
    {
        if (a < 0 || b < 0) return 0;
        return this->at(make_pair(a, b));
    }

    int &operator()(int a, int b)
    {
        return (*this)[make_pair(a, b)];
    }
} dp;

void solveRect(int x, int y, int width, int height)
{
    vector<int> bottom(width);
    vector<int> right(height);

    // to solve the top to bottom
    {
        vector<int> A(width);
        for (int i = 0; i < width; i++)
            A[i] = dp(y-1, x+i);
        vector<int> B(width);
        for (int i = 0; i < width; i++)
            B[i] = combine(height-1+i, i);
        multiply(A, B);
        for (int i = 0; i < width; i++)
            bottom[i] += A[i];
    }

    // to solve the top to right
    {
        vector<int> A(width);
        for (int i = 0; i < width; i++)
            A[i] = 1ll * dp(y-1, x+i) * invs[width - 1 - i] % MOD;
        vector<int> B(width + height - 1);
        for (int i = 0; i < width + height - 1; i++)
            B[i] = fac[i];
        multiply(A, B);
        for (int i = 0; i < height; i++)
            right[i] += 1ll * A[i+width-1] * invs[i] % MOD;
    }

    // to solve the left to right
    {
        vector<int> A(height);
        for (int i = 0; i < height; i++)
            A[i] = dp(y+i, x-1);
        vector<int> B(height);
        for (int i = 0; i < height; i++)
            B[i] = combine(width-1+i, i);
        multiply(A, B);
        for (int i = 0; i < height; i++)
            right[i] += A[i];
    }

    // to solve the left to bottom
    {
        vector<int> A(height);
        for (int i = 0; i < height; i++)
            A[i] = 1ll * dp(y+i, x-1) * invs[height - 1 - i] % MOD;
        vector<int> B(height + width - 1);
        for (int i = 0; i < height + width - 1; i++)
            B[i] = fac[i];
        multiply(A, B);
        for (int i = 0; i < width; i++)
            bottom[i] += 1ll * A[i+height-1] * invs[i] % MOD;
    }

    for (int i = 0; i < width; i++)
        dp(y+height-1, x+i) = bottom[i] % MOD;
    for (int i = 0; i < height; i++)
        dp(y+i, x+width-1) = right[i] % MOD;
}

const int MAXN = 2e5+5;
int a[MAXN];

void solve(int l, int r, int b)
{
    while (l <= r && a[l] < b) l++;
    if (l > r) return;

    int m = (l + r) >> 1;
    solve(l, m-1, b);
    printf("solveRect(%d, %d, %d, %d);\n", b, m, a[m]-b+1, r-m+1);
    solveRect(b, m, a[m]-b+1, r-m+1);
    solve(m+1, r, a[m]+1);
}

int main()
{
    init();
    int T, n;
    scanf("%d", &T);

    while (T--)
    {
        dp.clear();
        dp(0, 1) = 1;
        scanf("%d", &n);
        for (int i = 1; i <= n; i++)
            scanf("%d", &a[i]);
        a[n+1] = n;
        solve(1, n+1, 1);
        printf("%d\n", dp(n+1, n));
    }

    return 0;
}

2019ICPC西安邀请赛 B. Product

题目链接:计蒜客

对于给定的 n, m, p,我们要计算

\prod_{i=1}^n \prod_{j=1}^n \prod_{k=1}^n m^{\gcd(i,j) [k|\gcd(i,j)]} \mod p

首先,我们可以把乘积扔到指数上,变成求和;然后由于 p 是素数,我们求其指数模 p-1

\sum_{i=1}^n \sum_{j=1}^n \sum_{k=1}^n \gcd(i,j) [k|\gcd(i,j)] \mod p-1

考虑先对 k 求和,则此式可以变成

f(n) = \sum_{k=1}^n k \left(\sum_{i=1}^{\lfloor n/k\rfloor} \sum_{j=1}^{\lfloor n/k\rfloor} \gcd(i,j) \right)

考虑这样的一个函数,并使用莫比乌斯反演、适当的求和变换和狄利克雷卷积的结论

\begin{aligned} g(N) & = \sum_{i=1}^N \sum_{j=1}^N \gcd(i,j) \\ & = \sum_{s=1}^N s \sum_{i=1}^N \sum_{j=1}^N [\gcd(i,j) = s] \\ & = \sum_{s=1}^N s \sum_{t=1}^{\lfloor N/s \rfloor} \left\lfloor\frac{N}{st}\right\rfloor^2 \mu(t) \\ & = \sum_{d=1}^N \left\lfloor\frac{N}{d}\right\rfloor^2 \sum_{st=d} \mu(t) s \\ & = \sum_{d=1}^N \left\lfloor\frac{N}{d}\right\rfloor^2 \varphi(d) \end{aligned}

回到原式,再使用求和变换

\begin{aligned} f(n) & = \sum_{k=1}^n k \sum_{t=1}^{\lfloor n/k \rfloor} \left\lfloor\frac{n}{kt}\right\rfloor^2 \varphi(k) \\ & = \sum_{d=1}^n \left\lfloor\frac{n}{d}\right\rfloor^2 \sum_{kt=d} \varphi(k)t \end{aligned}

如果后面那个 \varphi * id 可以用杜教筛计算出来,那么再预处理部分前缀和,我们就能 O(n^{2/3}) 计算出结果啦!

既然杜教筛要配一个函数使得其狄利克雷卷积好计算,那么…我们好像只能配上 1 了。

(1 * \varphi * id)(n) = (1 * \mu * id * id)(n) = (id * id)(n) = n \times d(n)

好吧这个式子好像不太好求前缀和?不过考虑一下 d(n) 求前缀和的方法,我们好像还是可以求的。

\sum_{t=1}^n t d(t) = \sum_{t=1}^n \sum_{sk=t} sk = \sum_{s=1}^n s \sum_{k=1}^{\lfloor n/s \rfloor}k

行吧,三个数论分块,写吧…馍馍片。

#include <bits/stdc++.h>
typedef long long lld;
using namespace std;
const size_t MAXN = 1e6+5;
bool isnp[MAXN];
int pcn, pri[MAXN];
lld g[MAXN], q[MAXN], s[MAXN];
lld _P1[MAXN];
map<lld,lld> _P2;
int mod;

void init_prime()
{
    _P1[1] = g[1] = q[1] = 1;

    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i])
        {
            g[i] = s[i] = 2*i-1;
            q[i] = i-1;
            pri[pcn++] = i;
        }

        for (int j = 0; j < pcn; j++)
        {
            int k = i * pri[j];
            if (k >= MAXN) break;
            isnp[k] = true;

            if (i % pri[j] == 0)
            {
                s[k] = (s[i] + q[i] + 0ll) * pri[j];
                q[k] = q[i] * pri[j];
                g[k] = g[i] / s[i] * s[k];
                break;
            }
            else
            {
                q[k] = q[pri[j]];
                s[k] = s[pri[j]];
                g[k] = g[i] * s[k];
            }
        }

        _P1[i] = _P1[i-1] + g[i];
    }
}

inline lld sumn(lld n)
{
    return n*(n+1)/2%mod;
}

lld sum_xdx(lld n)
{
    lld ans = 0;

    for (lld l = 1, r; l <= n; l = r+1)
    {
        r = n / (n / l);
        (ans += (sumn(r) - sumn(l-1)) * sumn(n/l) % mod) %= mod;
    }

    return ans;
}

lld G(lld n)
{
    if (n < MAXN) return _P1[n] % mod;
    else if (_P2.count(n)) return _P2[n];

    lld ans = sum_xdx(n);
    for (lld k = 2, fk; k <= n; k = fk+1)
    {
        fk = n / (n/k);
        ans -= (fk-k+1) * G(n/k) % mod;
    }

    ans = (ans % mod + mod) % mod;
    return _P2[n] = ans;
}

lld fpow(lld b, lld n, lld MOD)
{
    lld a = 1;

    while (n)
    {
        if (n & 1) a = a * b % MOD;
        b = b * b % MOD;
        n >>= 1;
    }

    return a;
}

inline lld sumn2(lld n)
{
    __int128 s = n*2+1;
    s *= n * (n+1);
    s /= 6;
    return lld(s % mod);
}

int main()
{
    lld n, m, p;
    scanf("%lld %lld %lld", &n, &m, &p);
    mod = p-1;
    init_prime();
    G(n);

    lld fff = 0;
    for (lld l = 1, r; l <= n; l = r+1)
    {
        r = n / (n / l);
        (fff += (G(r) - G(l-1)) * ((n/l) * (n/l) % mod)) %= mod;
    }

    if (fff < 0) fff += mod;
    printf("%lld\n", fpow(m, fff, p));
    return 0;
}

上面式子中,为了求某函数前缀和,写了一堆数组,在这里解释一下:

首先考虑 g(p^n),则有

\begin{aligned} g(p^n) & = \sum_{ij=p^n} \varphi(i) j = \sum_{i=0}^n \varphi(p^i) p^{n-i} \\ & = p^n + n\varphi(p^n) = (n+1)p^n – np^{n-1} \end{aligned}

于是不停的利用最小因子和它的积性去线性筛就好啦~

  • g[n] 表示 g(n)
  • q[n] 表示 n 最小因子 p 以及其对应指数 e 所对应的 (e+1)p^e – e p^{e-1}
  • s[n] 表示 n 最小因子 p 以及其对应指数 e 所对应的 g(p^e)

好麻烦啊QAQ

2019ICPC西安邀请赛 I. Cracking Password

题目链接:计蒜客

题目意思是,已有 a, p 两个数字,给出 \lbrace a^k \mod p \rbrace 这个数列中的连续几项,问是否能够求出 a, p,是不存在,一解,还是多解?

首先先来枚举 p 咯。我们可以发现,p | b_{i-1}b_{i+1}-b_i^2。于是枚举这些gcd的质因数,然后判断是否有对应的a存在。

如果 n 只有一个,显然是unsure的。

如果这个数列完全为不相等的等比数列,那么由于模数任取,总可以找到一个模数符合要求并且其原根的某个幂次为 a,显然我们能给出无穷多组构造,是unsure的。

如果这个数列全部相等,如果为全 1 则unsure,如果全 非1 则error。

如果 n 有两个,那么也可以根据上面两条的得出类似unsure的结论。

最后的坑点在于,是否已知的数字在那个真实的数列里。例如 n=3, b=\lbrace3, 6, 5\rbrace,直接计算可以得到a=2, p=7,但是数列却是 \lbrace 2, 4, 1 \rbrace,所以也应当是error,利用离散对数判断一下即可。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MAXN = 1e5+5;
lld b[MAXN], n, maxb = 0, minb = 100000000000;
int pri[MAXN], pcn;
bool isnp[MAXN];

lld gcd(lld a, lld b) { return b ? gcd(b, a % b) : a; }

lld mul(lld a, lld b, lld p)
{
    __int128 aa = a;
    aa *= b;
    aa %= p;
    return lld(aa);
}

lld fpow(lld a, lld n, lld p)
{
    lld b = 1;
    assert(n >= 0);

    while (n)
    {
        if (n & 1) b = mul(b, a, p);
        a = mul(a, a, p);
        n >>= 1;
    }

    return b;
}

lld bsgs(lld a, lld b, lld p)
{
    map<lld,lld> treemap;
    lld m = ceil(sqrt(p+0.5));
    a %= p, b %= p;
    lld inva = fpow(a, m, p);

    for (lld i = 0, ans = b; i <= m; i++)
    {
        ans = i==0 ? b : mul(ans, a, p);
        treemap[ans] = i;
    }

    for (lld i = 1, ans = 1; i <= m; i++)
    {
        ans = mul(ans, inva, p);
        if (treemap.count(ans))
        {
            lld ret = i * m - treemap[ans];
            ret = (ret % p + p) % p;
            return ret;
        }
    }

    return -1;
}

pair<lld,lld> solve(lld p)
{
    lld f = fpow(b[0], p-2, p);
    lld a = mul(b[1], f, p);
    for (int i = 2; i < n; i++)
        if (b[i] != mul(b[i-1], a, p))
            return make_pair(-1, -1);
    if (bsgs(a, b[0], p) == -1)
        return make_pair(-1, -1);
    return make_pair(a, p);
}

void init_prime()
{
    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i]) pri[pcn++] = i;
        for (int j = 0; j < pcn; j++)
        {
            int k = i * pri[j];
            if (k >= MAXN) break;
            isnp[k] = true;
            if (i % pri[j] == 0) break;
        }
    }
}

int main()
{
    init_prime();
    scanf("%lld", &n);

    for (int i = 0; i < n; i++)
    {
        scanf("%lld", &b[i]);
        maxb = max(maxb, b[i]);
        minb = min(minb, b[i]);
    }

    if (n == 1)
    {
        printf("unsure\n");
    }
    else if (n == 2)
    {
        if (b[0] == b[1] && b[0] != 1)
        {
            // for the sequence should be part of a^n,
            // so this will not be legal.
            printf("error\n");
        }
        else
        {
            printf("unsure\n");
        }
    }
    else
    {
        set<lld> possibleP;
        lld wtf = 0;
        for (int i = 2; i < n; i++)
            wtf = abs(gcd(b[i] * b[i-2] - b[i-1] * b[i-1], wtf));

        for (int i = 0; i < pcn; i++)
        {
            if (1ll * pri[i] * pri[i] > wtf) break;

            if (wtf % pri[i] == 0)
            {
                if (pri[i] > maxb) possibleP.insert(pri[i]);
                while (wtf % pri[i] == 0) wtf /= pri[i];
            }
        }

        if (wtf > maxb) possibleP.insert(wtf);

        if (possibleP.empty())
        {
            if (!wtf && minb == maxb)
            {
                // all equal, look at 1
                printf(minb != 1 ? "error\n" : "unsure\n");
            }
            else
            {
                // does all gcd equals to 0?
                printf(wtf ? "error\n" : "unsure\n");
            }
        }
        else
        {
            vector<pair<lld,lld>> ans;
            for (auto i : possibleP)
            {
                auto res = solve(i);
                if (res.first != -1)
                    ans.push_back(res);
            }

            if (ans.empty())
            {
                // seems that p exists but no a exists
                printf("error\n");
            }
            else if (ans.size() > 1)
            {
                // many pairs of answers
                printf("unsure\n");
            }
            else
            {
                printf("%lld %lld\n", ans[0].first, ans[0].second);
            }
        }
    }

    return 0;
}

HDU6481 A Math Problem

题目需要求的数字是将 2n 个球放进 n 个盒子中,每个盒子中 2 个球。

考虑盒子的顺序,考虑盒子内球的顺序,就是先将所有的球作一个排列,有 (2n)! 种方法;考虑盒子的顺序,但将每个盒子中两球的位置忽略,那么每种方案会重复 2^k 次;不考虑盒子的顺序,那么每种真实方案还会再重复 n! 次。于是,最后答案就是

F(n) = \frac{(2n)!}{n!2^k} = \prod_{i=1}^n (2i-1)

我们可以构造这样一个多项式:

f_n(x) = \prod_{i=1}^n (2x + 2i-1)

F(n) = f_n(0)。由于取模是 2^{64},那么上述多项式至多只有 64 项。然后我们可以发现,

f_{2n}(x) = f_n(x) f_n(x+n)

于是也可以很方便算出 f_{2n+1}(x),然后就可以愉快的进行倍增啦!

这其中涉及到了一个问题:多项式变量线性变换。记 F(x) = \sum f_i x^i, F(x+n) = G(x) = \sum g_i x^i 有如下关系。

g_i = \sum_{j=i}^\infty C_j^i f_j n^{j-i}

这道题就到这里啦~

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
lld C[64][64], jc[64];

struct poly
{
    lld x[64];

    poly() { for (int i = 0; i < 64; i++) x[i] = 0; }
    poly(lld c1, lld c0) : poly() { x[1] = c1, x[0] = c0; }
    lld operator()(int k) const { return x[k]; }
    lld &operator()(int k) { return x[k]; }

    poly operator*(const poly &b) const
    {
        poly c;
        for (int i = 0; i < 64; i++)
            for (int j = 0; j <= i; j++)
                c(i) += x[j] * b(i-j);
        return c;
    }

    poly moveto(lld n)
    {
        poly ans;
        for (int i = jc[0] = 1; i < 64; i++)
            jc[i] = jc[i-1] * n;
        for (int i = 0; i < 64; i++)
            for (int j = i; j < 64; j++)
                ans(i) += x[j] * jc[j-i] * C[j][i];
        return ans;
    }
};

map<lld, poly> cache;

poly solve(lld n)
{
    if (cache.count(n)) return cache[n];
    poly half = solve(n / 2);
    half = half * half.moveto(n / 2);
    if (n & 1) half = half * poly(2, n*2-1);
    return cache[n] = half;
}

void init()
{
    poly Px(0, 1);
    cache[0] = Px;

    C[0][0] = C[1][0] = C[1][1] = 1;
    for (int n = 2; n < 64; n++)
    {
        C[n][0] = C[n][n] = 1;
        for (int m = 1; m < n; m++)
            C[n][m] = C[n-1][m-1] + C[n-1][m];
    }
}

int main()
{
    init();
    int T; lld n;
    scanf("%d", &T);

    while (T--)
    {
        scanf("%lld", &n);
        printf("%llu\n", solve(n)(0));
    }

    return 0;
}

2.20 JLU 寒假模拟赛

Problem A – 写文案

与此题相似的还有最近 CodeForces Round #540 (Div, 3) – D2Educational Codeforces Round 60 (Div. 2) – C

二分解法

本题需要求 v 的最小值,并且要求 f(v) \ge n

f(v) = \sum_{i=0}^\infty \left\lfloor\frac{v}{k^i}\right\rfloor

考虑这个函数的单调性,非常容易发现,f(v) 随着 v 的增大而增大。直接二分第一个 f(v) \ge n时的 v 即可。可以尝试用这种思路去做上面两题。

#include <bits/stdc++.h>
typedef long long lld;
using namespace std;
int n, k;

int check(int i)
{
    lld sa = 0, sb = i;
    while (sb)
        sa += sb, sb /= k;
    return sa >= n ? 1 : -1;
}

int solve()
{
    if (check(1)==1) return 1;
    int l = 1, r = 1e9+7, m;

    while (r-l>1)
    {
        m = (l+r)>>1;
        if (check(m)==1)
            r = m;
        else
            l = m;
    }

    return r;
}

int main()
{
    scanf("%d %d", &n, &k);
    printf("%d\n", solve());
    return 0;
}

小范围暴力做法

可以发现,答案一定大于等于 n – \lfloor n / k\rfloor – 1,剩余答案大约在 O(\log n) 的范围内,暴力判断即可。但是此做法似乎无法通过提到的两题。

Problem B – 摆石子

动态规划解法

dp[n][m] 表示加入前 n 种颜色后,长度为 m 的排列有多少种方案。那么有

dp[i][j] = \sum_{k=0}^{j} C_j^k dp[i-1][j-k]

其中 1 \le i \le n, 0 \le j \le \sum_{p=1}^ia_p,如果预处理了组合数,那么总的时间复杂度为 O(\prod_{1 \le i \lt j \le n}a_i a_j) \le 10^8,并且是多组测试样例,15秒,是可以做的。

生成函数解法

由于所求的是排列型生成函数,所以考虑排列型生成函数。对于第 i 种颜色,对于长度为 k 的排列,方案数总为 1,其生成函数为 C_i(x),总生成函数为 B(x),即有。

B(x) = \prod_{i=1}^n C_i(x) = \prod_{i=1}^n \sum_{k=0}^{a_i} \frac{x^k}{k!} = \sum_{k=0}^m b_k \frac{x^k}{k!}

那么最后答案就是 \sum_{i=1}^k b_k。如果采用直接暴力多项式乘法,复杂度同上,但是常数好像小一点,维护起来也比较简单;另外如果遇到这种1e9+7的取模多项式乘法,可以考虑用拆系数FFT来做。

#include <vector>
#include <cstdio>
#include <algorithm>
typedef long long lld;
using namespace std;
const int MOD = 1e9+7;
int a[110];
const int MAXN = 1e4+5;
int inv[MAXN], invs[MAXN], fac[MAXN];

void init()
{
    inv[0] = invs[0] = fac[0] = 1;
    inv[1] = invs[1] = fac[1] = 1;

    for (int i = 2; i < MAXN; i++)
    {
        inv[i] = 1LL*(MOD-MOD/i)*inv[MOD%i]%MOD;
        invs[i] = 1LL*invs[i-1]*inv[i]%MOD;
        fac[i] = 1LL*fac[i-1]*i%MOD;
    }
}

vector<int> conv(vector<int> &a, vector<int> &b)
{
    vector<int> ans(a.size()+b.size()-1, 0);

    for (int i = 0; i < a.size(); i++)
    {
        for (int j = 0; j < b.size(); j++)
        {
            ans[i+j] = int((ans[i+j] + 1LL * a[i] * b[j]) % MOD);
        }
    }

    return ans;
}

int main()
{
    int n, cas = 0; init();

    while (scanf("%d", &n) == 1)
    {
        for (int i = 0; i < n; i++)
            scanf("%d", &a[i]);
        sort(a, a+n);
        vector<int> fx(1,1), gx(1,1);

        for (int i = 0; i < n; i++)
        {
            while (gx.size() <= a[i])
                gx.push_back(invs[gx.size()]);
            vector<int> fgx = conv(fx, gx);
            fx = fgx;
        }

        int ans = 0;
        for (int i = 1; i < fx.size(); i++)
            ans = (ans + 1LL * fx[i] * fac[i]) % MOD;
        printf("Case %d: %d\n", ++cas, ans);
    }

    return 0;
}

Problem C – 算数字

简单筛选出2~10000的素数,保存到一个连续数组中。暴力求解;或者维护一个长度为10000的数组,表示对应标号的方案数,然后枚举连续数组的开头 i 和结尾 j 以计数。

#include <cstdio>
typedef long long lld;
using namespace std;
const int MAXN = 1e4+5;
int pri[MAXN], pcn, qzh[MAXN];
int cnt[MAXN];
bool isnp[MAXN];

void init_prime()
{
    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i]) pri[pcn++] = i;
        for (int j = 0; j < pcn; j++)
        {
            if (i*pri[j] >= MAXN) break;
            isnp[i*pri[j]] = 1;
            if (i%pri[j]==0) break;
        }
    }

    for (int i = 1; i <= pcn; i++)
        qzh[i] = qzh[i-1] + pri[i-1];
    for (int i = 0; i < pcn; i++)
        for (int j = i+1; j <= pcn; j++)
            if (qzh[j]-qzh[i] < MAXN)
                cnt[qzh[j]-qzh[i]]++;
            else break;
}

int main()
{
    init_prime();
    int n;
    while (scanf("%d", &n) == 1 && n)
        printf("%d\n", cnt[n]);
    return 0;
}

Problem D – 复制玩具

可以发现,y-1 就是原始玩具投入机器的次数。那么再简单判断一下克隆玩具的情况即可。

#include <bits/stdc++.h>
typedef long long lld;
using namespace std;

int main()
{
    int x, y;
    scanf("%d %d", &x, &y);
    bool checked = true;
    if (y<1) checked = false;
    if ((x-y+1)%2) checked = false;
    if (x < y-1) checked = false;
    if (y == 1 && x != 0) checked = false;
    printf(checked ? "Yes\n" : "No\n");
    return 0;
}

Problem E – 蛙跳

k 的奇偶,计算左跳右跳次数即可。

#include <cstdio>
typedef long long lld;
using namespace std;

int main()
{
    int t; lld a, b, k, ans;
    scanf("%d", &t);

    while (t--)
    {
        scanf("%lld %lld %lld", &a, &b, &k);
        ans = (a-b)*(k/2) + a*(k%2);
        printf("%lld\n", ans);
    }

    return 0;
}

Problem F – 地铁

这是一个图论题啦,感兴趣的同学可以去学一下图论中的最短路和最小生成树。

两点之间线段最短,然后再把同一条地铁线上地铁站之间的距离除以4加入图中,跑一遍最短路,最后将路程的米转换为分钟。

由于这是非负权完全图,直接使用最原始的 dijkstra 算法,复杂度为 O(V^2)。输入的时候要注意技巧,利用scanf的返回值。

#include <map>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
typedef long long lld;
using namespace std;
const int MAXV = 220;
double c[MAXV][MAXV];
double dis[MAXV];
bool vis[MAXV];
int n;

struct point
{
    int x, y;
    point(int _x = 0, int _y = 0) { x = _x, y = _y; }
    bool operator==(const point &b) const { return x==b.x && y==b.y; }
    bool operator<(const point &b) const { return (x==b.x && y<b.y) || (x<b.x); }

    double operator-(const point &b)
    {
        int xx = x - b.x;
        int yy = y - b.y;
        return sqrt(1.0*xx*xx + 1.0*yy*yy);
    }
};

map<point,int> V;
map<int,point> V2;

int get_id(const point &wtf)
{
    if (V.count(wtf)) return V[wtf];
    V[wtf] = ++n; V2[n] = wtf; return n;
}

bool input_station()
{
    int x1, y1, x2, y2;
    if (scanf("%d %d", &x1, &y1) != 2) return false;

    while (scanf("%d %d", &x2, &y2) == 2 && ~x2 && ~y2)
    {
        point p1(x1,y1), p2(x2,y2);
        int i1 = get_id(p1), i2 = get_id(p2);
        double dis = min((p1-p2)/4, c[i1][i2]);
        c[i1][i2] = c[i2][i1] = dis;
        x1 = x2, y1 = y2;
    }

    return true;
}

void dijkstra()
{
    dis[1] = 0;

    for (int i = 1; i <= n; i++)
    {
        double minc = c[0][0];
        int p = -1;
        for (int j = 1; j <= n; j++)
            if (!vis[j] && minc > dis[j])
                minc = dis[j], p = j;
        vis[p] = true;
        for (int j = 1; j <= n; j++)
            if (!vis[j] && dis[j] > dis[p] + c[p][j])
                dis[j] = dis[p] + c[p][j];
    }
}

int main()
{
    int xh, yh, xs, ys;
    scanf("%d %d %d %d", &xh, &yh, &xs, &ys);
    memset(c, 0x55, sizeof(c));
    memset(dis, 0x55, sizeof(dis)); // double INF in memory
    point h(xh,yh), s(xs,ys);
    V[h] = ++n; V2[n] = h;
    V[s] = ++n; V2[n] = s;
    while (input_station());

    for (int i = 1; i <= n; i++)
    {
        for (int j = i+1; j <= n; j++)
        {
            double dis = min(c[i][j], V2[i]-V2[j]);
            c[i][j] = c[j][i] = dis;
        }
    }

    dijkstra();
    printf("%.0f\n", dis[2]*6/1000+1e-8);
    return 0;
}

Problem G – 半素数

先筛选出1-3000的所有素数,保存到一个连续数组中,然后对1-3000的每个数字,将每个素数进行试除,如果满足素因子个数仅为2的条件,那么记录下来,最后统计一下,结束。

#include <bits/stdc++.h>
typedef long long lld;
using namespace std;
const int MAXN = 1e4+5;
int pri[MAXN], pcn, qzh[MAXN];
int cnts[MAXN];
bool isnp[MAXN];

void init_prime()
{
    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i]) pri[pcn++] = i;
        for (int j = 0; j < pcn; j++)
        {
            if (i*pri[j] >= MAXN) break;
            isnp[i*pri[j]] = 1;
            if (i%pri[j]==0) break;
        }
    }

    for (int i = 1; i <= 3000; i++)
    {
        int cnt = 0;
        for (int j = 0; j < pcn; j++)
            if (i % pri[j] == 0)
                cnt++;
        if (cnt == 2)
             cnts[i]++;
    }

    for (int i = 1; i <= 3000; i++)
        cnts[i] = cnts[i]+cnts[i-1];
}


int main()
{
    init_prime();
    int t;
    scanf("%d", &t);
    printf("%d\n", cnts[t]);
    return 0;
}

Problem H – 骨牌摆放

一个长度为 n 的排列,可以由长度为 n-1 的排列后面添加一个 |,或者由长度为 n-2 的排列后面添加 = 得到。那么 F_n = 2F_{n-2} + F_{n-1},递推计算就好了。

由于没有取模,而且 long long 也存不下,所以考虑手动写个十进制加法算法或者使用 JAVA 的大数类即可。

import java.math.BigInteger;
import java.util.Scanner;

public class Main {

    public static void main(String[] args) {
        BigInteger[] f = new BigInteger[251];
        f[0] = f[1] = BigInteger.ONE;
        for (int i = 2; i < 251; i++)
            f[i] = f[i-1].add(f[i-2]).add(f[i-2]);
        Scanner cin = new Scanner(System.in);
        while (cin.hasNextInt())
            System.out.println(f[cin.nextInt()]);
    }

}

Problem I – 几点了

根据每个字符自己的特点,进行超级大的分类讨论好了。(大雾)

#include <cstdio>
typedef long long lld;
using namespace std;
char screen[22][22];

int check(int x)
{
    if (screen[1][x+3] == '.' && screen[4][x] == '.')
        return 5;
    if (screen[1][x+3] == '.' && screen[4][x] == 'X')
        return 6;
    if (screen[0][x+1] == '.' && screen[3][x+1] == '.')
        return 1;
    if (screen[0][x+1] == '.' && screen[3][x+1] == 'X')
        return 4;
    if (screen[0][x+1] == 'X' && screen[3][x+1] == '.' && screen[6][x+1] == '.')
        return 7;
    if (screen[3][x+1] == '.')
        return 0;
    if (screen[4][x+3] == '.')
        return 2;
    if (screen[1][x] == '.')
        return 3;
    if (screen[4][x] == '.')
        return 9;
    return 8;
}

void solve()
{
    for (int i = 0; i < 7; i++)
        scanf("%s", screen[i]);
    printf("%d", check(0));
    printf("%d", check(5));
    printf(":");
    printf("%d", check(12));
    printf("%d", check(17));
    printf("\n");
}

int main()
{
    int t; scanf("%d", &t);
    while (t--) solve();
    return 0;
}

Problem J – 匹配字

这题是询问,有多个模板字符串,他们中有多少个在目标字符串中出现。

这个算法叫做 Aho-Corasick automaton,中文名叫做 AC自动机,推荐阅读 这篇文章 来了解它的算法内容。

Sparse Table

稀疏表,又称ST表,预处理 O(n \log n),查询 O(1)

st[i][j] 表示在数组中,[i,i+2^j-1] 区间的最值。

当我们需要查询 [L,R] 的最值时,由于 R-L+1 可能不是正好 2^j 形式,那么我们可以考虑用某个 j,查询 [L,L+2^j-1][R-2^j+1,R],并使这两个区间交叠。

int a[1010];
int st[1010][12];

void init(int n)
{
    for (int i = 0; i < n; i++)
        st[i][0] = a[i];
    for (int j = 1; (1<<j) <= n; j++)
        for (int i = 0; i+(1<<j)-1 < n; i++)
            st[i][j] = min(st[i][j-1], st[i+(1<<(j-1))][j-1]);
}

int search(int l, int r)
{
    int k = int(log2(r-l+1.0));
    return min(st[l][k], st[r-(1<<k)+1][k]);
}

并且由代码可以看出,它还可以处理区间最大值/最值坐标等。