Generating Functions - OGF

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// NTT for mod 998244353 (primitive root 3)
5
struct NTT {
6
    static constexpr uint32_t MOD = 998244353;
7
    static constexpr uint32_t G = 3;
8
9
    uint32_t mod_pow(uint32_t a, uint32_t e) const {
10
        uint64_t r = 1, x = a;
11
        while (e) {
12
            if (e & 1) r = r * x % MOD;
13
            x = x * x % MOD;
14
            e >>= 1;
15
        }
16
        return (uint32_t)r;
17
    }
18
19
    void ntt(vector<uint32_t> &a, bool invert) const {
20
        int n = (int)a.size();
21
        static vector<int> rev;
22
        static vector<uint32_t> roots{0, 1};
23
        if ((int)rev.size() != n) {
24
            int k = __builtin_ctz(n);
25
            rev.assign(n, 0);
26
            for (int i = 0; i < n; i++) {
27
                rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (k - 1));
28
            }
29
        }
30
        if ((int)roots.size() < n) {
31
            int k = __builtin_ctz(roots.size());
32
            roots.resize(n);
33
            while ((1 << k) < n) {
34
                uint32_t z = mod_pow(G, (MOD - 1) >> (k + 1));
35
                for (int i = 1 << (k - 1); i < (1 << k); i++) {
36
                    roots[2 * i] = roots[i];
37
                    roots[2 * i + 1] = (uint64_t)roots[i] * z % MOD;
38
                }
39
                k++;
40
            }
41
        }
42
        for (int i = 0; i < n; i++) if (i < rev[i]) swap(a[i], a[rev[i]]);
43
        for (int len = 1; len < n; len <<= 1) {
44
            for (int i = 0; i < n; i += 2 * len) {
45
                for (int j = 0; j < len; j++) {
46
                    uint32_t z = roots[len + j];
47
                    uint32_t u = a[i + j];
48
                    uint32_t v = (uint64_t)a[i + j + len] * z % MOD;
49
                    uint32_t x = u + v;
50
                    a[i + j] = (x >= MOD ? x - MOD : x);
51
                    a[i + j + len] = (u + MOD - v);
52
                    if (a[i + j + len] >= MOD) a[i + j + len] -= MOD;
53
                }
54
            }
55
        }
56
        if (invert) {
57
            reverse(a.begin() + 1, a.end());
58
            uint32_t inv_n = mod_pow(n, MOD - 2);
59
            for (int i = 0; i < n; i++) a[i] = (uint64_t)a[i] * inv_n % MOD;
60
        }
61
    }
62
63
    vector<uint32_t> multiply(vector<uint32_t> a, vector<uint32_t> b, int need = -1) const {
64
        if (a.empty() || b.empty()) return {};
65
        int n = 1, as = (int)a.size(), bs = (int)b.size();
66
        int tot = as + bs - 1;
67
        if (need != -1) tot = min(tot, need);
68
        while (n < as + bs - 1) n <<= 1;
69
        a.resize(n); b.resize(n);
70
        ntt(a, false); ntt(b, false);
71
        for (int i = 0; i < n; i++) a[i] = (uint64_t)a[i] * b[i] % MOD;
72
        ntt(a, true);
73
        a.resize(as + bs - 1);
74
        if (need != -1 && (int)a.size() > need) a.resize(need);
75
        return a;
76
    }
