DoorKickers的博客

博客

[AH2017/HNOI2017] 礼物

2020-04-19 12:22:55 By DoorKickers

首先, 两个手环增加非负整数亮度, 等价于其中一个增加一个整数亮度.

答案应该是 ni=1(xiyi+c)2 展开 ni=1(x2i+y2i2xiyi+2xic2yic+c2)=ni=1(x2i+y2i+c2)+2cni=1(xiyi)2ni=1xiyi 然后发现我们只需要算最后一项. ni=1xiyix0=y0=0

考虑将x翻转, 然后左移一位 原式等于 ni=0xniyi 是个卷积.

考虑将x 破环成链(主要是卷积的性质)

假设

1 2 3 1 2 3 0 :x

0 4 5 6 0 0 0 :y

设结果为z(x) zi=ik=0xikyk 发现卷起来以后n ~ 2n - 1取个max就是答案. 代码:

#include <bits/stdc++.h>

using namespace std;

#define int long long

const int g = 3;
const int N = 1e6 + 10;
const int mod = 1004535809;

int omega[N], omegaInv[N], a[N], b[N], rev[N];

inline int fp(int a, int n, int mod)
{
    int res = 1;
    while (n)
    {
        if (n & 1) res = (res * a) % mod;
        n >>= 1;
        a = (a * a) % mod;
    }
    return res;
}

inline int inv(int x)
{
    return fp(x, mod - 2, mod);
}

inline void init(const int n)
{
    int k = 0;
    while ((1 << k) < n) k++;
    int x = fp(g, (mod - 1) / n, mod);
    omega[0] = omegaInv[0] = 1;
    for (int i = 1; i < n; i++)
    {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << k - 1);
        omega[i] = (omega[i - 1] * x) % mod;
        omegaInv[i] = inv(omega[i]);
    }
}

inline void ntt(int *a, const int n, int *omega)
{
    for (int i = 0; i < n; i++)
    {
        if (i > rev[i]) continue;
        swap(a[i], a[rev[i]]);
    }
    for (int l = 2; l <= n; l *= 2)
    {
        int m = l / 2;
        for (int *p = a; p != a + n; p += l)
        {
            for (int i = 0; i < m; i++)
            {
                int t = (omega[n / l * i] * p[m + i]) % mod;
                p[m + i] = (p[i] - t + mod) % mod;
                p[i] = (p[i] + t) % mod;
            }
        }
    }
}

inline void dft(int *a, const int n)
{
    ntt(a, n, omega);
}

inline void idft(int *a, const int n)
{
    ntt(a, n, omegaInv);
    int x = inv(n);
    for (int i = 0; i < n; i++)
        a[i] = (a[i] * x) % mod;
}

inline void multiply(int *a1, const int n1, int *a2, const int n2)
{
    int n = 1;
    while (n < n1 + n2 + 2) n *= 2;
    init(n);
    dft(a1, n);
    dft(a2, n);
    for (int i = 0; i < n; i++)
        a1[i] = (a1[i] * a2[i]) % mod;
    idft(a1, n);
    idft(a2, n);
}

inline void squared(int *a, int n1)
{
    int n = 1;
    while (n < n1 * 2 + 2) n *= 2;
    init(n);
    dft(a, n);
    for (int i = 0; i < n; i++)
        a[i] = (a[i] * a[i]) % mod;
    idft(a, n);
}

inline void poly_print(int *a, int n)
{
    for (int i = 0; i <= n; i++)
        cout << a[i] << ' ';
    cout << endl;
}

inline void test()
{
    double m;
    cin >> m;
    int t = m;
    cout << "fk " << t << endl;
    int n = floor(m);
    cout << n << endl;
    // n = ceil(m);
    cout << n << endl;
    exit(0);
}


