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]);
}

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

ZOJ4069 Sub-cycle Graph

是时候回顾一波青岛站打铁时没想出来的题了QAQ。

题目传送门:ZOJ4069

求由 n 个点组成 k 条链,形成的所有图的方案数,对 1e9+7 取模。

首先,我们先考虑由 m 个点形成 1 条链的不同方案个数。显然对任意排列,去除倒过来的那个就行了。将方案数记为 l_m,则可以得到

l_m = \begin{cases} 1, & m=1 \\ \frac{m!}{2}, & m>1 \end{cases}

那么似乎接下来就是求它的生成函数了。在这里我们选择它的指数型生成函数。

L(x) = \sum_{m=1}^\infty l_m \frac{x^m}{m!} = \frac{1}{2}\left(2x+x^2+x^3+\cdots\right) = x\left(1-\frac{x}{2}\right)\frac{1}{1-x}

那么现在取 k 条链,我们可以计算出

A_k(x) = \sum_{n=1}^\infty a_{n,k} \frac{x^n}{n!} = \left[L(x)\right]^k = x^k \left(1-\frac{x}{2}\right)^k \left(\frac{1}{1-x}\right)^k

由于这 k 条链之间是不考虑先后关系的,所以

b_{n,k} = \frac{a_{n,k}}{k!}

就是最终的答案。

考虑一下题目是多测的,而且测试样例数高达 T = 2 \times 10^4,是不能进行FFT的。

我们可以发现

\left(1-\frac{x}{2}\right)^k = \sum_{s=0}^k C_k^s \left(-\frac{1}{2}\right)^s x^s

是可以用二项式定理做出来的。

然后考虑一下,后面的 1/(1-x)^k 与前面这项进行卷积,其实就相当于对前面的数列进行了 k 次前缀和。

emmmmm,看到网上有个dalao说

求k次前缀和的第m项,想必现在是个人都应该会了

突然自闭.jpg。然后推了推。假设这个数列原来是

\lbrace x_0, x_1, x_2, \cdots \rbrace

对它进行 k+1 次前缀和,它会变成

\lbrace C_k^k x_0, C_{k+1}^k x_0 + C_k^k x_1, C_{k+2}^k x_0 + C_{k+1}^k x_1 + C_k^k x_2, \cdots \rbrace

即第 n 项为

\sum_{i=0}^n C_{k+i}^k x_{n-i}

单次询问复杂度是 O(m) 的,那么现在可以开始写代码啦!

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MAXN = 1e5+5;
const int MOD = 1e9+7;
lld fac[MAXN], inv[MAXN], invs[MAXN];

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

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

inline lld C(int n, int k)
{
    if (n < k) return 0;
    return fac[n] * invs[k] % MOD * invs[n-k] % MOD;
}

void solve()
{
    int n, m, k;
    scanf("%d %d", &n, &m);

    if (m > n)
    {
        printf("0\n");
    }
    else if (m == n)
    {
        printf("%lld\n", fac[n-1]*inv[2]%MOD);
    }
    else
    {
        k = n - m;
        lld ans = 0, inv2s = 1;

        for (int i = m; i >= 0; i--)
        {
            ans += C(k-1+i,i) * C(k, m-i) % MOD * inv2s % MOD;
            inv2s = inv2s * (MOD-inv[2]) % MOD;
        }

        ans = (ans % MOD + MOD) % MOD;
        ans = ans * fac[n] % MOD * invs[k] % MOD;
        printf("%lld\n", ans);
    }
}

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

Camp Day7 Div2

Problem G – 抢红包机器人

一个简单的暴力做法。

枚举第 i 个人为机器人,则根据聊天记录查找所有的机器人,例如floyd-warshall算法或者直接拿bitset维护一下,就可以统计出,在他为机器人的时候,有多少人是机器人。

Problem C – 斐波那契数列

首先可以发现,x \& (x-1)=x-lowbit(x),那么

ans = \sum_{i=1}^R F_i-\sum_{i=1}^Rlowbit(F_i)

首先斐波那契数列的求和可以通过矩阵快速幂计算。

