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