signed main()
{
    // freopen("input.txt", "r", stdin);
    // test();
    int n, m;
    cin >> n >> m;
    int X = 0;
    int Y = 0;
    int XX = 0;
    int YY = 0;
    for (int i = 0; i < n; i++)
    {
        cin >> a[i];
        X += a[i];
        XX += fp(a[i], 2, mod);
        a[i + n] = a[i];
    }
    reverse(a, a + n * 2);
    // for (int i = 0; i < n * 2; i++)
    // {
    //     cout << a[i] << ' ';
    // }
    // cout << endl;
    for (int i = 1; i <= n; i++)
    {
        cin >> b[i];
        Y += b[i];
        YY += fp(b[i], 2, mod);
    }
    int ans1 = ((Y - X) / n) + 1;
    // cout << "fk " << (double)((Y - X) / n);
    int ans2 = (Y - X) / n;
    int ans3 = (Y - X) / n - 1;
    // // cout << endl;
    // cout << "ans1 " << ans1 << endl;
    // cout << "ans2 " << ans2 << endl;
    // cout << "ans3 " << ans3 << endl;
    // cout << "x " << (Y - X) << endl;
    // cout << "n " << n << endl;
    // cout << fixed << setprecision(50) << (double)((Y - X) / n) << endl;
    // double t1 = Y - X;
    // double t2 = n;
    // cout << t1 / t2 << endl;
    // cout << "double " << (double)((Y - X) / n) << endl;
    // cout << (int)(ceil((double)((Y - X) / n))) << endl;
    int res1 = XX + YY;
    int res2 = min(ans1 * 2 * (X - Y) + n * ans1 * ans1, ans2 * 2 * (X - Y) + n * ans2 * ans2);
    res2 = min(res2, ans3 * 2 * (X - Y) + n * ans3 * ans3);
    int len = n * 2 - 1;
    multiply(a, len, b, len);
    int res3 = 0;
    for (int i = n; i < 2 * n; i++)
        res3 = max(res3, a[i]);
    cout << res1 + res2 - 2 * res3 << endl;
    return 0;
}

[SDOI2017]序列计数

2020-04-19 08:31:39 By DoorKickers

而且, 这题用NTT + FFT貌似不可行.

NTT模数不对, FFT爆精度.

首先想一个暴力算法:

f[i][j] 表示选i个, 在模p 意义下是j的方案数. f[i][j]=f[i1][k]×l(k+l)mod 复杂度O(nm)

然后发现一个优化, 可以只记录一下每个数模p的结果, 记录一下cnt

然后可以优化到差不多O(np)

然后又发现一个性质 f[k] * f[j] = f[j + k] 然后类比 x^a * x^b = x^{a + b} 如何快速计算x^{a + b}

快速幂就好了.

同理, 计算f 也可以快速幂.

然后数组可以用滚动数组优化.

最终复杂度是O(\log n * p) 再算上一些奇奇怪怪的预处理.

#include <bits/stdc++.h>

using namespace std;

// #define int long long
#define ll long long

const int N = 2e2 + 10;
const int L = 2e7 + 10;
const int mod = 20170408;

ll F[N], G[N], f[N], g[N], n, m, p;

int v[L], q[L / 10], cnt[N];


inline void multiply(ll *a, ll *b)
{
    static ll res[N];
    for (int i = 0; i < p; i++)
        for (int j = 0; j < p; j++)
        {
            res[i + j] = (a[i] * b[j] % mod + res[i + j]) % mod;
        }
    for (int i = 0; i < p; i++)
    {
        a[i] = (res[i] + res[i + p]) % mod;
        res[i] = res[i + p] = 0;
    }
}

signed main()
{
    cin >> n >> m >> p;
    F[1] = G[1] = f[1] = g[1] = 1;
    for (int i = 2; i <= m; i++)
    {
        f[i % p]++;
        g[i % p]++;
        F[i % p]++;
        G[i % p]++;
        if (v[i] == 0)
        {
            v[i] = i;
            // cout << i << ' ';
            G[i % p]--;
            g[i % p]--;
            q[++q[0]] = i;
        }
        for (int j = 1; j <= q[0]; j++)
        {
            if (q[j] > v[i] || q[j] * i > m) break;
            v[i * q[j]] = q[j];
        }
    }
    n--;
    while (n)
    {
        if (n & 1)
        {
            multiply(F, f);
            multiply(G, g);
        }
        n >>= 1;
        multiply(f, f);
        multiply(g, g);
    }
    cout << (F[0] - G[0] + mod) % mod << endl;
    return 0;
}

多项式入门推式子好题

2020-04-18 16:05:01 By DoorKickers

题目大意

给定n个数q_1, q_2 \dots q_n 定义 F_j = \sum_{i = 1}^{j - 1}\frac{q_i \times q_j}{(i - j)^ 2} - \sum_{i = j + 1}^{n}\frac{q_i \times q_j}{(i - j) ^ 2}

同时 E_i = \frac{F_i}{q_i}

你的任务是对于所有的1 \leq i \leq n, 求出 E_i


Solution

推下公式: E_j = \frac{F_j}{q_j} = \sum_{i = 1}^{j - 1}{\frac{q_i}{(i - j)^2}} - \sum_{i = j + 1}^{n}{\frac{q_i}{(i - j) ^2}} 然后如果要用FFT加速的话, 要把这个推成卷积式