77
} ntt;
78
79
// Multiply a list of polynomials using divide-and-conquer; truncate to degree <= max_deg
80
vector<uint32_t> multiply_all(vector<vector<uint32_t>> &polys, int l, int r, int max_deg) {
81
    if (l + 1 == r) {
82
        auto res = polys[l];
83
        if ((int)res.size() > max_deg + 1) res.resize(max_deg + 1);
84
        return res;
85
    }
86
    int m = (l + r) / 2;
87
    auto L = multiply_all(polys, l, m, max_deg);
88
    auto R = multiply_all(polys, m, r, max_deg);
89
    auto P = ntt.multiply(L, R, max_deg + 1);
90
    if ((int)P.size() > max_deg + 1) P.resize(max_deg + 1);
91
    return P;
92
}
93
94
int main() {
95
    ios::sync_with_stdio(false);
96
    cin.tie(nullptr);
97
98
    // Problem: Given coin denominations d_1..d_m (unbounded), count the number of multisets summing to S.
99
    // OGF: F(x) = prod_{i=1..m} (1 + x^{d_i} + x^{2 d_i} + ...), answer is [x^S] F(x).
100
101
    int m, S;
102
    cin >> m >> S;
103
    vector<int> d(m);
104
    for (int i = 0; i < m; i++) cin >> d[i];
105
106
    vector<vector<uint32_t>> polys;
107
    polys.reserve(m);
108
    for (int i = 0; i < m; i++) {
109
        vector<uint32_t> p(S + 1, 0);
110
        for (int e = 0; e <= S; e += d[i]) p[e] = 1; // 1 + x^{d} + x^{2d} + ... up to degree S
111
        polys.push_back(move(p));
112
    }
113
114
    auto res = multiply_all(polys, 0, (int)polys.size(), S);
115
    uint32_t ans = (S < (int)res.size() ? res[S] : 0);
116
    cout << ans << "\n";
117
    return 0;
118
}
119

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Compute nth term of k-term linear recurrence:
5
// a_n = c1*a_{n-1} + c2*a_{n-2} + ... + ck*a_{n-k}
6
// given initial terms a_0..a_{k-1}. Uses Kitamasa O(k^2 log n).
7
8
using int64 = long long;
9
const int64 MOD = 1000000007LL;
10
11
vector<int64> combine(const vector<int64>& a, const vector<int64>& b, const vector<int64>& coef) {
12
    int k = (int)coef.size();
13
    vector<int64> t(2*k - 1, 0);
14
    for (int i = 0; i < k; i++)
15
        for (int j = 0; j < k; j++)
16
            t[i + j] = (t[i + j] + a[i] * b[j]) % MOD;
17
    for (int i = 2*k - 2; i >= k; i--) {
18
        for (int j = 1; j <= k; j++) {
19
            t[i - j] = (t[i - j] + t[i] * coef[j - 1]) % MOD;
20
        }
21
    }
22
    t.resize(k);
23
    return t;
24
}
25
26
int64 nth_term(vector<int64> init, vector<int64> coef, long long n) {
27
    int k = (int)coef.size();
28
    if (n < (long long)init.size()) return (init[n] % MOD + MOD) % MOD;
29
    vector<int64> pol(k), e(pol); // pol encodes current linear combination, e is like x^1
30
    pol[0] = 1; // x^0
31
    e[1 % k] = 1; // x^1
32
    long long m = n;
33
    while (m) {
34
        if (m & 1) pol = combine(pol, e, coef);
35
        e = combine(e, e, coef);
36
        m >>= 1;
37
    }
38
    int64 ans = 0;
39
    for (int i = 0; i < k; i++) ans = (ans + pol[i] * init[i]) % MOD;
40
    return (ans + MOD) % MOD;
41
}
42
43
int main() {
44
    ios::sync_with_stdio(false);
45
    cin.tie(nullptr);
46
47
    // Example: Fibonacci F_0=0, F_1=1, F_n = F_{n-1} + F_{n-2}
48
    // OGF A(x) = x / (1 - x - x^2), rational -> fast nth computation.
49
50
    long long n; cin >> n;
51
    vector<int64> init = {0, 1};
52
    vector<int64> coef = {1, 1}; // c1=1, c2=1
53
    cout << nth_term(init, coef, n) << "\n";
54
    return 0;
55
}
56