\begin{aligned} & F_n = F_{n-1} + F_{n-2} \\ \Leftrightarrow\ & S_n = 2S_{n-1} – S_{n-3} \end{aligned}

\left[\begin{matrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 2 & 0 & -1\end{matrix}\right]^{n-2}\left[\begin{matrix} S_{0} \\ S_{1} \\ S_{2} \end{matrix}\right] = \left[\begin{matrix} S_{n-1} \\ S_{n} \\ S_{n+1} \end{matrix}\right]

关于 F_n 的lowbit的计算。首先我们可以打表注意到,当 n\not=3k 时,F_n 为奇数。那么考虑 n=3k 时,考虑 F_n 的2的幂次。

\begin{aligned} F_n & = \frac{1}{\sqrt5}\left((\frac{1+\sqrt5}{2})^{3k}+(\frac{1-\sqrt5}{2})^{3k}\right) \\ & = \frac{1}{\sqrt5}\left((2+\sqrt5)^k+(2-\sqrt5)^k\right) \\ & = \frac{1}{\sqrt5}\left(\sum_{i=1}^k C_k^i 2^i \left(\sqrt5^{k-i}-(-\sqrt5)^{k-i}\right)\right) \\ & = \sum_{i=1}^k C_k^i 2^i \left(\sqrt5^{k-i-1}+(-\sqrt5)^{k-i-1}\right) \end{aligned}

k 为奇数,i 为奇数时,显然和式项对2的幂次无贡献。当 k 为奇数,i 为偶数的时候,考虑 i=1 即可得到 F_n 中2的幂次为 1。当 k 为偶数时,考虑 i=2,可以得到 lowbit(4k)|F_n

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int mod = 998244353;
lld SS[64] = { 0, 2 };

struct matrix
{
    lld V[3][3];
    matrix() noexcept { memset(V, 0, sizeof(V)); }
    lld& operator()(int a, int b) { return V[a][b]; }

    matrix& operator*=(const matrix b)
    {
        static lld tmp[3][3];
        memset(tmp, 0, sizeof(tmp));

        for (int i = 0; i < 3; i++)
        for (int j = 0; j < 3; j++)
        {
            for (int k = 0; k < 3; k++)
                tmp[i][j] += V[i][k] * b.V[k][j] % mod;
            tmp[i][j] %= mod;
        }

        memcpy(V, tmp, sizeof(tmp));
        return *this;
    }
};

lld SumFib(lld n)
{
    if (n < 3) return n;
    matrix A, B;
    A(0,0) = A(1,1) = A(2,2) = 1;
    B(0,1) = B(1,2) = 1;
    B(2,0) = -1; B(2,2) = 2;

    n -= 2;
    while (n)
    {
        if (n & 1) A *= B;
        B *= B;
        n >>= 1;
    }

    return (A(2,1) + 2 * A(2,2)) % mod;
}

lld Sum001002(lld R)
{
    lld r6 = R / 6;
    lld ans = (r6 % mod) * 6 % mod;
    lld m6 = R % 6;
    ans += m6;
    if (m6 > 2) ans += 1;
    return ans % mod;
}

int log2(lld ans)
{
    int a = 0;
    while (ans)
        ans >>= 1, a++;
    return a;
}

lld Sum1213(lld R)
{
    if (R == 0) return 0;
    lld R2 = R;
    while ((R2&(-R2)) != R2)
        R2 ^= R2&(-R2);
    return (SS[log2(R2)] + Sum1213(R-R2)) % mod;
}

void solve()
{
    lld R; scanf("%lld", &R);
    lld s1 = SumFib(R);
    lld s2 = Sum001002(R);
    lld s3 = Sum1213(R/6)*4%mod;
    printf("%lld\n", ((s1-s2-s3) % mod + mod) % mod);
}

int main()
{
    for (int i = 2; i < 64; i++)
        SS[i] = (SS[i-1]*2 + (1LL<<(i-1))) % mod;
    int T; scanf("%d", &T);
    while (T--) solve();
    return 0;
}

Problem D – 二次函数

首先可以注意到,将 F(x)=x^2+ax+b 进行水平移动 F(x-\lfloor a/2\rfloor),不改变对应相同处的函数值。在平移后,所有的二次函数都可以被规约为 A(x)=x^2+C_AB(x)=x^2+x+C_B。那么,进行一下分类讨论好了。

  • 当同时存在有 A(x)B(x) 时,令 A(x)=B(x),即有 x=C_A – C_B 处两者相交。

  • 当三者均为 A(x) 型时

    • C_1C_2 不同奇偶性,令 A_1(x+1)=A_2(x),则有 2x+1+C_1-C_2=0
    • C_1C_2 同奇偶性且 C_1 \equiv C_2 \mod 4,令 A_1(x+2)=A_2(x),则有 4x+4+C_1-C_2=0
    • 可以证明,如果三个函数各不相同,根据抽屉原理,必定存在上述两种情况之一
  • 当三者均为 B(x) 型时
    • C_1C_2 同奇偶性,令 B_1(x+1)=B_2(x),则有 2x+2+C_1-C_2=0
    • 可以证明,如果三个函数各不相同,根据抽屉原理,必定是有两个函数同奇偶性

可以证明,所有函数一定可以规约到以上规范形式。

Problem J – 强壮的排列

和前几天的置置置换是很像的,都是 p[2n]>max\lbrace p[2n-1], p[2n+1]\rbrace,但是敦哥哥给了个更适合这题的做法。记 f_n 为长度为 n 的排列符合题目要求的概率,则 n!f_n 就是本题的答案。

考虑到在此序列中,n 为最大数字在 i 的位置,那么我们可以将此序列一切两半,得到这样的一个递推公式。

f_{n} = \frac{1}{n}\sum_{i=1}^{n}f_{n-i}f_{i-1} = \frac{1}{n} \sum_{i+j=n-1} f_i f_j

这个式子可以考虑使用分治FFT来做,复杂度为 O(n \log^2 n)

进一步考虑,令 F(x)=\sum_{n=0}^\infty f_n x^n,那么有

\begin{aligned} nf_n &= \sum_{i=1}^{n}f_{n-i}f_{i-1} \\ \Leftrightarrow\ nf_nx^{n-1} &= \sum_{i+j=n-1}f_i x^i \times f_j x^j \\ \Leftrightarrow\ F'(x) &= F(x)^2 + f_1 \end{aligned}

其中 f_0=1,根据观察(大雾)可以发现 F(x) = \tan x = \frac{\sin x}{\cos x} 是这个微分方程的一个解。

利用麦克劳林展开,得到 \sin x\cos x 的组合型生成函数表示,然后使用NTT来进行多项式求逆和乘法,即可得到我们的答案。预处理阶乘和阶乘逆元 O(n),生成两者多项式 O(n),处理逆元 O(n \log n),多项式乘法 O(n\log n),所以总复杂度是 O(n \log n)

#include <bits/stdc++.h>
typedef long long lld, num_t;
using namespace std;
const int mod = 998244353, G=3;
const size_t max_len = 19;
const int MAXN = 524289;
int fac[MAXN], inv[MAXN], invs[MAXN];

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

inline num_t fpow(lld a, lld k)
{
    lld ans = 1;
    while (k)
    {
        if (k & 1) ans = ans * a % mod;
        a = a * a % mod;
        k >>= 1;
    }
    return ans;
}

class Fourier
{
public:
    int len, N;

private:
    num_t wmk[1 << max_len];
    int rev[1 << max_len];

    inline num_t create(int m)
    {
        int k = (mod-1)/m;
        if (k < 0) k += mod-1;
        return fpow(G, k);
    }

    void dft(num_t a[], int DFT)
    {
        for (int i = 0; i < N; i++)
            if (i < rev[i])
                swap(a[i], a[rev[i]]);

        for (int m = 2; m <= N; m <<= 1)
        {
            int m2 = m >> 1;
            num_t wm = create(DFT * m);
            wmk[0] = 1;
            for (int j = 1; j < m2; j++)
                wmk[j] = wmk[j-1] * wm % mod;
            num_t t, u;

            for (int k = 0; k < N; k += m)
            {
                for (int j = 0; j < m2; j++)
                {
                    t = wmk[j] * a[k+j+m2] % mod;
                    u = a[k+j];
                    a[k+j] = (u + t) % mod;
                    a[k+j+m2] = ((u - t) % mod + mod) % mod;
                }
            }
        }

        if (DFT == -1)
            for (int i = 0; i < N; i++)
                a[i] = (a[i] * inv[N] % mod + mod) % mod;
    }

public:
    void init(int _len)
    {
        len = _len, N = 1 << _len;
        for (int i = 0; i < N; i++)
            rev[i] = (rev[i>>1]>>1)|((i&1)<<len-1);
    }

    void DFT(num_t a[]) { dft(a, 1); }
    void IDFT(num_t a[]) { dft(a, -1); }

    int log2_ceil2(int len)
    {
        int ans = 0;
        while (len) len >>= 1, ans++;
        return ans;
    }
} trans;

num_t x[MAXN], y[MAXN], dp[MAXN];

void cdqFFT(int L, int R)
{
    if (L == R)
    {
        if (L == 0) dp[L] = 0;
        else if (L == 1) dp[L] = 1;
        else dp[L] = (dp[L] * inv[L] % mod + mod) % mod;
        return;
    }

    int mid = (L + R) >> 1;
    cdqFFT(L, mid);

    int len = R - L + 1;
    trans.init(trans.log2_ceil2(len));

    // x = dp[0], dp[1], dp[2], ..., dp[len], 0, 0, ...
    // y = dp[L], dp[L+1], dp[L+2], ..., dp[mid], 0, 0, ...
    for (int i = 0; i < trans.N; i++)
    {
        x[i] = i < len ? dp[i] : 0;
        y[i] = (i <= mid-L) ? dp[i+L] : 0;
    }

    trans.DFT(x); trans.DFT(y);
    for (int i = 0; i < trans.N; ++i)
        x[i] = x[i] * y[i] % mod;
    trans.IDFT(x);

    for (int i = mid + 1; i <= R; ++i)
        dp[i] = (dp[i] + x[i-L-1] * (L?2:1)) % mod;
    cdqFFT(mid + 1, R);
}

void get_inv(num_t a[], num_t b[], int len)
{
    static num_t A[MAXN], B[MAXN];

    if (len == 1)
    {
        b[0] = fpow(a[0], mod-2);
        return;
    }

    get_inv(a, b, len>>1);
    trans.init(trans.log2_ceil2(len));
    for (int i = 0; i < len; i++)
        A[i] = a[i], B[i] = b[i], A[i+len] = 0, B[i+len] = 0;
    trans.DFT(A); trans.DFT(B);
    for (int i = 0; i < trans.N; i++)
        A[i] = (A[i] * B[i] % mod) * B[i] % mod;
    trans.IDFT(A);

    for (int i = 0; i < len; i++)
        b[i] = (2LL * b[i] % mod - A[i] + mod) % mod;
}

int main()
{
    init_factory();

#ifdef CDQFFT
    // f_n = 1/n * \sum {i+j=n-1} f_i*f_j
    cdqFFT(0, 131071);
#else
    // y for cosx then x for secx
    // then y for sinx then dp for secx+tanx
    for (int i = 0; i < 262144; i++)
        y[i] = (i&1) ? 0 : (i&2) ? mod-invs[i] : invs[i];
    get_inv(y, x, 262144);
    for (int i = 1; i < 262144; i++)
        y[i] = (i&1) ? (i&2) ? mod-invs[i] : invs[i] : 0;
    trans.init(19);
    trans.DFT(x); trans.DFT(y);
    for (int i = 0; i < trans.N; i++)
        dp[i] = x[i] * y[i] % mod;
    trans.IDFT(dp);
#endif

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

    while (T--)
    {
        scanf("%d", &n);
        printf("%lld\n", dp[n] * fac[n] % mod);
    }

    return 0;
}

Camp Day5 Div2

Problem A – Gactus Draw

由于有要求边不能有公共点,那么我们可以让每颗子树都相离较远。如果我们将树的 dfs序 作为横坐标,depth 作为纵坐标,那么构造出来的图形一定是符合要求的。

Problem C – Division

题中的操作相当于将数字进行一次右移,相当于减去了某个某个指定的数值。

如果我们将会被减去的数记录下来,至多 32n 个,那么先把和保存下来,然后减去的数保存下来排个序,把最前面至多 k 个数字加进和即可。

Problem F – Kropki

考虑做一个状压dp,记为 dp[S][n],表示以状态 S 排布的数字以 n 为结尾的方案数,其中 S 中已选上的元素数量为 c(S)。另外点的01串存在 stat[] 中。那么有

dp[S][n] = \begin{cases} \sum_{u=2n,n/2} dp[S-\lbrace n \rbrace][u], & stat[c(S)]=1 \\ \sum_{u\not= 2n,n/2}dp[S-\lbrace n\rbrace][u], & stat[c(S)]=0 \end{cases}

然后计算即可。可以写成搜索的形式,用起bitset,然后记忆化,这样写法似乎会舒服一点(大雾)。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int LEN = 1 << 15;
const int mod = 1e9+7;
lld dp[LEN][16], n, dpt;
char st[20];

lld dfs(int S, int end)
{
    bitset<15> wtf(S);
    size_t len = wtf.count();
    if (!wtf[end-1]) return 0;
    if (len == 1) return 1;
    if (dp[S][end] != -1) return dp[S][end];

    lld ans = 0;
    wtf.flip(end-1);
    dpt++;
    for (int i = 1; i <= n; i++)
    {
        if (!wtf[i-1]) continue;

        if (end == i * 2)
        {
            if (st[len-2] == '1')
                ans += dfs(wtf.to_ulong(), i);
        }
        else if (i == end * 2)
        {
            if (st[len-2] == '1')
                ans += dfs(wtf.to_ulong(), i);
        }
        else
        {
            if (st[len-2] == '0')
                ans += dfs(wtf.to_ulong(), i);
        }
    }

    return dp[S][end] = ans % mod;
}

int main()
{
    memset(dp, -1, sizeof(dp));
    scanf("%lld %s", &n, st);
    lld ans = 0;
    for (int i = 1; i <= n; i++)
        ans = (ans + dfs((1<<n)-1, i)) % mod;
    printf("%lld\n", ans);
    return 0;
}

很迷啊,不会写最原始的状压dp,所以只能写写这种奇怪的东西了。

Problem H – Nested Tree

考虑一棵树的每对点间距离和的一个性质:对于某个边 e(u,v),在距离和中的贡献是 u 侧的节点数乘以 v 侧节点数。那么,我们可以利用一次dfs,对以某节点为根的子树进行搜索,统计出子树中的节点个数 x,则边另一侧的节点个数为 nm-x,所以答案就可以很轻松的维护出来。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int mod = 1e9+7;
vector<int> G[1000010];
bool vis[1000010];
int n, m;
long long ans;

int dfs(int u)
{
    if (vis[u]) return -1;
    vis[u] = true;

    int cnt = 1;
    for (auto v : G[u])
    {
        int x = dfs(v);
        if (x == -1) continue;
        ans += 1LL * x * (n*m - x) % mod;
        cnt += x;
    }

    return cnt;
}

int main()
{
    int a, b, u, v;
    scanf("%d %d", &n, &m);

    for (int i = 1; i < n; i++)
    {
        scanf("%d %d", &u, &v);

        for (int j = 0; j < m; j++)
        {
            G[j*n+u].push_back(j*n+v);
            G[j*n+v].push_back(j*n+u);
        }
    }

    for (int i = 1; i < m; i++)
    {
        scanf("%d %d %d %d", &a, &b, &u, &v);
        G[(a-1)*n+u].push_back((b-1)*n+v);
        G[(b-1)*n+v].push_back((a-1)*n+u);
    }

    assert(dfs(1) == n*m);
    printf("%lld\n", (ans % mod + mod) % mod);
    return 0;
}

Problem J – Special Judge

线段相交的判断。