卷积式就是C_k = \sum_{0}^{k}{A[i] * B[k - i]}

继续推 E_j = \sum_{i = 1}^{j}{\frac{q_i}{(i - j)^2}} - \sum_{i = j}^{n}{\frac{q_i}{(i - j) ^ 2}} \\ f[i] = q_i, g[i] = \frac{1}{i ^ 2}

原式变成 E_j = \sum_{i = 1}^{j}{f[i] * g[j - i]} - \sum_{i = j}^{n}{f[i] * g[i - j]} 然后设f[0] = g[0] = 0 E_j = \sum_{i = 0}^{j}f[i] * g[j - i] - \sum_{i = j}^{n}f[i] * g[i - j] 发现左边已经是一个卷积式子了, 继续推右边的.

将他展开, 发现 \sum_{i = j}^{n}f[i] * g[i - j] = \\ f[j] * g[0] + f[j + 1] * g[1] + \dots + f[j + (n - j)] * g[n - j] \\ = \sum_{i = 0}^{n - j}f[j + i] * g[i] 然后发现f还不太行, 翻转一下.

f'[i] = f[n - i] \sum_{i = 0}^{n - j}f[j + i] * g[i] = \sum_{i = 0}^{n - j}f'[n - (j + i)] * g[i] \\ = \sum_{i = 0}^{n - j}f'[n - j - i] * g[i] \\ 然后设t = n - j

原式 \sum_{i = 0}^{t}f'[t - i] * g[i] 答案就是 E_j = \sum_{i = 0}^{j}{f[i] * g[j - i]} - \sum_{i = 0}^{t}f[t - i] * g[i]\quad (t = n - j) 然后设A(x) = \sum_{i = 0}^{n}f[i], B(x) = \sum_{i = 0}^{n} g[i] , C(x) = \sum_{i = 0}^{n}f'[i]

D(x) = A(x) * B(x), E(x) = B(x) * C(x)

然后 E_j = D(j) - E(n - j)

就做完了

代码:

#include <bits/stdc++.h>

using namespace std;

#define int long long

const long double PI = 3.14159265358978323846264338;
const int N = 1e6 + 10;

complex<double> omega[N], omegaInv[N];
int rev[N];
double f[N], g[N], ff[N], res1[N], res2[N];

inline void init(const int n)
{
    int k = 0;
    while ((1 << k) < n) k++;
    for (int i = 0; i < n; i++)
    {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));
        omega[i] = complex<double>(cos(2 * PI / n * i), sin(2 * PI / n * i));
        omegaInv[i] = conj(omega[i]);
    }
    rev[0] = 0;
}

void fft(complex<double> *a, const int n, const complex<double> *omega)
{
    for (int i = 0; i < n; i++)
    {
        if (i > rev[i]) continue;
        swap(a[i], a[rev[i]]);
    }
    for (int l = 2; l <= n; l *= 2)
    {
        int m = l / 2;
        for (complex<double> *p = a; p != a + n; p += l)
            for (int i = 0; i < m; i++)
            {
                complex<double> t = omega[n / l * i] * p[m + i];
                p[m + i] = p[i] - t;
                p[i] += t;
            }
    }
}

inline void dft(complex<double> *a, const int n)
{
    fft(a, n, omega);
}

inline void idft(complex<double> *a, const int n)
{
    fft(a, n, omegaInv);
    for (int i = 0; i < n; i++)
    {
        a[i] /= n;
    }
}

inline void squared(int *a1, const int n1)
{
    int n = 1;
    while (n < n1 + n1 + 2) n *= 2;
    static complex<double> c1[N];
    for (int i = 0; i <= n1; i++) c1[i].real(a1[i]);
    init(n);
    dft(c1, n);
    for (int i = 0; i < n; i++)
        c1[i] *= c1[i];
    idft(c1, n);
    for (int i = 0; i <= n1 * 2; i++)
        a1[i] = static_cast<int>(floor(c1[i].real() + 0.5));
}

inline void multiply(double *a1, const int n1, double *a2, const int n2, double *res)
{
    int n = 1;
    while (n < n1 + n2 + 2) n *= 2;
    static complex<double> c1[N], c2[N];
    memset(c1, 0, sizeof(c1));
    memset(c2, 0, sizeof(c2));
    for (int i = 0; i <= n1; i++) c1[i].real(a1[i]);
    for (int i = 0; i <= n2; i++) c2[i].real(a2[i]);
    init(n);
    dft(c1, n);
    dft(c2, n);
    for (int i = 0; i < n; i++)
        c1[i] *= c2[i];
    idft(c1, n);
    // idft(c2, n);
    for (int i = 0; i <= n1 + n2; i++)
        res[i] = c1[i].real();
}