1
#include <bits/stdc++.h>
2
using namespace std;
3
using int64 = long long;
4
const int MOD = 1000000007;
5
6
int addmod(int a, int b){ a+=b; if(a>=MOD) a-=MOD; return a; }
7
int submod(int a, int b){ a-=b; if(a<0) a+=MOD; return a; }
8
int mulmod(long long a, long long b){ return int((a*b)%MOD); }
9
int modpow(int a, long long e){ long long r=1,x=a; while(e){ if(e&1) r=r*x%MOD; x=x*x%MOD; e>>=1; } return (int)r; }
10
11
int main(){
12
    ios::sync_with_stdio(false);
13
    cin.tie(nullptr);
14
15
    // Using [x^n] (1 - x)^{-k} = C(n + k - 1, n) (stars and bars)
16
    int T; cin >> T;
17
    vector<pair<int,int>> qs(T);
18
    int maxn = 0, maxk = 0;
19
    for(int i=0;i<T;i++){ cin >> qs[i].first >> qs[i].second; maxn=max(maxn, qs[i].first + qs[i].second - 1); }
20
    // Precompute factorials up to maxN = max(n + k - 1)
21
    int N = maxn;
22
    vector<int> fact(N+1), ifact(N+1);
23
    fact[0]=1; for(int i=1;i<=N;i++) fact[i]=mulmod(fact[i-1], i);
24
    ifact[N]=modpow(fact[N], MOD-2);
25
    for(int i=N;i>0;i--) ifact[i-1]=mulmod(ifact[i], i);
26
27
    auto C = [&](int n, int r)->int{
28
        if(r<0||r>n) return 0;
29
        return mulmod(fact[n], mulmod(ifact[r], ifact[n-r]));
30
    };
31
32
    for(auto [n,k]: qs){
33
        // answer = C(n + k - 1, n)
34
        cout << C(n + k - 1, n) << "\n";
35
    }
36
    return 0;
37
}
38

1
#include <bits/stdc++.h>
2
using namespace std;
3
using int64 = long long;
4
const int MOD = 1000000007;
5
6
int addmod(int a, int b){ a+=b; if(a>=MOD) a-=MOD; return a; }
7
int submod(int a, int b){ a-=b; if(a<0) a+=MOD; return a; }
8
int mulmod(long long a, long long b){ return int((a*b)%MOD); }
9
int modpow(int a, long long e){ long long r=1,x=a; while(e){ if(e&1) r=r*x%MOD; x=x*x%MOD; e>>=1; } return (int)r; }
10
11
int main(){
12
    ios::sync_with_stdio(false);
13
    cin.tie(nullptr);
14
15
    // Given P and Q with Q(0)=1, compute A(x)=P(x)/Q(x) up to x^N coefficients.
16
    // Coefficient relation: sum_{i=0}^{degQ} q_i a_{n-i} = p_n (with q_0=1), so
17
    // a_n = p_n - sum_{i=1}^{min(n,degQ)} q_i * a_{n-i}.
18
19
    int degP, degQ, N; // degrees and how many coefficients to output (0..N)
20
    cin >> degP >> degQ >> N;
21
    vector<int> P(degP+1), Q(degQ+1);
22
    for (int i = 0; i <= degP; i++) cin >> P[i];
23
    for (int i = 0; i <= degQ; i++) cin >> Q[i];
24
25
    if (Q[0] != 1) {
26
        // Normalize to make Q[0]=1: divide P and Q by Q[0]
27
        int invq0 = modpow(Q[0], MOD-2);
28
        for (auto &x : P) x = mulmod(x, invq0);
29
        for (auto &x : Q) x = mulmod(x, invq0);
30
    }
31
32
    vector<int> A(N+1, 0);
33
    for (int n = 0; n <= N; n++) {
34
        long long val = (n <= degP ? P[n] : 0);
35
        int lim = min(n, degQ);
36
        for (int i = 1; i <= lim; i++) {
37
            val -= 1LL * Q[i] * A[n - i] % MOD;
38
            if (val < 0) val += MOD;
39
        }
40
        A[n] = (int)(val % MOD);
41
    }
42
43
    // Example input to test quickly:
44
    // For Fibonacci, A(x) = x / (1 - x - x^2): degP=1, P=[0,1]; degQ=2, Q=[1, (MOD-1+1)%MOD, (MOD-1+1)%MOD] -> Q=[1,MOD-1,MOD-1] for -x - x^2.
45
46
    for (int i = 0; i <= N; i++) cout << A[i] << (i==N?'\n':' ');
47
    return 0;
48
}
49