分类: OI笔记

  • CFGYM103415A CCPC2021广州A Math Ball

    题意: 有$n$种球,每种有无限个,同时第$i$种球有一个代价$c_i$,你要拿不超过$w$个球。如果最后第$i$种球你拿了$k_i$个,那么你会获得$\prod_{1\leq i\leq n}{k_i^{c_i}}$的权值,求所有合法方案的权值和。$n\leq 1e5,\sum{c_i}\leq 1e5,w\leq 10^{18}$

    $$  \text{考虑对于价值是}c_i\text{的球,构造生成函数}  \\  F_{c_i}\left( x \right) =\sum_{n\ge 0}{n^{c_i}x^n}  \\  \text{这样}\frac{\prod_i{F_{c_i}\left( x \right)}}{1-x}\text{的}w\text{次项即为答案}  \\  \text{设}F_k\left( x \right) =\sum_{n\ge 0}{n^kx^n},\text{显然可得}F_k\left( x \right) =x\frac{\mathrm{d}F_{k-1}\left( x \right)}{\mathrm{d}x}  \\  \text{进一步递推可得}, F_k\left( x \right) =\frac{\sum_{0\le i\le k}{T\left( k,i \right) x^i}}{\left( 1-x \right) ^{k+1}},\text{其中}T\left( k,i \right) \text{表示欧拉数}  \\  \text{考虑如何快速计算欧拉数}  \\  \text{首先由具体数学可得}  \\  \sum_{0\le i\le k}{T\left( k,i \right) \times \left( z+1 \right) ^i}=\sum_{0\le i\le k}{S_2\left( k,i \right) \times i!\times z^{k-i}}  \\  \text{进一步推导得欧拉数的通项公式}T\left( n,k \right) =\sum_{0\le i\le k}{\begin{array}{c}    \binom{n+1}{i}\left( -1 \right) ^i\left( k+1-i \right) ^n\\  \end{array}}  \\  \text{构造卷积形式}T\left( n,k \right) =\left( n+1 \right) !\sum_{0\le i\le k}{\frac{\left( -1 \right) ^i}{i!\left( n+1-i \right) !}}\times \left( k+1-i \right) ^n  \\  \text{卷积即可求出一行欧拉数。}  \\  \text{现在我们可以对每一个}c_i\text{计算出}F_{c_i}\left( x \right) \text{的分子,并且得到他们分子的乘积,设为}F\left( x \right)   \\  \text{设他们分母的乘积,再乘个}1-x\text{为}\left( 1-x \right) ^k  \\  \text{则答案即为}\frac{F\left( x \right)}{\left( 1-x \right) ^k}\text{的}w\text{次项系数}  \\  \text{即}\sum_{0\le i\le n}{\left[ x^i \right] F\left( x \right) \binom{w-i+k-1}{w-i}}\left( \text{广义二项式定理展开分母} \right)   \\  \text{令}a_i=\left[ x^i \right] F\left( x \right)   \\  \text{即}\sum_{0\le i\le n}{a_i\binom{w-i+k-1}{k-1}}  \\  \text{后者可以通过递推快速转移。}  \\    \\    \\  \frac{x^a}{\left( 1-x \right) ^b}=x^a\left( \sum_{n\ge 0}{\binom{-b}{n}\left( -1 \right) ^n}x^n \right)   \\  =x^a\left( \sum_{n\ge 0}{\frac{\left( -b \right) ^{\underline{n}}}{n!}\left( -1 \right) ^n}x^n \right)   \\  =x^a\left( \sum_{n\ge 0}{\frac{\left( -b-0 \right) \left( -b-1 \right) \left( -b-2 \right) ..\left( -b-\left( n-1 \right) \right)}{n!}\left( -1 \right) ^n}x^n \right)   \\  =x^a\left( \sum_{n\ge 0}{\frac{\left( b+0 \right) \left( b+1 \right) \left( b+2 \right) ..\left( b+\left( n-1 \right) \right)}{n!}}x^n \right)   \\  =x^a\left( \sum_{n\ge 0}{\frac{\left( b+n-1 \right) ^{\underline{n}}}{n!}}x^n \right)   \\  =\sum_{n\ge 0}{\frac{\left( b+n-1 \right) ^{\underline{n}}}{n!}x^{n+a}}  \\  =\sum_{n\ge 0}{\binom{b+n-1}{n}x^{n+a}}  \\  =\sum_{n\ge a}{\binom{b+n-a-1}{n-a}x^n}  \\  =\sum_{n\ge a}{\binom{b+n-a-1}{b-1}x^n}  \\  =\sum_{n\ge a}{\begin{array}{c}     \frac{\left( b-a+n-1 \right) ^{\underline{b-1}}}{\left( b-1 \right) !}x^n\\  \end{array}}  \\  \text{令}t=b-1  \\  \text{则}  \\  \sum_{n\ge a}{\begin{array}{c}      \frac{\left( t-a+n \right) ^{\underline{t}}}{t!}x^n\\  \end{array}}  \\  \text{从}a\text{推进到}a+1  \\  \frac{\left( t-a+n \right) ^{\underline{t}}}{\left( t-a+n-1 \right) ^{\underline{t}}}=\frac{t-a+n}{n-a}  \\    \\  10*9*8/\left( 9*8*7 \right) =10/7  \\    \\  \left( \frac{x^a}{\left( 1-x \right) ^b}\text{的}n\text{次项} \right)   \\    \\  \binom{n+b-1}{n}=\binom{n+b-1}{b-1}  \\  \binom{n+b-1}{n}\rightarrow \binom{n+b}{n+1}  \\    \\  \frac{\binom{n+b}{n+1}}{\binom{n+b-1}{n}}=\frac{n+b}{\left( n+1 \right)}  \\    \\  \frac{\left( n-2+b-i \right) !\left( b-1 \right) !\left( n-i \right) !}{\left( b-1 \right) !\left( n-i-1 \right) !\left( n-i+b-1 \right) !}  \\    \\  \frac{\left( n-i+b-2 \right) !\left( n-i \right)}{\left( n-i+b-1 \right) !}  \\  \frac{\left( n-i \right)}{\left( n-i+b-1 \right)}  \\    \\    $$
    #include <cinttypes>
    #include <cstring>
    #include <fstream>
    #include <iostream>
    #include <random>
    using int_t = int;
    using std::cin;
    using std::cout;
    using std::endl;
    
    const int_t LARGE = 2e6;
    
    using i64 = int64_t;
    #ifdef NTTCNT
    std::ofstream nttcnt("nttcnt.txt");
    #endif
    const int_t mod = 998244353;
    const int_t g = 3;
    
    int_t power(int_t b, int_t i) {
        int_t r = 1;
        if (i < 0)
            i = ((i64)i % (mod - 1) + mod - 1) % (mod - 1);
        while (i) {
            if (i & 1)
                r = (i64)r * b % mod;
            b = (i64)b * b % mod;
            i >>= 1;
        }
        return r;
    }
    
    void makeflip(int_t* arr, int_t size2) {
        int_t len = (1 << size2);
        arr[0] = 0;
        for (int_t i = 1; i < len; i++) {
            arr[i] = (arr[i >> 1] >> 1) | ((i & 1) << (size2 - 1));
        }
    }
    
    int_t upper2n(int_t x) {
        int_t r = 0;
        while ((1 << r) < x)
            r++;
        return r;
    }
    template <int_t arg = 1>
    void transform(int_t* A, int_t size2, int_t* flip) {
        // #define int_t i64
        int_t len = (1 << size2);
    #ifdef NTTCNT
        nttcnt << len << endl;
    #endif
        for (int_t i = 0; i < len; i++) {
            int_t r = flip[i];
            if (r > i)
                std::swap(A[i], A[r]);
        }
        for (int_t i = 2; i <= len; i *= 2) {
            int_t mr = power(g, (i64)arg * (mod - 1) / i);
            for (int_t j = 0; j < len; j += i) {
                int_t curr = 1;
                for (int_t k = 0; k < i / 2; k++) {
                    int_t u = A[j + k], v = (i64)curr * A[j + k + i / 2] % mod;
                    A[j + k] = ((i64)u + v) % mod;
                    A[j + k + i / 2] = ((i64)u - v + mod) % mod;
                    curr = (i64)curr * mr % mod;
                }
            }
        }
        // #undef int_t
    }
    void poly_mul(const int_t* A, int_t n, const int_t* B, int_t m, int_t* C) {
        /*
           计算n次多项式A与m次多项式B的乘积
        */
        int_t size2 = upper2n(n + m + 1);
        int_t len = 1 << size2;
        static int_t T1[LARGE], T2[LARGE];
        for (int_t i = 0; i < len; i++) {
            if (i <= n)
                T1[i] = A[i];
            else
                T1[i] = 0;
            if (i <= m)
                T2[i] = B[i];
            else
                T2[i] = 0;
        }
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        transform(T1, size2, fliparr);
        transform(T2, size2, fliparr);
        for (int_t i = 0; i < len; i++)
            T1[i] = (i64)T1[i] * T2[i] % mod;
        transform<-1>(T1, size2, fliparr);
        int_t inv = power(len, -1);
        for (int_t i = 0; i <= n + m; i++)
            C[i] = (i64)T1[i] * inv % mod;
    }
    
    const int_t FLARGE = 2e5 + 10;
    
    int_t fact[FLARGE + 1], inv[FLARGE + 1], factInv[FLARGE + 1];
    void euler_num(int_t n, int_t* out) {
        static int_t P1[LARGE];
        static int_t P2[LARGE];
        for (int_t i = 0; i <= n; i++) {
            P1[i] = (i64)((i % 2) ? (mod - 1) : 1) * factInv[i] % mod *
                    factInv[n + 1 - i] % mod;
            P2[i] = power(i, n);
        }
        poly_mul(P1, n, P2, n, out);
        for (int_t i = 0; i <= n; i++)
            out[i] = (i64)out[i] * fact[n + 1] % mod;
    }
    
    using poly_t = std::vector<int_t>;
    poly_t poly_dcmul(const poly_t* P, int_t left, int_t right) {
        if (left == right)
            return P[left];
        int_t mid = (left + right) / 2;
        auto lret = poly_dcmul(P, left, mid);
        auto rret = poly_dcmul(P, mid + 1, right);
        poly_t ret;
        ret.resize(lret.size() - 1 + rret.size() - 1 + 1);
        poly_mul(&lret[0], lret.size() - 1, &rret[0], rret.size() - 1, &ret[0]);
        while (!ret.empty() && ret.back() == 0)
            ret.pop_back();
        return ret;
    }
    int_t C(int_t n, int_t m) {
        return (i64)fact[n] * factInv[m] % mod * factInv[n - m] % mod;
    }
    i64 getC(i64 n, int_t m) {
        if (n < 0 || n < m)
            return 0;
        i64 prod = 1;
        for (int_t i = 0; i < m; i++)
            prod = prod * (n % mod - i + mod) % mod;
        return prod * factInv[m] % mod;
    }
    int main() {
        std::ios::sync_with_stdio(false);
        cin.tie(0), cout.tie(0);
        {
            fact[0] = fact[1] = inv[1] = factInv[0] = factInv[1] = 1;
            for (int_t i = 2; i <= FLARGE; i++) {
                fact[i] = (i64)fact[i - 1] * i % mod;
                inv[i] = (i64)(mod - mod / i) * inv[mod % i] % mod;
                factInv[i] = (i64)factInv[i - 1] * inv[i] % mod;
            }
        }
        int_t n;
        i64 w;
        cin >> n >> w;
    
        static poly_t up[int_t(1e5 + 10)];
        static poly_t cache[int_t(1e5 + 10)];
        int_t sum = 0;
        for (int_t i = 1; i <= n; i++) {
            int_t c;
            cin >> c;
            sum += c + 1;
            poly_t& curr = up[i];
            if (!cache[c].empty()) {
                curr = cache[c];
                continue;
            }
            curr.resize(2 * c + 3);
            euler_num(c, &curr[0]);
            curr.resize(c + 1);
            cache[c] = curr;
    #ifdef DEBUG
            cout << "up " << i << " = ";
            for (auto t : curr)
                cout << t << " ";
            cout << endl;
    #endif
        }
        sum++;
        // for (int_t i = 0; i <= sum; i++) {
        //     down[i] = (i64)((i % 2) ? (mod - 1) : 1) * C(sum, i) % mod;
        // }
        auto upprod = poly_dcmul(up, 1, n);
        // cout << upprod.size() << endl;
        upprod.resize(sum + 1);
    #ifdef DEBUG
        cout << "up prod = ";
        for (auto t : upprod)
            cout << t << " ";
        cout << endl;
        cout << "sum = " << sum << endl;
    #endif
        // cout << "down size = " << sum << endl;
        // w -= cutcount;
        // if (w < 0) {
        //     cout << 0 << endl;
        //     return 0;
        // }
        // int_t ret = poly_divat(&upprod[0], down, upprod.size() - 1, sum, w);
        i64 ret = 0;
        {
            i64 b = sum;
            i64 n = w;
            i64 curr = getC(n - 0 + b - 1, b - 1);
    // C(n-i+b-1,b-1)到C(n-(i+1)+b-1,b-1)
    // getC(n - i + b - 1, b - 1)
    #ifdef DEBUG
            cout << "init curr = " << curr << endl;
    #endif
            for (i64 i = 0; i < upprod.size(); i++) {
                ret = (ret + upprod[i] * curr % mod) % mod;
                curr = curr * (n % mod - i + mod) % mod *
                       power((n % mod - i + b - 1 + (i64)3 * mod) % mod, -1) % mod;
                if (curr == 0 || n <= i - n + 1)
                    break;
    #ifdef DEBUG
                cout << "curr = " << curr << " at i = " << i << " with n = " << n
                     << endl;
    #endif
            }
            // int_t curr = getC(w - b + b - 1, b - 1);
            // for (i64 i = std::max<i64>(0, w - sum); i < upprod.size(); i++) {
            //     ret = (ret + curr * upprod[i] % mod) % mod;
            //     curr = curr * (n % mod + n % mod - sum % mod + mod + i) % mod *
            //            power(n + 1, -1) % mod;
            // }
            // i64 n = w;
            // i64 t = sum - 1;
            // i64 prod = 1;
            // // i64 curr = (t + n) % mod;
            // for (int_t i = 0; i < t; i++)
            //     prod = (prod * (t + n % mod - i) % mod);
            // for (int_t a = 0; a < upprod.size(); a++) {
            //     if (n >= a) {
            //         ret = (ret + prod * factInv[t] % mod * upprod[a] % mod) %
            //         mod;
            //     }
            //     prod = prod * (t - a + n % mod) % mod *
            //            power((n % mod - a + mod) % mod, -1) % mod;
            // }
        }
        cout << ret << endl;
        return 0;
    }

     

  • 跑的比较快的多项式板子

    #include <cstring>
    #include <iostream>
    #include <random>
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    
    const int_t LARGE = 2e6;
    
    const int_t mod = 998244353;
    const int_t g = 3;
    int_t power(int_t b, int_t i) {
        int_t r = 1;
        if (i < 0)
            i = (i % (mod - 1) + mod - 1) % (mod - 1);
        while (i) {
            if (i & 1)
                r = r * b % mod;
            b = b * b % mod;
            i >>= 1;
        }
        return r;
    }
    
    void makeflip(int_t* arr, int_t size2) {
        int_t len = (1 << size2);
        arr[0] = 0;
        for (int_t i = 1; i < len; i++) {
            arr[i] = (arr[i >> 1] >> 1) | ((i & 1) << (size2 - 1));
        }
    }
    
    int_t upper2n(int_t x) {
        int_t r = 0;
        while ((1 << r) < x)
            r++;
        return r;
    }
    template <int_t arg = 1>
    void transform(int_t* A, int_t size2, int_t* flip) {
        int_t len = (1 << size2);
        for (int_t i = 0; i < len; i++) {
            int_t r = flip[i];
            if (r > i)
                std::swap(A[i], A[r]);
        }
        for (int_t i = 2; i <= len; i *= 2) {
            int_t mr = power(g, arg * (mod - 1) / i);
            for (int_t j = 0; j < len; j += i) {
                int_t curr = 1;
                for (int_t k = 0; k < i / 2; k++) {
                    int_t u = A[j + k], v = curr * A[j + k + i / 2] % mod;
                    A[j + k] = (u + v) % mod;
                    A[j + k + i / 2] = (u - v + mod) % mod;
                    curr = curr * mr % mod;
                }
            }
        }
    }
    void poly_inv(const int_t* A, int_t n, int_t* result) {
        /*
        这里的n是模x^n
        计算B(x)*A(x) = 1 mod x^n, 其中A(x)已知
        假设已知A(x)*B(x) = 1 mod x^{ceil(n/2)}
        假设C(x)*A(x) = 1 mod x^n
        (A(x)B(x)-1)^2 = A^2(x)B^2(x)-2A(x)B(x)+1= 0
        A(x)B^2(x)-2B(x)+C(x) = 0
        C(x) = 2B(x)-A(x)B^2(x)
        */
        // #ifdef DEBUG
        //     cout << "at " << n << endl;
        // #endif
        if (n == 1) {
            result[0] = power(A[0], -1);
            return;
        }
        int_t next = n / 2 + n % 2;
        poly_inv(A, next, result);
        //次数不要选错了,应该用n次的A和B去卷
        int_t size2 = upper2n(n + 2 * next + 1);
        static int_t X[LARGE];
        static int_t Y[LARGE];
        int_t len = (1 << size2);
        // 写错设置范围了!
        memset(X + n, 0, sizeof(X[0]) * (len - n));
        memset(Y + next, 0, sizeof(Y[0]) * (len - next));
        memcpy(X, A, sizeof(A[0]) * n);
        memcpy(Y, result, sizeof(result[0]) * next);
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        transform<1>(X, size2, fliparr);
        transform<1>(Y, size2, fliparr);
    
        for (int_t i = 0; i < len; i++) {
            X[i] = (2 * Y[i] - X[i] * Y[i] % mod * Y[i] % mod + mod) % mod;
        }
        transform<-1>(X, size2, fliparr);
        const int_t inv = power(len, -1);
        for (int_t i = 0; i < n; i++)
            result[i] = X[i] * inv % mod;
    #ifdef DEBUG
        cout << "poly inv at " << n << endl;
        cout << "result = ";
        for (int_t i = 0; i < next; i++)
            cout << result[i] << " ";
        cout << endl;
    
    #endif
    }
    int_t poly_divat(const int_t* F, const int_t* G, int_t n, int_t k) {
        /*
            n次多项式F和G
            计算F(x)/G(x)的k次项前系数
            考虑F(x)*G(-x)/G(x)*G(-x),分母只有偶数次项,写为C(x^2);分子写成xA(x^2)+B(x^2),如果k是奇数,那么递归(A,C,n,k/2),如果k是偶数,那么递归(B,C,n,k/2)
            到k<=n时直接计算
        */
        int_t size2 = upper2n(2 * n + 1);
        int_t len = 1 << size2;
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        static int_t T1[LARGE], T2[LARGE], T3[LARGE];
        for (int_t i = 0; i < len; i++) {
            if (i <= n)
                T1[i] = F[i], T2[i] = G[i];
            else
                T1[i] = T2[i] = T3[i] = 0;
        }
        const int_t inv = power(len, -1);
        while (k >= n) {
    #ifdef DEBUG
            cout << "curr at k = " << k << endl;
    #endif
            for (int_t i = 0; i < len; i++) {
                if (i <= n) {
                    T3[i] = T2[i] * (i % 2 ? (mod - 1) : 1);
                } else
                    T3[i] = 0;
            }
            transform(T1, size2, fliparr);
            transform(T2, size2, fliparr);
            transform(T3, size2, fliparr);
            for (int_t i = 0; i < len; i++) {
                T1[i] = T1[i] * T3[i] % mod;
                T2[i] = T2[i] * T3[i] % mod;
            }
            transform<-1>(T1, size2, fliparr);
            transform<-1>(T2, size2, fliparr);
    #ifdef DEBUG
            cout << "prod T1 = ";
            for (int_t i = 0; i < len; i++)
                cout << T1[i] * inv % mod << " ";
            cout << endl;
            cout << "prod T2 = ";
            for (int_t i = 0; i < len; i++)
                cout << T2[i] * inv % mod << " ";
            cout << endl;
    
    #endif
            for (int_t i = 0; i < len; i++) {
                if (i * 2 < len) {
                    T2[i] = T2[i * 2] * inv % mod;
                } else
                    T2[i] = 0;
            }
            int_t b = k % 2;
            for (int_t i = 0; i < len; i++) {
                if (i % 2 == b) {
                    T1[i / 2] = T1[i] * inv % mod;
    #ifdef DEBUG
                    cout << "at " << i << " assgin " << T1[i] * inv << " to "
                         << i / 2 << endl;
    #endif
                }
    #ifdef DEBUG
                cout << "assign 0 to " << i << endl;
    #endif
                if (i > 0)
                    T1[i] = 0;
            }
    #ifdef DEBUG
            cout << "finished T1 = ";
            for (int_t i = 0; i < len; i++)
                cout << T1[i] % mod << " ";
            cout << endl;
            cout << "finished T2 = ";
            for (int_t i = 0; i < len; i++)
                cout << T2[i] % mod << " ";
            cout << endl;
    
    #endif
            k >>= 1;
        }
    
        poly_inv(T2, k + 1, T3);
        // #ifdef DEBUG
    
        //     cout << "finished k = " << k << endl;
        //     cout << "T1 = ";
        //     for (int_t i = 0; i < len; i++)
        //         cout << T1[i] % mod << " ";
        //     cout << endl;
        //     cout << "T2 = ";
        //     for (int_t i = 0; i < len; i++)
        //         cout << T2[i] % mod << " ";
        //     cout << endl;
        //     cout << "T2 inv = ";
        //     for (int_t i = 0; i < len; i++)
        //         cout << T3[i] % mod << " ";
        //     cout << endl;
    
        // #endif
        int_t result = 0;
        //计算结果的k次项
        for (int_t i = 0; i <= k; i++) {
            result = (result + T1[i] * T3[k - i] % mod) % mod;
        }
        return result;
    }
    void poly_mul(const int_t* A, int_t n, const int_t* B, int_t m, int_t* C) {
        /*
           计算n次多项式A与m次多项式B的乘积
        */
        int_t size2 = upper2n(n + m + 1);
        int_t len = 1 << size2;
        static int_t T1[LARGE], T2[LARGE];
        for (int_t i = 0; i < len; i++) {
            if (i <= n)
                T1[i] = A[i];
            else
                T1[i] = 0;
            if (i <= m)
                T2[i] = B[i];
            else
                T2[i] = 0;
        }
        static int_t fliparr[LARGE];
        makeflip(fliparr, size2);
        transform(T1, size2, fliparr);
        transform(T2, size2, fliparr);
        for (int_t i = 0; i < len; i++)
            T1[i] = T1[i] * T2[i] % mod;
        transform<-1>(T1, size2, fliparr);
        int_t inv = power(len, -1);
        for (int_t i = 0; i <= n + m; i++)
            C[i] = T1[i] * inv % mod;
    }
    int_t poly_linear_rec(const int_t* A0, const int_t* F0, int_t n, int_t k) {
        /*
            计算线性递推
            F[1],F[2]...F[k] 递推系数
            A[0],A[1]...A[k-1] 首项
    
            A[m]=A[m-1]*F[1]+A[m-2]*F[2]+...+A[m-k]*F[k]
        */
    
        static int_t T1[LARGE], T2[LARGE];
        static int_t Ax[LARGE], Fx[LARGE];
        Fx[0] = 0;
        Ax[k] = 0;
        for (int i = 1; i <= k; i++) {
            Fx[i] = F0[i];
            Ax[i - 1] = A0[i - 1];
        }
        poly_mul(Ax, k, Fx, k, T1);
        T1[0] = Ax[0];
        for (int i = 1; i <= k - 1; i++) {
            T1[i] = (Ax[i] - T1[i] + mod) % mod;
        }
        for (int i = k; i <= 2 * k; i++)
            T1[i] = 0;
        T2[0] = 1;
        for (int i = 1; i <= k; i++)
            T2[i] = (mod - Fx[i]) % mod;
    
        // #ifdef DEBUG
        //     cout << "T1 = ";
        //     for (int_t i = 0; i <= k; i++) {
        //         cout << T1[i] << " ";
        //     }
        //     cout << endl;
        //     cout << "T2 = ";
        //     for (int_t i = 0; i <= k; i++) {
        //         cout << T2[i] << " ";
        //     }
        //     cout << endl;
    
        // #endif
        return poly_divat(T1, T2, k, n);
    }
    
    void poly_log(const int_t* A, int_t n, int_t* out) {
        /*
            计算log(A(x)), A(x)为n次多项式
            DlogF(x) =DF(x)/F(x)
        */
        static int_t Ad[LARGE];
        static int_t Finv[LARGE];
        static int_t R[LARGE];
        const int_t m = n - 1;
        for (int_t i = 0; i <= m; i++) {
            Ad[i] = A[i + 1] * (i + 1) % mod;
        }
        Ad[n] = 0;
        poly_inv(A, n + 1, Finv);
        poly_mul(Ad, m, Finv, n, R);
        for (int_t i = 1; i <= n; i++) {
            out[i] = R[i - 1] * power(i, -1) % mod;
        }
    }
    void poly_exp(const int_t* A, int_t n, int_t* out) {
        /*
            计算exp(A(x)), A(x)为n次多项式
            H(x)=G(x)(1-logG(x)+A(x))
            H(x)为当次递推项,G(x)为上一次递推项
        */
        if (n == 1) {
            out[0] = 1;
            return;
        }
        int_t r = n / 2 + n % 2;
        poly_exp(A, r, out);
        static int_t G[LARGE];
        static int_t G2[LARGE];
        static int_t logG[LARGE];
        static int_t R[LARGE];
        for (int_t i = 0; i < r; i++)
            G[i] = out[i];
        for (int_t i = r; i < n; i++)
            G[i] = 0;
        poly_log(G, n - 1, logG);
        // for (int_t i = r; i < n; i++)
        //     logG[i] = 0;
        for (int_t i = 0; i < n; i++) {
            G2[i] = (mod - logG[i] + A[i]) % mod;
        }
        G2[0] = (G2[0] + 1) % mod;
        poly_mul(G, n - 1, G2, n - 1, R);
        for (int_t i = 0; i < n; i++)
            out[i] = R[i];
    #ifdef DEBUG
        cout << "at " << n << endl;
        cout << "A = ";
        for (int_t i = 0; i < n; i++)
            cout << A[i] << " ";
        cout << endl;
        cout << "G = ";
        for (int_t i = 0; i < n; i++)
            cout << G[i] << " ";
        cout << endl;
        cout << "logG = ";
        for (int_t i = 0; i < n; i++)
            cout << logG[i] << " ";
        cout << endl;
        cout << "G2 = ";
        for (int_t i = 0; i < n; i++)
            cout << G2[i] << " ";
        cout << endl;
        cout << "R = ";
        for (int_t i = 0; i < n; i++)
            cout << R[i] << " ";
        cout << endl;
        cout << endl;
    #endif
    }
    int main() {
        std::ios::sync_with_stdio(false);
        cin.tie(0), cout.tie(0);
    
        int_t n;
        cin >> n;
        static int_t A[LARGE], B[LARGE];
        static int_t C[LARGE];
        for (int_t i = 0; i < n; i++)
            cin >> A[i];
        poly_exp(A, n, B);
        for (int_t i = 0; i < n; i++)
            cout << B[i] << " ";
    
        return 0;
    }

     

  • CFGYM103104I-WHUPC I Sequence

    题意:有个$n\times n,n\leq 5000$的矩阵,挖掉$m,m\leq 10^6$个格子,问不包括这些格子的子矩阵个数。

    做法很容易考虑:枚举一行,这一行从左往右扫,依次计算每个格子为右下角的子矩阵个数。

    扫的时候维护一个单调栈,里面存(元素,以这个元素为高度的格子最长向右延伸了多少格)

    为什么要存第二个项呢?我们维护的单调栈里面,每个东西实际上表示的是一个矩形,含义是:在我们处理完当前这个列后,这些矩形里的元素都可以作为合法的子矩阵左上顶点。

    扫的时候维护一个东西:单调栈里所有矩形的大小和$sum$,即合法的左上顶点个数。

    每次加入一个元素$x$,分以下三种情况讨论:

    • $x$大于单调栈最后的元素:直接加入。长度记为1,sum加上$x$(显然多了这么多合法的格子)
    • $x$等于单调栈最后的元素:单调栈最后的元素的长度加1,sum加上$x$
    • $x$小于单调栈最后的元素:弹出单调栈最后的元素,并把其长度加到倒数第二个元素上,然后重复这三个check

    最后我们每插入一个元素后,所维护的sum就是以这个点为右下角的答案。

    另外有一个坑:输入1625 0,答案输出正确,但是从1626 0开始,答案越来越小,看起来像是溢出了,但是开-fsanitize=undefined没报任何问题,结果最后查出来了,算一行的答案时使用的是std::accumlate,初始值传的是int类型,而后累加变量自动推导为int,导致这个函数里面溢出了,然后由于ubsan并不会对包含的其他文件的代码做检查,于是挂掉了。

     

    #include <assert.h>
    #include <algorithm>
    #include <iostream>
    #include <numeric>
    #include <utility>
    #include <vector>
    using int_t = long long;
    using std::cin;
    using std::cout;
    using std::endl;
    /**
     * 1624 0 -> R
     * 1625 0 -> R
     * 1626 0 -> E 正确1749670208001 本程序 1736785306113
     * 1627 0 -> E
     */
    bool block[5001][5001];
    int_t upmost[5001][5001];
    int_t n, m;
    int main() {
        std::ios::sync_with_stdio(false);
        cin >> n >> m;
        for (int_t i = 1; i <= m; i++) {
            int_t r, c;
            cin >> r >> c;
            block[r][c] = true;
        }
        for (int_t i = 1; i <= n; i++) {  //列
            int_t lastpos = 0;
            for (int_t j = 1; j <= n; j++) {
    #ifdef DEBUG
                // cout << "block " << j << " " << i << " = " << block[j][i] <<
                // endl;
    #endif
                if (block[j][i])
                    lastpos = j;
                upmost[j][i] = lastpos;
    #ifdef DEBUG
                cout << "upmost " << j << " " << i << " = " << upmost[j][i] << endl;
    #endif
            }
        }
        int_t result = 0;
        //枚举底边行号
        for (int_t i = 1; i <= n; i++) {
            static int_t answer[5001];  // j->以(i,j)为右下角的矩形个数
            //以当前元素为右下角的矩形个数
            int_t sum = 0;
            // first为元素,second为出现次数
            std::vector<std::pair<int_t, int_t>> stack;
            stack.emplace_back(0, 0);
            for (int_t j = 1; j <= n; j++) {
                int_t x = i - upmost[i][j];  //向上延伸的空格子数(包括i,j)
                int_t count = 1;
                while (x < stack.back().first) {
                    count += stack.back().second;
                    sum -= stack.back().first * stack.back().second;
                    stack.pop_back();
                }
                assert(sum >= 0);
                if (x == stack.back().first) {
                    stack.back().second += count;
                } else {
                    stack.emplace_back(x, count);
                }
                sum += x * count;
                answer[j] = sum;
                // cout << "answer col " << j << " = " << answer[j] << endl;
    #ifdef DEBUG
                cout << "pushed col " << j << " val " << x << " sum = " << sum
                     << " answer = " << answer[j] << endl;
    #endif
            }
            int_t currrow = std::accumulate(answer + 1, answer + 1 + n, (int_t)0);
            result += currrow;
            // cout << "answer row " << i << " = " << currrow << " " << currrow / i
            //      << " " << currrow % i << endl;
    #ifdef DEBUG
            cout << "answer row " << i << " = " << currrow << endl;
    #endif
        }
        cout << result << endl;
        return 0;
    }

     

  • CF1521C Nastia and a Hidden Permutation

    真是个沙比提,比赛的时候也想出来了,但是因为太迷糊了没写完

    首先我们先考虑如何确定下最大值所在的位置。

    执行询问$max(min(n-1,p_i),min(n,p_{i+1}))$,这个询问实际上等价于$max(min(n-1,p_i),p_{i+1})$。

    • 如果结果是n,那么说明$p_{i+1}=n$(很显然,max里第一项不会超过n-1)。
    • 如果结果比n-1要小,那么说明$p_i$和$p_{i+1}$都不是n。
    • 如果结果是n-1,那么就要讨论下了。这个时候n可能出现在这两个数里,也可能没出现。

    我们再构造另一个询问:$min(max(n-1,p_i),max(n,p_{i+1}))$,显然这个询问等价于$max(n-1,p_i)$,所以我们可以利用这个询问来检查$p_i$的值是不是n。

    在找到最大值的位置,设它为$maxpos$,那么我们对于每一个$i\neq maxpos$,构造询问$min(max(1,p_i),max(2,p_{maxpos}))$,这个询问等价于$p_i$,然后就可以求出来每个位置的值了。

    #include <algorithm>
    #include <iostream>
    using std::cin;
    using std::cout;
    using std::endl;
    using int_t = long long;
    int_t result[int_t(1e5)];
    int main() {
        std::ios::sync_with_stdio(false);
        int_t T;
        cin >> T;
        while (T--) {
            int_t n;
            cin >> n;
            int_t maxpos = -1;
            const auto check = [&](int_t i, int_t j) {
                cout << "? 1 " << i << " " << j << " " << n - 1 << endl;
                int_t x;
                cin >> x;
                if (x == n) {
                    maxpos = j;
                    return;
                }
                if (x < n - 1)
                    return;
                // x==n-1
                for (int_t _ = 1; _ <= 2; _++) {
                    cout << "? 2 " << i << " " << j << " "<< n - 1 << endl;
                    cin >> x;
                    if (x == n) {
                        maxpos = i;
                    }
                    std::swap(i, j);
                }
            };
            for (int_t i = 1; i <= n && i + 1 <= n && maxpos == -1; i += 2) {
                int_t j = i + 1;
    
                check(i, j);
            }
            if (maxpos == -1) {
                maxpos = n;
            }
            result[maxpos] = n;
            for (int_t i = 1; i <= n; i++) {
                if (i == maxpos)
                    continue;
                cout << "? 2 " << i << " " << maxpos << " 1" << endl;
                ;
                int_t x;
                cin >> x;
                result[i] = x;
            }
            cout << "! ";
            for (int_t i = 1; i <= n; i++) {
                cout << result[i] << " ";
            }
            cout << endl;
        }
        return 0;
    }
    

     

  • CFGYM102394A CCPC2019哈尔滨A Artful Paintings

    一眼看过去是差分约束板子题,实际上也就是差分约束板子题。

    首先这个题可以对答案进行二分。如果铺x个可行,那么x+1个显然可行,所以我们先二分确定一个x。

    然后我们的题目转变成了,能不能在恰好铺x个的情况下构造一组合法方案,根据最短路的松弛条件可知,不等式$x_a+v\leq x_b$等价于从a到b,边权为$-v$的边(最短路存在时,令$x_a$表示源点到a的距离,$x_b$同理,那么$x_a$和$x_b$在作为未知变量的时候一定是满足上述不等式的。)

    所以我们可以得到,需要建立以下六种不等式:

    1. $x_{n-1}+0\leq x_n$(前缀和不减)
    2. $x_n+(-1)\leq x_{n-1}$前缀和要么加1要么不变
    3. $x_{l-1}+x\leq x_r$ 区间染色数下限
    4. $x_r +(x-M) \leq  x_{l-1}$ 区间外染色数下限
    5. $x_0+M\leq x_n$ 一共至少染M个
    6. $x_n+(-M)\leq x_0$ 一共至多染M个(与5一起保证恰好染M个)

    然后据此建边,其中$M$是我们二分的铺的个数。建出来的图跑一个从n到0的最短路,如果最短路存在(即没有负环),那说明M可行,答案可以考虑进一步缩小。如果存在负环,则说明M不可行。

    但是这里的SPFA要注意优化,大概需要一个SLF(Small Label First),即如果待入队节点的已知最短距离小于队列头部节点,则把他扔到队列头部。
    然后别用long long,跑的有点慢。
    #include <algorithm>
    #include <deque>
    #include <iostream>
    #include <vector>
    using std::cin;
    using std::cout;
    using std::endl;
    
    using int_t = int;
    
    /*
    1. x_{n-1}+0<=x_n(前缀和不减)
    2. x_n+(-1)<=x_{n-1}前缀和要么加1要么不变
    3. x_{l-1}+x<=x_r 区间染色数下限
    4. x_r +(x-M) <= x_{l-1} 区间外染色数下限
    5. x_0+M<=x_n 一共至少染M个
    6. x_n+(-M)<=x_0 一共至多染M个(与5一起保证恰好染M个)
    
    */
    
    const int_t LARGE = 3010;
    const int_t INF = 0x3fffffff;
    int_t dis[LARGE + 1];
    int_t inqN[LARGE];
    bool inqueue[LARGE + 1];
    struct Edge {
        int_t to, weight;
        Edge(int_t v, int_t w) : to(v), weight(w) {}
    };
    struct Limitation {
        int_t left, right, val;
    } cond1[LARGE + 1], cond2[LARGE + 1];
    
    int_t n, m1, m2;
    //不等式A+X<=B变为边(A,B),长度-X
    std::vector<Edge> graph[LARGE + 1];
    
    void clear() {
        std::fill(dis, dis + 1 + n, INF);
        std::fill(inqueue, inqueue + 1 + n, false);
        std::fill(inqN, inqN + 1 + n, 0);
        for (int_t i = 0; i <= n; i++)
            graph[i].clear();
    }
    
    void pushEdge(int_t u, int_t v, int_t w) {
    #ifdef DEBUG
        cout << "edge " << u << " to " << v << " val " << w << endl;
    #endif
        graph[u].push_back(Edge(v, w));
    }
    //添加不等式x_a+v<=x_b
    void pushLEQ(int_t a, int_t b, int_t v) {
    #ifdef DEBUG
        cout << "cond x_" << a << " + " << v << " <= x_" << b << endl;
    #endif
        pushEdge(b, a, -v);
    }
    
    bool SPFA(int_t src, int_t target) {
        std::deque<int_t> queue;
        dis[src] = 0;
        inqueue[src] = true;
        queue.push_back(src);
        while (!queue.empty()) {
            int_t front = queue.front();
            queue.pop_front();
            inqueue[front] = false;
            inqN[front]++;
            if (inqN[front] > n)
                return false;
            for (const auto& edge : graph[front]) {
                if (dis[edge.to] > dis[front] + edge.weight) {
                    dis[edge.to] = dis[front] + edge.weight;
    #ifdef DEBUG
                    cout << "dis " << edge.to << " to " << dis[edge.to] << endl;
    #endif
                    if (!inqueue[edge.to]) {
                        if (!queue.empty() && dis[edge.to] < dis[queue.front()])
                            queue.push_front(edge.to);
                        else
                            queue.push_back(edge.to);
                        inqueue[edge.to] = true;
                    }
                }
            }
        }
        return true;
    }
    
    bool check(int_t x) {
    #ifdef DEBUG
        cout << "checking " << x << endl;
    #endif
        clear();
        for (int_t i = 1; i <= n; i++) {
            pushLEQ(i - 1, i, 0);
        }
        for (int_t i = 1; i <= n; i++) {
            pushLEQ(i, i - 1, -1);
        }
        for (int_t i = 1; i <= m1; i++) {
            const auto& ref = cond1[i];
            pushLEQ(ref.left - 1, ref.right, ref.val);
        }
        for (int_t i = 1; i <= m2; i++) {
            const auto& ref = cond2[i];
            pushLEQ(ref.right, ref.left - 1, ref.val - x);
        }
        pushLEQ(0, n, x);
        pushLEQ(n, 0, -x);
        auto ret = SPFA(n, 0);
    #ifdef DEBUG
        cout << "check " << x << " = " << ret << endl;
    #endif
        return ret;
    }
    
    int main() {
        std::ios::sync_with_stdio(false);
        int_t T;
        cin >> T;
        while (T--) {
            cin >> n >> m1 >> m2;
            for (int_t i = 1; i <= m1; i++) {
                auto& ref = cond1[i];
                cin >> ref.left >> ref.right >> ref.val;
            }
            for (int_t i = 1; i <= m2; i++) {
                auto& ref = cond2[i];
                cin >> ref.left >> ref.right >> ref.val;
            }
            int_t left = 0, right = n;
            int_t result = INF;
            while (left + 1 < right) {
                int_t mid = (left + right) / 2;
                if (check(mid)) {
                    result = std::min(result, mid);
                    right = mid - 1;
                } else
                    left = mid + 1;
            }
            for (int_t i = left; i <= right; i++)
                if (check(i))
                    result = std::min(result, i);
            cout << result << endl;
        }
    
        return 0;
    }

     

  • CFGYM101982F (2018 ICPC Pacific Northwest Regional Contest – Rectangles)

    题意:给定n(1e5)个矩形,边与坐标轴平行,四个点的坐标都是整数(1e9级别),问被矩形覆盖了奇数次的面积和。

    首先把所有的矩形拆成上边和下边两条边,这样我们就得到了若干条与x轴平行的线段。然后我们把这些线段按照y坐标从小到大排个序。

    现在我们维护一棵线段树,下标表示的是(离散化后的)整个横轴的覆盖情况。以及为了方便起见,整个线段树上所表示的区间都是左闭右开的(一条线段所能表示的实际长度是右端点减左端点)。

    现在我们按照y坐标从小到大枚举每一条横线,对于一条横线,我们把它在线段树上所覆盖的的区域异或上1,同时面积取反(即用线段长度减掉区域内被标记过的线段长度和,此操作即为异或),由于我们存储的都是左闭右开区间,所以直接用右端点对应的数减掉左端点对应的数即可。

    最后每插入一条线段,我们就查一下现在线段树上被覆盖的区间长度,然后乘上下一条线段的高度减掉当前线段的高度(此方法是可以正常处理有相同高度的线段的,因为他们高度相同,所以算出来都是0)。

    #include <algorithm>
    #include <iostream>
    #include <vector>
    
    using int_t = long long int;
    
    using std::cin;
    using std::cout;
    using std::endl;
    
    struct Segment {
        int_t left, right, height;
    };
    struct Node {
        Node *left = nullptr, *right = nullptr;
        int_t begin, end;
        bool flag = false;
        int_t sum = 0;
        const std::vector<int_t>& vals;
        Node(int_t begin, int_t end, const std::vector<int_t>& vals)
            : begin(begin), end(end), vals(vals) {}
        void swap() {
            flag ^= 1;
            sum = vals[end] - vals[begin] - sum;
        }
        void pushdown() {
            if (flag) {
                left->swap();
                right->swap();
                flag = 0;
            }
        }
        void maintain() { sum = left->sum + right->sum; }
        void insert(int_t begin, int_t end) {  //插入一条线段
    #ifdef DEBUG
            cout << "insert " << begin << "," << end << " to " << this->begin << ","
                 << this->end << endl;
    #endif
            if (this->begin >= begin && this->end <= end) {
                this->swap();
                return;
            }
            if (this->begin >= end || this->end <= begin) {
                return;
            }
            int_t mid0 = (this->begin + this->end) / 2;
            pushdown();
            left->insert(begin, end);
            right->insert(begin, end);
            maintain();
        }
    };
    Node* build(int_t begin, int_t end, const std::vector<int_t>& vals) {
        Node* node = new Node(begin, end, vals);
        if (begin + 1 != end) {
            int_t mid = (begin + end) / 2;
            node->left = build(begin, mid, vals);
            node->right = build(mid, end, vals);
        }
        return node;
    }
    int main() {
        std::ios::sync_with_stdio(false);
        std::vector<Segment> segs;
        std::vector<int_t> nums;
        int_t n;
        cin >> n;
        for (int_t i = 1; i <= n; i++) {
            int_t x1, y1, x2, y2;
            cin >> x1 >> y1 >> x2 >> y2;
            segs.push_back(Segment{x1, x2, y1});
            segs.push_back(Segment{x1, x2, y2});
            nums.push_back(x1);
            nums.push_back(x2);
        }
        std::sort(nums.begin(), nums.end());
        nums.resize(std::unique(nums.begin(), nums.end()) - nums.begin());
    
        std::sort(segs.begin(), segs.end(), [](const Segment& a, const Segment& b) {
            return a.height < b.height;
        });
        int_t result = 0;
        Node* root = build(0, nums.size() - 1, nums);
        const auto rank = [&](int_t x) {
            return std::lower_bound(nums.begin(), nums.end(), x) - nums.begin();
        };
    #ifdef DEBUG
        for (const auto& s : segs) {
            cout << "seg " << s.left << " " << s.right << " " << s.height << endl;
        }
    #endif
        for (int_t i = 0; i < segs.size() - 1; i++) {
            //先插入后统计结果
            const auto& curr = segs[i];
            root->insert(rank(curr.left), rank(curr.right));
            result += (segs[i + 1].height - curr.height) * root->sum;
        }
        cout << result << endl;
        return 0;
    }
    
  • CFGYM 102823 CCPC2018 桂林 B题 Array Modify

    题倒是还行,达标找规律找出来结论,然后就想到了一个多项式快速幂的做法。

    一开始写的两个log的快速幂,TLE。

    然后尝试改成exp+log,跑得看起来快了,但是遇到多组数据(n的大小递增的时候)就挂掉了?

    仔细调查后发现源于自己一直以来exp函数内一直存在的错误:

    exp内调用了log,log内调用了inv,exp内给G0(用以存储log返回值的数组)预留了upper2n(2*n)的空间,但是inv内却使用了upper2n(3*n+1)的空间!

    (调试方法:数组切换为动态分配并开启内存检查)

    对于只计算exp一次而言,这个错误不会造成影响,但是计算多次的时候由于空间污染会导致错误结果。

    解决方法:在exp内也预留upper2n(3*n)的空间,通过本题。

    另外小测试了一下,两个log的多项式快速幂跑本题极限数据,NTT变换的多项式总长度为“27262976”,而exp+log实现的快速幂,NTT变换的多项式总长度为 19922408,我原本以为因为常数原因,exp+log会跑得更慢,没想到事实并非如此。

     

    // luogu-judger-enable-o2
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    const int mod = 998244353;
    const int g = 3;
    const int LARGE = 1 << 20;
    int revs[21][LARGE + 1];
    
    int power(int base, int index);
    void transform(int* A, int len, int arg);
    void poly_inv(const int* A, int n, int* result);
    void poly_log(const int* A, int n, int* result);
    void poly_exp(const int* A, int n, int* result);
    
    std::vector transform_count;
    #ifdef TDEBUG
    #define TRANSFORM_DEBUG(n) \
        { transform_count.push_back(n); }
    #else
    #define TRANSFORM_DEBUG(n)
    #endif
    int bitRev(int base, int size2) {
        int result = 0;
        for (int i = 1; i < size2; i++) {
            result |= (base & 1);
            base >>= 1;
            result <<= 1;
        }
        result |= (base & 1);
        return result;
    }
    int upper2n(int x) {
        int result = 1;
        while (result < x)
            result *= 2;
        return result;
    }
    int main() {
        std::ios::sync_with_stdio(false);
    
        for (int i = 0; (1 << i) <= LARGE; i++) {
            for (int j = 0; j < LARGE; j++) {
                revs[i][j] = bitRev(j, i);
            }
        }
        static int F[LARGE + 1], G[LARGE + 1];
        static int arr[LARGE + 1];
        int T;
        cin >> T;
        for (int_t _i = 1; _i <= T; _i++) {
            memset(F, 0, sizeof(F));
            memset(G, 0, sizeof(G));
            memset(arr, 0, sizeof(arr));
            int_t n, L, m;
            cin >> n >> L >> m;
    
            for (int_t i = 1; i <= n; i++)
                cin >> arr[i];
            int size2 = upper2n(4 * (n + 1));
            for (int i = 0; i < size2; i++) {
                if (i < L)
                    F[i] = 1;
                else
                    F[i] = 0;
                G[i] = 0;
            }
    #ifdef DEBUG
            cout << "init = " << endl;
            cout << "F = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << F[i] << " ";
            cout << endl;
    
            cout << "G = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << G[i] << " ";
            cout << endl;
    #endif
    
            poly_log(F, n, G);
    #ifdef DEBUG
            cout << "after log = " << endl;
            cout << "G = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << G[i] << " ";
            cout << endl;
    #endif
            for (int i = 0; i < n; i++)
                G[i] = G[i] * m % mod;
            for (int i = 0; i < size2; i++)
                F[i] = 0;
    #ifdef DEBUG
            cout << "before exp = " << endl;
            cout << "G = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << G[i] << " ";
            cout << endl;
    #endif
            poly_exp(G, n, F);
    #ifdef DEBUG
            cout << "after exp = " << endl;
            cout << "F = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << F[i] << " ";
            cout << endl;
    #endif
    
            for (int i = 0; i < size2; i++) {
                if (i >= n)
                    F[i] = 0;
                if (i >= n && i < 2 * n)
                    G[i] = arr[i - n + 1];
                else
                    G[i] = 0;
            }
            std::reverse(G + n, G + 2 * n);
    
    #ifdef DEBUG
            cout << "F = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << F[i] << " ";
            cout << endl;
    
            cout << "G = ";
            for (int_t i = 0; i < 2 * n; i++)
                cout << G[i] << " ";
            cout << endl;
    #endif
    
            transform(F, size2, 1);
            transform(G, size2, 1);
            for (int i = 0; i < size2; i++)
                F[i] = (int_t)F[i] * G[i] % mod;
            transform(F, size2, -1);
            int inv = power(size2, -1);
            std::reverse(F + n, F + 2 * n);
            cout << "Case " << _i << ": ";
    
            for (int i = n; i < 2 * n; i++) {
                cout << (int_t)F[i] * inv % mod << " ";
            }
            cout << endl;
    #ifdef DEBUG
            for (int i = 0; i < n; i++) {
                cout << "F " << i << " = " << F[i] << endl;
            }
    #endif
    #ifdef TDEBUG
            for (auto x : transform_count) {
                cout << x << " ";
            }
            cout << endl;
    #endif
        }
        return 0;
    }
    
    int power(int base, int index) {
        const int phi = mod - 1;
        index = (index % phi + phi) % phi;
        base = (base % mod + mod) % mod;
        int result = 1;
        while (index) {
            if (index & 1)
                result = (int_t)result * base % mod;
            base = (int_t)base * base % mod;
            index >>= 1;
        }
        return result;
    }
    void transform(int* A, int len, int arg) {
        TRANSFORM_DEBUG(len);
        int size2 = int(log2(len) + 0.5);
        for (int i = 0; i < len; i++) {
            int x = revs[size2][i];
            if (x > i)
                std::swap(A[i], A[x]);
        }
        for (int i = 2; i <= len; i *= 2) {
            int mr = power(g, (int_t)arg * (mod - 1) / i);
            for (int j = 0; j < len; j += i) {
                int curr = 1;
                for (int k = 0; k < i / 2; k++) {
                    int u = A[j + k];
                    int t = (int_t)A[j + k + i / 2] * curr % mod;
                    A[j + k] = (u + t) % mod;
                    A[j + k + i / 2] = (u - t + mod) % mod;
                    curr = (int_t)curr * mr % mod;
                }
            }
        }
    }
    //计算多项式A在模x^n下的逆元
    // C(x)<-2B(x)-A(x)B^2(x)
    void poly_inv(const int* A, int n, int* result) {
        if (n == 1) {
            result[0] = power(A[0], -1);
            return;
        }
        poly_inv(A, n / 2 + n % 2, result);
        static int temp[LARGE + 1];
        int size2 = upper2n(3 * n + 1);
        for (int i = 0; i < size2; i++) {
            if (i < n)
                temp[i] = A[i];
            else
                temp[i] = 0;
        }
        transform(temp, size2, 1);
        transform(result, size2, 1);
        for (int i = 0; i < size2; i++) {
            result[i] =
                ((int_t)2 * result[i] % mod -
                 (int_t)temp[i] * result[i] % mod * result[i] % mod + 2 * mod) %
                mod;
        }
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++) {
            if (i < n) {
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            } else {
                result[i] = 0;
            }
        }
    }
    void poly_log(const int* A, int n, int* result) {
        int size2 = upper2n(3 * n + 1);
        static int Ad[LARGE + 1];
        for (int i = 0; i < size2; i++) {
            if (i < n) {
                Ad[i] = (int_t)(i + 1) * A[i + 1] % mod;
            } else {
                Ad[i] = 0;
            }
            result[i] = 0;
        }
        transform(Ad, size2, 1);
        poly_inv(A, n, result);
        transform(result, size2, 1);
        for (int i = 0; i < size2; i++) {
            result[i] = (int_t)result[i] * Ad[i] % mod;
        }
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++)
            if (i < n)
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            else
                result[i] = 0;
        for (int i = n - 1; i >= 1; i--) {
            result[i] = (int_t)result[i - 1] * power(i, -1) % mod;
        }
        result[0] = 0;
    }
    void poly_exp(const int* A, int n, int* result) {
        if (n == 1) {
            assert(A[0] == 0);
            result[0] = 1;
            return;
        }
        poly_exp(A, n / 2 + n % 2, result);
        static int Ax[LARGE + 1];
        // memset(G0, 0, sizeof(G0));
        // memset(Ax,0,sizeof(Ax));
        int size2 = upper2n(3 * n + 1);
        // int* G0 = new int[size2];
        static int G0[LARGE + 1];
        for (int_t i = 0; i < size2; i++) {
            if (i < n)
                Ax[i] = A[i];
            else
                Ax[i] = 0;
            G0[i] = 0;
        }
        poly_log(result, n, G0);
    
        transform(Ax, size2, 1);
        transform(G0, size2, 1);
        transform(result, size2, 1);
        for (int i = 0; i < size2; i++)
            result[i] =
                (result[i] - (int_t)result[i] * (G0[i] - Ax[i] + mod) % mod + mod) %
                mod;
        transform(result, size2, -1);
        for (int i = 0; i < size2; i++)
            if (i < n)
                result[i] = (int_t)result[i] * power(size2, -1) % mod;
            else
                result[i] = 0;
        // delete[] G0;
    }
    
  • CF1009F Dominant Indices 长剖入门题 复习

    复习长链剖分。

    首先划分长链,走到每个长链顶就给这个长链开个数组拿来记答案,数组的长度(从0开始算)是整个长链上点的个数,下标为k的位置所表示的意义是到长链顶距离为k的点的个数,沿着重链走的时候,直接指针偏移一位就行了(毕竟状态在深度上连续)。

    然后DFS,走到每一个点,先给这个点的DP数组指针(设为arr),下标为0的位置加个1,表示这个点自己的答案。然后如果这个点有重子节点,就先沿着重子节点走,走的时候DP数组直接偏移上1,然后这个点的答案先记录为重子节点的答案+1(我们要考虑所有子节点的答案,一会拿这玩意来更新)。

    对于非重子节点的点,开个新的DP数组,然后直接传下去往下递归。

    非重子节点返回后尝试更新答案。枚举非重子节点DP数组上的每一个状态,把他的值加到当前节点的DP值上(即合并子链状态,把状态合并到重链上),然后顺带更新一波结果(如果某个深度的值出现的比当前的多,就直接更新,如果相当就比一下深度),合并晚状态后直接把非重子节点的内存回收掉即可。

    显然复杂度是$\theta(n)$的,每条重链会有一个DP数组,这个DP数组的长度是重链的长度,每条重链的答案只会被合并一次,而所有重链的长度之和是$n$。

    #include <algorithm>
    #include <cstdio>
    #include <cstring>
    #include <iostream>
    #include <vector>
    using int_t = long long int;
    
    const int_t LARGE = 1e6;
    
    using std::cin;
    using std::cout;
    using std::endl;
    
    int_t results[LARGE + 1];
    
    int_t depth[LARGE + 1];
    int_t maxDepth[LARGE + 1];
    int_t maxChd[LARGE + 1];
    
    std::vector<int_t> graph[LARGE + 1];
    int_t n;
    
    void DFS1(int_t vtx, int_t from = -1) {
        maxChd[vtx] = -1;
        maxDepth[vtx] = depth[vtx];
        for (auto to : graph[vtx])
            if (to != from) {
                depth[to] = depth[vtx] + 1;
                DFS1(to, vtx);
                if (maxChd[vtx] == -1 || maxDepth[to] > maxDepth[maxChd[vtx]]) {
                    maxChd[vtx] = to;
                }
                maxDepth[vtx] = std::max(maxDepth[vtx], maxDepth[to]);
            }
    }
    // arr 距离为x的点的个数
    void DFS2(int_t vtx, int_t from = -1, int_t* arr = nullptr) {
        arr[0]++;
        auto& result = results[vtx];
        if (maxChd[vtx] != -1) {
            DFS2(maxChd[vtx], vtx, arr + 1);
            result = results[maxChd[vtx]] + 1;  //答案初始值
        }
        for (auto to : graph[vtx])
            if (to != from && to != maxChd[vtx]) {
                int_t arrlen = maxDepth[to] - depth[vtx] + 1;
                int_t* next = new int_t[arrlen];
                std::fill(next, next + arrlen, 0);
                // next[0] = 1;
                DFS2(to, vtx, next + 1);
                for (int_t i = 1; i < arrlen; i++) {
                    arr[i] += next[i];
                    if (arr[i] > arr[result])
                        result = i;
                    else if (arr[i] == arr[result] && i < result)
                        result = i;
                }
                delete[] next;
            }
        if (arr[result] == 1)
            result = 0;  //最小的距离
    }
    int main() {
        std::ios::sync_with_stdio(false);
        cin >> n;
        for (int_t i = 1; i <= n - 1; i++) {
            int_t u, v;
            cin >> u >> v;
            graph[u].push_back(v);
            graph[v].push_back(u);
        }
        depth[1] = 1;
        DFS1(1);
    #ifdef DEBUG
    
        for (int_t i = 1; i <= n; i++) {
            cout << "maxchd " << i << " = " << maxChd[i] << endl;
        }
    #endif
        static int_t arr[LARGE + 1];
        DFS2(1, -1, arr);
        for (int_t i = 1; i <= n; i++) {
            cout << results[i] << endl;
        }
    
        return 0;
    }

     

  • CCPC2019 秦皇岛 A Angle Beats

    傻逼HDU

    原题内存1G,贵题内存128M,真的nb

    傻逼题,先对于一个点i,枚举所有点j,然后计算i和j构成的直线的斜率,记在i的map里(即对于每个点i,统计i作为线段一个端点的情况下,不同的斜率取值所对应的点的个数)

    然后对于每个询问,分$A_i$是直角端点和非直角端点的情况讨论。

    $A_i$是直角端点的情况下,枚举下n个点,算出来斜率的取值情况然后存起来。然后再枚举一个点,钦定这是第一条直角边,然后算出来另一条直角边的斜率,去刚刚存的map里反查有多少种情况,最后除个2(一种情况会被算两次)

    $A_i$不是直角端点的情况下,枚举一个点,连起来构成一条直角边,然后算一下另一条直角边的斜率,然后去一开始预处理的数组里去查下能有多少种即可。(这个不需要除2,因为不会算重)

    如何存斜率?别用浮点数,写个有理数类。如何记斜率不存在?分母写0,分子写1(分子不要写别的数,不方便统计)。如何记斜率为0?分子写1,分母写0。

    另外为了统计方便,我们所存储的有理数均为最简分数,并且正负号全都在分子上(0和不存在除外,这时不需要考虑符号)

    要做直接去CF交吧,HDU太傻逼了。

    要用unordered_map并且自己写个hasher,map太慢了。

    https://codeforces.com/gym/102361/problem/A

    #include <assert.h>
    #include <functional>
    #include <iostream>
    #include <map>
    #include <unordered_map>
    using int_t = long long int;
    
    using std::cin;
    using std::cout;
    using std::endl;
    
    int_t gcd(int_t a, int_t b) {
        a = std::abs(a);
        b = std::abs(b);
        if (b == 0)
            return a;
        return gcd(b, a % b);
    }
    
    struct R {
        //有理数a/b
        int_t a = 0, b = 0;
        R simplify() const {
            int_t d = gcd(a, b);
            return R{a / d, b / d};
        }
        // R operator*(const R& rhs) const {
        //     if (b == 0 || rhs.b == 0)
        //         return R{1, 0};
        //     return R{a * rhs.a, b * rhs.b}.simplify();
        // }
        bool operator<(const R& rhs) const {
            return (long double)a / b < (long double)rhs.a / rhs.b;
        }
        bool operator==(const R& rhs) const { return a == rhs.a && b == rhs.b; }
        R(int_t a = 0, int_t b = 0) : a(a), b(b) {}
    };
    
    const auto my_hash = [](const R& x) -> size_t {
        return (((uint64_t)std::abs(x.a) << 31) | (uint64_t)std::abs(x.b)) %
               998244353;
    };
    
    struct my_hasher {
        size_t operator()(const R& x) const { return my_hash(x); }
    };
    //以某个点为一个端点的,斜率为R的点的个数
    //无序点对
    std::unordered_map<R, int, my_hasher> byS[2001];
    
    struct P {
        int x, y;
    } points[2001];
    std::ostream& operator<<(std::ostream& os, const R& r) {
        os << "R{" << r.a << "," << r.b << "}";
        return os;
    }
    int n, q;
    int main() {
        {
            (scanf("%d%d", &n, &q));
            // break;
            for (int i = 1; i <= n; i++) {
                auto& ref = points[i];
                scanf("%d%d", &ref.x, &ref.y);
                // byS[i].clear();
            }
            for (int i = 1; i <= n; i++) {
    #ifdef DEBUG
                cout << "pnt " << i << " = " << points[i].x << " " << points[i].y
                     << endl;
    #endif
                for (int j = 1; j <= n; j++) {
                    if (i == j)
                        continue;
                    int_t dy = points[j].y - points[i].y;
                    int_t dx = points[j].x - points[i].x;
    #ifdef DEBUG
                    cout << "point " << i << " " << j << " dy = " << dy
                         << " dx = " << dx << endl;
    
    #endif
                    if (dx < 0) {
                        dy *= -1, dx *= -1;
                    }
                    int_t d = gcd(dy, dx);
                    dy /= d, dx /= d;
                    if (dx == 0)
                        dy = 1;
                    else if (dy == 0)
                        dx = 1;
                    byS[i][R{dy, dx}]++;
                }
            }
    
            while (q--) {
                int x, y;
                scanf("%d%d", &x, &y);
                int_t result = 0;
    #ifdef DEBUG
                cout << "BEGIN " << x << " " << y << endl;
    #endif
                //作为直角顶点 先预处理
                {
                    int_t curr = 0;
                    static std::unordered_map<R, int, my_hasher> byS;
                    // for(int i=1;i<=n;i++) byS[i].clear();
                    byS.clear();
                    for (int i = 1; i <= n; i++) {
                        int_t dy = y - points[i].y;
                        int_t dx = x - points[i].x;
                        if (dx < 0) {
                            dy *= -1, dx *= -1;
                        }
                        int_t d = gcd(dy, dx);
                        dy /= d, dx /= d;
                        if (dx == 0)
                            dy = 1;
                        else if (dy == 0)
                            dx = 1;
                        byS[R{dy, dx}]++;
                    }
    #ifdef DEBUG
                    for (const auto& kvp : byS) {
                        cout << kvp.first << " = " << kvp.second << endl;
                    }
    #endif
    
    #ifdef DEBUG
    
                    cout << "TYPE 2" << endl;
    #endif
                    //枚举一条直角边 然后根据斜率查另一条
                    for (int i = 1; i <= n; i++) {
    #ifdef DEBUG
                        cout << "FOR PNT " << i << endl;
    
    #endif
                        int_t dy = y - points[i].y;
                        int_t dx = x - points[i].x;
                        dy *= -1;
                        if (dy < 0) {
                            dy *= -1, dx *= -1;
                        }
                        int_t d = gcd(dy, dx);
                        dy /= d, dx /= d;
    
                        R r0{dx, dy};
                        if (r0.b == 0) {
                            r0.a = 1;
                        } else if (r0.a == 0)
                            r0.b = 1;
                        if (byS.count(r0))
                            curr += byS[r0];
                    }
    #ifdef DEBUG
                    cout << "type1 curr=" << curr << " div2 " << curr / 2 << endl;
    #endif
                    // assert(curr % 2 == 0);
                    result += curr / 2;
                }
                //作为非直角顶点 枚举一个点,连上一条直角边,然后去查斜率
                {
                    int_t curr = 0;
                    for (int i = 1; i <= n; i++) {
                        int_t dy = y - points[i].y;
                        int_t dx = x - points[i].x;
                        dy *= -1;
                        if (dy < 0) {
                            dy *= -1, dx *= -1;
                        }
                        int_t d = gcd(dy, dx);
                        dx /= d, dy /= d;
                        R r0{dx, dy};
                        if (r0.b == 0)
                            r0.a = 1;
                        else if (r0.a == 0)
                            r0.b = 1;
                        if (byS[i].count(r0))
                            curr += byS[i][r0];
    #ifdef DEBUG
    
                        cout << "type 2 at pnt " << i << " req r = " << r0
                             << " count " << byS[i][r0] << endl;
    #endif
                    }
                    // assert(curr % 2 == 0);
    #ifdef DEBUG
                    cout << "type2 curr=" << curr << " div2 " << curr / 2 << endl;
    #endif
                    result += curr;
                }
                printf("%lld\n", result);
            }
        }
        return 0;
    }

     

  • 洛谷3327 约数个数和 复习

    众所周知,这题不是个人也会。

    $$ \text{首先可知}d\left( nm \right) =\sum_{a|n}{\sum_{b|m}{\left[ \mathrm{gcd}\left( a,b \right) =1 \right]}} \\ \sum_{1\le i\le n}{\sum_{1\le j\le m}{\sum_{a|i}{\sum_{b|j}{\left[ \mathrm{gcd}\left( a,b \right) =1 \right]}}}} \\ =\sum_{1\le i\le n}{\sum_{1\le j\le m}{\sum_{a|i}{\sum_{b|j}{\sum_{k|gcd\left( a,b \right)}{\mu \left( k \right)}}}}} \\ =\sum_{1\le k\le n}{\mu \left( k \right) \sum_{1\le i\le n}{\sum_{1\le j\le m}{\sum_{a|i}{\sum_{b|j}{\left[ k|\mathrm{gcd}\left( a,b \right) \right]}}}}} \\ =\sum_{1\le k\le n}{\mu \left( k \right) \sum_{1\le a\le n}{\sum_{1\le b\le m}{\sum_{a|i,i\le n}{\sum_{b|j,j\le m}{\left[ k|gcd\left( a,b \right) \right]}}}}} \\ a\gets ak,b\gets bk \\ =\sum_{1\le k\le n}{\mu \left( k \right) \sum_{1\le ak\le n}{\sum_{1\le bk\le m}{\left[ k|gcd\left( ak,bk \right) \right] \lfloor \frac{n}{ak} \rfloor \lfloor \frac{m}{bk} \rfloor}}} \\ =\sum_{1\le k\le n}{\mu \left( k \right) \sum_{1\le a\le \lfloor \frac{n}{k} \rfloor}{\sum_{1\le b\le \lfloor \frac{m}{k} \rfloor}{\lfloor \frac{n}{ak} \rfloor \lfloor \frac{m}{bk} \rfloor}}} \\ =\sum_{1\le k\le n}{\mu \left( k \right) \sum_{1\le a\le \lfloor \frac{n}{k} \rfloor}{\sum_{1\le b\le \lfloor \frac{m}{k} \rfloor}{\lfloor \frac{\lfloor \frac{n}{k} \rfloor}{a} \rfloor \lfloor \frac{\lfloor \frac{m}{k} \rfloor}{b} \rfloor}}} \\ =\sum_{1\le k\le n}{\mu \left( k \right) \left( \sum_{1\le a\le \lfloor \frac{n}{k} \rfloor}{\lfloor \frac{\lfloor \frac{n}{k} \rfloor}{a} \rfloor} \right) \left( \sum_{1\le b\le \lfloor \frac{m}{k} \rfloor}{\lfloor \frac{\lfloor \frac{m}{k} \rfloor}{b} \rfloor} \right)} \\ S\left( n \right) =\sum_{1\le i\le n}{\lfloor \frac{n}{i} \rfloor} \\ \text{原式}=\sum_{1\le k\le n}{\mu \left( k \right) S\left( \lfloor \frac{n}{k} \rfloor \right) S\left( \lfloor \frac{m}{k} \rfloor \right)} \\ $$

    #include <algorithm>
    #include <cstdio>
    #include <cstring>
    #include <iostream>
    using int_t = long long int;
    
    const int_t LARGE = 5e4;
    
    using std::cin;
    using std::cout;
    using std::endl;
    
    int_t S[LARGE + 1];
    bool isPrime[LARGE + 1];
    int_t factor[LARGE + 1];
    int_t mu[LARGE + 1];
    int_t muSum[LARGE + 1];
    
    int_t Sx(int_t n) {
        int_t result = 0;
        int_t i = 1;
        while (i <= n) {
            int_t next = n / (n / i);
            result += (next - i + 1) * (n / i);
            i = next + 1;
        }
        return result;
    }
    
    int_t query(int_t n, int_t m) {
        if (n > m)
            std::swap(n, m);
        int_t result = 0;
        int_t i = 1;
        while (i <= n) {
            int_t next = std::min(n / (n / i), m / (m / i));
            result += (muSum[next] - muSum[i - 1]) * S[n / i] * S[m / i];
            i = next + 1;
        }
        return result;
    }
    
    int main() {
        for (int_t i = 1; i <= LARGE; i++) {
            S[i] = Sx(i);
        }
        memset(isPrime, 1, sizeof(isPrime));
        for (int_t i = 2; i <= LARGE; i++) {
            if (isPrime[i]) {
                factor[i] = i;
                for (int_t j = i * i; j <= LARGE; j += i) {
                    isPrime[j] = false;
                    if (!factor[j])
                        factor[j] = i;
                }
            }
        }
        mu[1] = muSum[1] = 1;
        for (int_t i = 2; i <= LARGE; i++) {
            int_t factor = ::factor[i];
            if (i / factor % factor == 0)
                mu[i] = 0;
            else
                mu[i] = mu[i / factor] * -1;
            muSum[i] = muSum[i - 1] + mu[i];
        }
        int T;
        scanf("%d", &T);
        while (T--) {
            int n, m;
            scanf("%d%d", &n, &m);
            auto result = query(n, m);
            printf("%lld\n", result);
        }
        return 0;
    }