signed main()
{
    int n;
    cin >> n;
    for (int i = 1; i <= n; i++)
    {
        cin >> f[i];
        ff[n - i] = f[i];
        g[i] = (double)(1.0 / i / i);
    }
    // for (int i = 1; i <= n; i++)
    //     cin >> g[i];
    multiply(f, n, g, n, res1);
    // for (int i = 0; i <= n * 2; i++) 
    //     cout << fixed << setprecision(10) << res1[i] << ' ';
    // cout << endl;
    // exit(0);
    multiply(ff, n, g, n, res2);
    for (int i = 1; i <= n; i++)
        cout << fixed << setprecision(3) << res1[i] - res2[n - i] << endl;
    return 0;
}

SAM Bilion 老师的SAM 入门小题题

2020-04-10 16:15:55 By DoorKickers

首先 发现L满足单调性, 可以二分

考虑一个DP​

f_i 表示将前i个元素分成若干段最大有多少被匹配的

转移

f_i = \max(f_j + i - j, f_{i - 1}) \quad (i - len(i) \leq j \leq i - L)

len(i) 表示以i 结尾最长能被匹配到的子串的长度, L是当前二分的值

如何求出以i结尾最长能被匹配的子串长度呢?

一个朴素的方法是对每个位置暴力走, 然后能走多少就是多少, 但复杂度太高

发现走的路程都是重复的, 可以记录一些信息来节省时间.

我们可以从位置1开始暴力匹配, 如果走不下去, 就跳后缀链接再匹配

这样可以充分利用之前的信息, 减少复杂度.

然后观察一下转移方程, 发现状态不能优化, 只能优化转移, 发现一个性质.

i - len(i) \leq i + 1 - len(i + 1) \\i - L < i + 1 - L

发现决策的位置是递增的, 于是我们可以维护f_j - j 的单调递减队列.

总体思路就是对M个串建广义SAM, 然后对每个询问串, 先预处理len(i)

然后二分答案, 用dp check.

总复杂度为O(\sum{|S|} + \sum{|P|\log|P|})

代码:

#include <bits/stdc++.h>

using namespace std;

#define int long long

const int N = 1e6 + 10;

char s[N];
int tot, fa[N * 2], ch[N * 2][26], len[N * 2], mx[N * 2], f[N * 2];
pair<int, int> qu[N * 2];
inline int insert(int x, int las)
{
    if (ch[las][x] && len[las] + 1 == len[ch[las][x]])
        return ch[las][x];
    int cur = ++tot;
    int p = las;
    int clone = 0;
    int fl = 0;
    len[cur] = len[las] + 1;
    while (p && !ch[p][x])
    {
        ch[p][x] = cur;
        p = fa[p];
    }
    if (p == 0)
        fa[cur] = 1;
    else 
    {
        int q = ch[p][x];
        if (len[p] + 1 == len[q])
            fa[cur] = q;
        else 
        {
            if (p == las) fl = 1;
            clone = ++tot;
            len[clone] = len[p] + 1;
            fa[clone] = fa[q];
            memcpy(ch[clone], ch[q], sizeof(ch[q]));
            while (p && ch[p][x] == q)
            {
                ch[p][x] = clone;
                p = fa[p];
            }
            fa[q] = fa[cur] = clone;
        }
    }
    return fl ? clone : cur;
}

inline bool check(int k, int len)
{
    int res = 0;
    int l = 1, r = 0;
    qu[++r] = make_pair(0, 0);
    for (int i = 0; i < k; i++)
        f[i] = 0;
    for (int i = k; i <= len; i++)
    {
        while (l <= r && qu[l].second < i - mx[i])
            l++;
        if (l > r)
            f[i] = f[i - 1];
        else 
            f[i] = max(f[i - 1], qu[l].first + i);
        int t = i - k + 1;
        while (l <= r && qu[l].first <= f[t] - t)
            r--;
        qu[++r] = make_pair(f[t] - t, t);
        res = max(f[i], res);
    }
    double temp = res;
    double t = (double)(len) * 0.9;
    if (temp >= t)
        return true;
    else 
        return false;
}

signed main()
{
    int n, m;
    cin >> n >> m;
    tot = 1;
    for (int i = 1; i <= m; i++)
    {
        cin >> s + 1;
        int len = strlen(s + 1);
        int las = 1;
        for (int j = 1; j <= len; j++)
            las = insert(s[j] - '0', las);
    }
    for (int i = 1; i <= n; i++)
    {
        cin >> s + 1;
        int leng = strlen(s + 1);
        int cur = 1;
        int res = 0;
        for (int j = 1; j <= leng; j++)
        {
            int x = s[j] - '0';
            while (cur != 1 && !ch[cur][x])
            {
                cur = fa[cur];
                res = len[cur];
            }
            if (!ch[cur][x])
                mx[j] = 0;
            else 
            {
                res += 1;
                mx[j] = res;
                cur = ch[cur][x];
            }
        }
        int l = 0;
        int r = leng;
        while (l < r)
        {
            int mid = (l + r + 1) >> 1;
            if (check(mid, leng))
                l = mid;
            else 
                r = mid - 1;
        }
        cout << l << endl;
    }
    return 0;
}

SBzlt's Math note

2020-03-25 22:12:05 By DoorKickers

Binomial Coefficients

Let \begin{pmatrix}r \\ k \end{pmatrix} to denote that choose k elements from a set containing r terms without considering order.

With combinational interpretation. \begin{pmatrix} r \\ k \end{pmatrix} = \frac{n(n - 1) \dots (n - k + 1)}{k(k - 1)\dots 1} = \frac{n!}{(n - k)!\ k!}

There are some identities that hold when r is a positive real number and k is non-negative integer.

The first is symmetry identity {r \choose k} = {r \choose r - k} Next, we got *absorption identity* {r \choose k} = \frac{r}{k}{r - 1 \choose k - 1} If we multiply both sides by k, then we have k{r \choose k} = r{r - 1 \choose k - 1} There is a companion identity of the above one. \begin{align}(r - k){r \choose k} &= (r - k) {r \choose r - k} \\&=r{r - 1 \choose r - k - 1} \\&=r{r - 1 \choose k}\end{align} One of the other important formula is called *addition formula* {r \choose k} = {r - 1\choose k - 1} + {r - 1 \choose k} And it's easy to prove.

Beside these, we also have two identities of summation form. \sum_{k = -r}^{n}{{r + k \choose k}} = {r + n + 1 \choose n} and \sum_{k = 0}^{n}{k \choose m} = {n + 1 \choose m + 1} Notice that the summation 0 + 1 + 2 + \dots + n is a special case of former identity when m = 1 , and its result is equal to \frac{n(n+ 1)}{2} which is {n + 1 \choose 2} manipulated by that formula.

Here is one of the most important identity in binomial coefficients. (x + y)^r =\sum_{k}{{r \choose k}x^ky^{r - k}} And when x = y = 1 and r = n is nonnegative integer, we got {n \choose 0} + {n \choose 1} + {n \choose 2} + \dots + {n \choose n} = 2^n If we replace x by -1, then {n \choose 0} - {n \choose 1} + {n \choose 2} - {n \choose 3} + \dots + (-1)^n{n \choose n} = 0

Assumed that y = 1 and writing z instead of x to emphasize that an arbitrary complex number can be involved here. (1 + z)^r = \sum_{k}{{r \choose k}z^k}

二分图入门好题

2020-03-25 10:47:27 By DoorKickers

题目大意

给定一张有向无环图G, 求最长反链长大小

反链的定义为一个点集S, 其中任意点对(i,j) (i \neq j) 既不存在ij的路径, 也不存在ji的路径

最长反链是|S|最大的的反链

Solution

根据CGL第一零零八六定理

最长反链的大小 = 有向无环图最小路径可重复点覆盖


证明:

我们对G进行传递闭包, 得到的新图记为G'

问题转化为证明G'的最小路径点覆盖(路径不可重叠)为最长反链的大小

?藏身点集合为H, 路径集合为P, 显然, 每个路径中最多选取一个点, 所以|P| \geq |H|

接下来, 尝试构造一种情况使得|P| = |H|

?G_2'G' 的拆点二分图, 对于左侧非匹配点x_{out}, 说明这个点是路径的终点, 反复横跳后找到该路径的起点

?E to denote that the set of the end point. For each of them, we label all the point that current point can arrive, and denote as Nex.

If E \bigcap Nex = \emptyset , E is the answer we expected.

else we consider each point p in E \bigcap Nex, move along the path until current point p' \notin Nex, then remove p from E and add p' to E.

and repeat with new E until the condition is satisfied.

We can prove that we can find a legal p' in any time when the condition is not satisfied.

Assumed that we can not find that p', there must exsit a point q that can arrive whole path where includes p. And it also means that the P is not satisfied with the minimality which contradict with the premise.

BZOJ3771瞎写

2020-03-16 17:45:01 By DoorKickers

题目大意

我们讲一个悲伤的故事。 从前有一个贫穷的樵夫在河边砍柴。 这时候河里出现了一个水神,夺过了他的斧头,说: “这把斧头,是不是你的?” 樵夫一看:“是啊是啊!” 水神把斧头扔在一边,又拿起一个东西问: “这把斧头,是不是你的?” 樵夫看不清楚,但又怕真的是自己的斧头,只好又答:“是啊是啊!” 水神又把手上的东西扔在一边,拿起第三个东西问: “这把斧头,是不是你的?” 樵夫还是看不清楚,但是他觉得再这样下去他就没法砍柴了。 于是他又一次答:“是啊是啊!真的是!” 水神看着他,哈哈大笑道: “你看看你现在的样子,真是丑陋!” 之后就消失了。 樵夫觉得很坑爹,他今天不仅没有砍到柴,还丢了一把斧头给那个水神。 于是他准备回家换一把斧头。 回家之后他才发现真正坑爹的事情才刚开始。 水神拿着的的确是他的斧头。 但是不一定是他拿出去的那把,还有可能是水神不知道怎么偷偷从他家里拿走的。 换句话说,水神可能拿走了他的一把,两把或者三把斧头。 樵夫觉得今天真是倒霉透了,但不管怎么样日子还得过。 他想统计他的损失。 樵夫的每一把斧头都有一个价值,不同斧头的价值不同。总损失就是丢掉的斧头价值和。 他想对于每个可能的总损失,计算有几种可能的方案。 注意:如果水神拿走了两把斧头a和b,(a,b)和(b,a)视为一种方案。拿走三把斧头时,(a,b,c),(b,c,a),(c,a,b),(c,b,a),(b,a,c),(a,c,b)视为一种方案。

随笔瞎写

考虑一个生成函数 A = a_0 + a_1x + a_2x^2 + a_3x^3 + \cdots + a_nx^n 把系数看成方案数, 把指数看成价值

两个多项式的卷积为 C = A \times B = \sum_{i, j} a_i \cdot b_j 看一个简单的例子 A = 2x + x^2 \\B = 3x + 2x^2 \\C = 6x^2 + 4x^3 + 3x^3 + 2x^4 = 6x^2 + 7x^3+ 2x^4 显然, 卷积后的多项式可以代表从A, B 里各选出一个物品, 然后总价值达到某一个数的方案数

这个题中, 我们要从A里选出两个物品, 就可以自己对自己进行卷积 A = 2x + x^2 \\A \times A = 4x^2 + 4x^3 + x^4 但是存在一个物品被选中两次的情况

我们设B为每个物品取两次的情况 B = 2x^2 + x^4 显然, 最终答案为A \times A - B .

同理我们考虑从A中选3个物品, A \times A \times A = 发现可能会出现一个物品被选中多次的情况

枚举一下哪个物品被多选了, 设为x

发现可能会出现(x, x, x), (x, x, y), (y, x, x), (x, y, x) 这几种情况

我们设D = A \times B

发现可以得到(x, x, x), (y, x, x) 这两种情况

我们可以把结果乘三, 但发现(x, x, x) 减多了, 那就再加回来

C为选三个物品的情况, 再加回来两次就好了

最终答案是三种情况加起来, 又因为不同顺序算一种, 所以要除一下 A + \frac{A \times A - B}{2!} + \frac{A \times A \times A - A \times B \cdot 3 + 2 \cdot C }{3!}

答案就算出来了

用FFT或者NTT加速

复杂度是O(不会分析但好像能过把)

SB代码

/*
Date&Time: 2:37pm 3/16/2020]
*/

#include <bits/stdc++.h>

using namespace std;

#define int long long 

const int N = 3e6 + 10;
const int mod = 998244353;
const int g = 3;

int omega[N], omegaInv[N], rev[N], ta[N], tta[N], tab[N], a[N], b[N], c[N];
int n;

inline int fp(int a, int n, int mod) 
{
    int res = 1;
    while (n)
    {
        if (n & 1) res = (res * a) % mod;
        a = (a * a) % mod;
        n >>= 1;
    }
    return res % mod;
}

inline int inv(int x)
{
    return fp(x, mod - 2, mod);
}

inline void init(const int n)
{
    int k = 0;
    while ((1 << k) < n) k++;
    int x = fp(g, (mod - 1) / n, mod);
    omega[0] = omegaInv[0] = 1;
    for (int i = 1; i < n; i++)
    {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));
        omega[i] = (omega[i - 1] * x) % mod;
        omegaInv[i] = inv(omega[i]);
    }
}

inline void ntt(int *a, const int n, int *omega)
{
    for (int i = 0; i < n; i++)
    {
        if (i > rev[i]) continue;
        swap(a[i], a[rev[i]]);
    }
    for (int l = 2; l <= n; l *= 2)
    {
        int m = l / 2;
        for (int *p = a; p != a + n; p += l)
        {
            for (int i = 0; i < m; i++)
            {
                int t = (omega[n / l * i] * p[m + i]) % mod;
                p[m + i] = (p[i] - t + mod) % mod;
                p[i] = (p[i] + t) % mod;
            }
        }
    }
}

inline void dft(int *a, const int n)
{
    ntt(a, n, omega);
}

inline void idft(int *a, const int n)
{
    ntt(a, n, omegaInv);
    int x = inv(n);
    for (int i = 0; i < n; i++)
        a[i] = (a[i] * x) % mod;
}

inline void poly_multiply(int *a1, int n1, int *a2, int n2)
{
    int n = 1;
    while (n < n1 + n2 + 2) n *= 2;
    init(n);
    dft(a1, n), dft(a2, n);
    for (int i = 0; i < n; i++)
        a1[i] = (a1[i] * a2[i]) % mod;
    idft(a1, n);
    idft(a2, n);
}

inline void poly_squred(int *a, int n)
{
    int n1 = 1;
    while (n1 < n * 2 + 2) n1 *= 2;
    init(n1);
    dft(a, n1);
    for (int i = 0; i < n1; i++)
        a[i] = (a[i] * a[i]) % mod;
    idft(a, n1);
}

inline void poly_plus(int *a1, int n1 , int *a2, int n2, bool fl)
{
    int n = max(n1, n2);
    for (int i = 0; i <= n; i++)
        a1[i] = fl ? (a1[i] + a2[i]) % mod : (a1[i] - a2[i] + mod) % mod;
}

inline void poly_devide(int *a, int n, int v)
{
    int x = inv(v);
    for (int i = 0; i <= n; i++)
        a[i] = (a[i] * x) % mod;
}

inline void poly_times(int *a, int n, int v)
{
    for (int i = 0; i <= n; i++)
        a[i] = (a[i] * v) % mod;
}

int ans[N];

inline void print_poly(int* a, int n) 
{
    for (int i = 0; i <= n; i++)
    {
        cout << a[i] << ' ';
    }
    cout << endl;
}

signed main() 
{
    ios::sync_with_stdio(0);
    cin.tie(0);
    cin >> n;
    int siz = 0;
    for (int i = 1; i <= n; i++)
    {
        int x;
        cin >> x;
        siz = max(siz, x * 3);

        a[x]++, b[x * 2]++, c[x * 3]++;
        ta[x]++, tta[x]++, tab[x]++;

    }
    poly_squred(ta, siz);
    poly_multiply(tta, siz, ta, siz);
    poly_multiply(tab, siz, b, siz);
    poly_times(tab, siz, 3);
    poly_times(c, siz, 2);
    poly_plus(ans, siz, a, siz, 1);
    poly_plus(ta, siz, b, siz, 0);
    poly_devide(ta, siz, 2);
    poly_plus(ans, siz, ta, siz, 1);
    poly_plus(tta, siz, tab, siz, 0);
    poly_plus(tta, siz, c, siz, 1);
    poly_devide(tta, siz, 6);
    poly_plus(ans, siz, tta, siz, 1);
    for (int i = 0; i <= siz; i++)
    {
        if (ans[i])
            cout << i << ' ' << ans[i] << endl;
    }
    return 0;
}
/*
AC
*/

OrzCgl OrzSqy Orzxsj

后缀排序

2020-02-26 12:26:19 By DoorKickers

Pre-defination

S : a string named S that consists of several characters

|S| : the size of S

S[l:r] : a sub-string of S which consists of the characters S[l], ... , S[r]

Suf(i) : a suffix of S that means S[i:|S|]

Algorithm

Sqy Algorithm I is an algorithm that can sort the suffix of a given string S in O(nlogn) time, and the result of this algorithm will be stored in a array called Suffix Array.

Let me show the code to you first, and then I will explain some intricacy details

(Assumed that you have already got the Radix Sort which can sort an array in O(n) time.)

inline void SA_sort(char *s, int n) {
    for (int i = 1; i <= n; i++) ++c[x[i] = s[i]];
    for (int i = 2; i <= m; i++) c[i] += c[i - 1];
    for (int i = n; i >= 1; i--) sa[c[x[i]]--] = i;
    for (int k = 1; k <= n; k <<= 1) {
        int num = 0;
        for (int i = n - k + 1; i <= n; i++) y[++num] = i;
        for (int i = 1; i <= n; i++) {
            if (sa[i] > k)
                y[++num] = sa[i] - k;
        }
        for (int i = 1; i <= m; i++) c[i] = 0;
        for (int i = 1; i <= n; i++) ++c[x[i]];
        for (int i = 2; i <= m; i++) c[i] += c[i - 1];
        for (int i = n; i >= 1; i--) sa[c[x[y[i]]]--] = y[i], y[i] = 0;
        swap(x, y);
        x[sa[1]] = 1, num = 1;
        for (int i = 2; i <= n; i++) {
            x[sa[i]] = (y[sa[i]] == y[sa[i - 1]] && y[sa[i] + k] == y[sa[i - 1] + k]) ? num : ++num;
        }
        if (num == n) break;
        m = num;

    }
}

Ok, that's the code.

Above all, the most important part that can make you better understanding of that snippet is the meaning of each array.

To simplify, we can use the function concept instead of array concept.

A function f is a machine that can get input and return an output.

Good, looks like you understand it completely.

 

Next step : look through the arrays that mentioned in that snippet.

You will comprehend them via function concept (expect c[])

 

x[] : input a initial position i, and return the "value" of the first-key whose length is k and whose start position is i.

y[] : input a current partial rank k of second-keys, and return the initial position i of the first-key which corresponds to the No.k second-key.

sa[] : input a current global rank k of the combination of first-keys and it's corresponding second-key, and return the initial position i which represents for that rank k string.

Remember it!!!!!!!!!!!!!

Ok, the writter is fucked up.

So

To be continued.

最长公共上升子序列

2019-08-12 16:39:15 By DoorKickers

问题引入

给定两个序列A,B 求出其最长公共上升子序列

分析

本题是LIS,LCS问题的结合\require{enclose}\enclose{updiagonalstrike}{一看就是一道水题} 窝萌举个栗子

A_{1-4} = 2,2,1,3

B_{1-4} = 2,1,2,3

窝萌可以设A_0 = B_0 = -\infty

为什么呢(窝萌后面就知道了)

状态的定义:F_{ij} 表示 a以下标i结尾,b以下标j结尾,最长上升公共子序列的长度

如何得到状态转移方程呢

经过手算

F_{00} = 0 F_{01} = 0 F_{02} = 0 F_{03} = 0 F_{04} = 0

F_{10} = 0 \color{red}{F_{11} = 1} F_{12} = 1 F_{13} = 1 F_{14} = 1

F_{20} = 0 F_{21} = 1 F_{22} = 1 F_{23} = 1 F_{24} = 1

F_{30} = 0 F_{31} = 1 F_{32} = 1 F_{33} = 1 F_{34} = 1

F_{40} = 0 F_{41} = 1 F_{42} = 1 F_{43} = 1 \color{red}{F_{44} = 2}

观察红色部分

窝萌可以发现一个规律,F_{ij} 增加 的一个条件是 A_i = B_j , 否则就继承上一个状态

窝萌可以得到半个状态转移方程

if(A_i == B_j) F_{ij} = ......

else F_{ij} = F_{i-1j}

窝萌再去考虑另外一半方程

结合最长上升子序列的状态转移方程F_i=max\{F_{j} + 1\}(0 \leq j < i \mid A_j < A_i)

窝萌尝试推导一下最长公共上升子序列问题的状态转移方程

F_{ij}=max\{F_{i-1k} + 1\}(0 \leq k < j \mid B_k < B_i)

因为A_i 要和 B_j 配对 所以是F_{i-1k}

又因为窝萌把A_0 B_0 设为了-\infty

所以当A_i == B_j 时,能保证F_{ij} 至少为1 (想一想,为什么)

参考最长上升子序列,窝萌的目标是找到F_{ij}中最大的那个

最终算法复杂度为O(n^3)

本题还有一种优化方法,但对这个问题的完整讨论超过了本博客的范围,请读者朋友们自行上网上搜寻资料

到这里,本弱鸡第一篇博客就写完了

OrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOrzstOstOrzstOrzstOrzstOrzstOrzstO

共 9 篇博客