loj 3090 「BJOI2019」勘破神機 - 數學
阿新 • • 發佈：2019-04-27
二次 gcd 組合數 ole long long 。。 borde pac color
G‘_{n} &= 3G_{n - 1} + 2\sum_{i = 0}^{n - 2}G_{i} \\
&= 4G_{n - 1} - 3G_{n - 2} - 2\sum_{i = 0}^{n - 3} G_{i - 1} + 2\sum_{i =0}^{n - 2}G_i \\
&= 4G_{n - 1} - G_{n - 2}
\end{align}
$$
題目傳送門
傳送門
題目大意
設$F_{n}$表示用$1\times 2$的骨牌填$2\times n$的網格的方案數，設$G_{n}$$表示用$1\times 2$的骨牌填$3\times n$的網格的方案數.
- 給定$l, r, k$，求$\frac{1}{r - l + 1}\sum_{i = l}^{r} \binom{F_{i}}{k}$.
- 給定$l, r, k$，求$\frac{1}{r - l + 1}\sum_{i = l}^{r} \binom{G_{i}}{k}$.
之前好像在loj / uoj群上問過求Fibonacci求立方和。sad。。。
校內考的時候忘了求數列通項，sad。。。
首先用Stirling數把組合數拆成通常冪。
對於第一部分，直接用$Fibonacci$通項，由於$5$不是模$998244353$的二次剩余，所以把答案在運算中表示成$a + b\sqrt{5}$的形式。
第$n$項的$k$次冪大概是$(a^n + b^n)^{K}$，這個方便直接求和。用二項式定理展開就行了。
對於第二部分。考察選手是否上過高考數學課（bushi
顯然當$2 \nmid n$時$G_{n} = 0$，現在考慮$G‘_{n} = G_{2n}$。
註意到最後只有這幾種填法：
所以有
$$
G‘_{n} &= 3G_{n - 1} + 2\sum_{i = 0}^{n - 2}G_{i} \\
&= 4G_{n - 1} - 3G_{n - 2} - 2\sum_{i = 0}^{n - 3} G_{i - 1} + 2\sum_{i =0}^{n - 2}G_i \\
&= 4G_{n - 1} - G_{n - 2}
\end{align}
$$
用特征根法推通項（不會的話可以回教室了）。然後做完了。
（我被Jode錘爆了。瑟瑟發抖.jpg）
Code
#include <bits/stdc++.h>
using namespace std;
typedef bool boolean;
#define ll long long
void exgcd(int a, int b, int& x, int& y) {
if (!b) {
x = 1, y = 0;
} else {
exgcd(b, a % b, y, x);
y -= (a / b) * x;
}
}
int inv(int a, int n) {
int x, y;
exgcd(a, n, x, y);
return (x < 0) ? (x + n) : (x);
}
const int Mod = 998244353;
template <const int Mod = :: Mod>
class Z {
public:
int v;
Z() : v(0) { }
Z(int x) : v(x){ }
Z(ll x) : v(x % Mod) { }
friend Z operator + (const Z& a, const Z& b) {
int x;
return Z(((x = a.v + b.v) >= Mod) ? (x - Mod) : (x));
}
friend Z operator - (const Z& a, const Z& b) {
int x;
return Z(((x = a.v - b.v) < 0) ? (x + Mod) : (x));
}
friend Z operator * (const Z& a, const Z& b) {
return Z(a.v * 1ll * b.v);
}
friend Z operator ~(const Z& a) {
return inv(a.v, Mod);
}
friend Z operator - (const Z& a) {
return Z(0) - a;
}
Z& operator += (Z b) {
return *this = *this + b;
}
Z& operator -= (Z b) {
return *this = *this - b;
}
Z& operator *= (Z b) {
return *this = *this * b;
}
friend boolean operator == (const Z& a, const Z& b) {
return a.v == b.v;
}
};
Z<> qpow(Z<> a, int p) {
Z<> rt = Z<>(1), pa = a;
for ( ; p; p >>= 1, pa = pa * pa) {
if (p & 1) {
rt = rt * pa;
}
}
return rt;
}
typedef Z<> Zi;
template <const int I>
class ComplexTemp {
public:
Zi r, v;
ComplexTemp() : r(0), v(0) { }
ComplexTemp(Zi r) : r(r), v(0) { }
ComplexTemp(Zi r, Zi v) : r(r), v(v) { }
friend ComplexTemp operator + (const ComplexTemp& a, const ComplexTemp& b) {
return ComplexTemp(a.r + b.r, a.v + b.v);
}
friend ComplexTemp operator - (const ComplexTemp& a, const ComplexTemp& b) {
return ComplexTemp(a.r - b.r, a.v - b.v);
}
friend ComplexTemp operator - (const ComplexTemp& a, const int& b) {
return ComplexTemp(a.r - b, a.v);
}
friend ComplexTemp operator * (const ComplexTemp& a, const ComplexTemp& b) {
return ComplexTemp(a.r * b.r + a.v * b.v * I, a.r * b.v + a.v * b.r);
}
friend ComplexTemp operator * (const ComplexTemp& a, const Zi& x) {
return ComplexTemp(a.r * x, a.v * x);
}
friend ComplexTemp operator / (const ComplexTemp& a, const ComplexTemp& b) {
ComplexTemp c = b.conj();
return a * c * ~((b * c).r);
}
ComplexTemp conj() const {
return ComplexTemp(r, -v);
}
boolean operator == (ComplexTemp b) {
return r == b.r && v == b.v;
}
};
const int Kmx = 510;
Zi s1[Kmx][Kmx];
Zi comb[Kmx][Kmx];
Zi fac[Kmx], _fac[Kmx];
void prepare(int n) {
fac[0] = 1;
for (int i = 1; i <= n; i++) {
fac[i] = fac[i - 1] * i;
}
_fac[n] = ~fac[n];
for (int i = n; i; i--) {
_fac[i - 1] = _fac[i] * i;
}
s1[0][0] = 1;
for (int i = 1; i <= n; i++) {
s1[i][0] = 0, s1[i][i] = 1;
for (int j = 1; j < i; j++) {
s1[i][j] = s1[i - 1][j - 1] + s1[i - 1][j] * (i - 1);
}
}
comb[0][0] = 1;
for (int i = 1; i <= n; i++) {
comb[i][0] = comb[i][i] = 1;
for (int j = 1; j < i; j++) {
comb[i][j] = comb[i - 1][j - 1] + comb[i - 1][j];
}
}
}
template <typename T>
T qpow(T a, ll p) {
T rt = Zi(1), pa = a;
for ( ; p; p >>= 1, pa = pa * pa) {
if (p & 1) {
rt = rt * pa;
}
}
return rt;
}
namespace subtask1 {
typedef ComplexTemp<5> Complex;
const Zi inv2 ((Mod + 1) >> 1);
const Complex q1 (inv2, inv2), q2 (inv2, -inv2);
Complex pwq1[Kmx], pwq2[Kmx], pwqq[Kmx][Kmx];
Complex get_sum(int k1, int k2, ll n) {
Complex x = pwq1[k1] * pwq2[k2];
if (x == Zi(1))
return Zi(n);
return (pwqq[k1][k2] - 1) / (x - 1);
}
inline void init() {
pwq1[0] = pwq2[0] = Zi(1);
for (int i = 1; i < Kmx; i++) {
pwq1[i] = pwq1[i - 1] * q1;
pwq2[i] = pwq2[i - 1] * q2;
}
}
Zi work(ll n, int K) {
Complex coef = qpow(Complex(0, ~Zi(5)), K);
Complex ret (0, 0);
for (int k = 0; k <= K; k++) {
// Complex tmp = get_sum(pwq1[k] * pwq2[K - k], n) * comb[K][k];
Complex tmp = get_sum(k, K - k, n) * comb[K][k];
if ((K - k) & 1) {
ret = ret - tmp;
} else {
ret = ret + tmp;
}
}
ret = ret * coef;
assert(ret.v.v == 0);
// cerr << n << " " << K << " " << ret.r.v << ‘\n‘;
return ret.r;
}
Zi solve(ll n, int K) {
Zi ans = 0;
Complex q1n = qpow(q1, n + 1);
Complex q2n = qpow(q2, n + 1);
pwqq[0][0] = Zi(1);
for (int i = 0; i <= K; i++) {
for (int j = (i == 0); i + j <= K; j++) {
pwqq[i][j] = ((!i) ? (pwqq[i][j - 1] * q2n) : (pwqq[i - 1][j] * q1n));
}
}
for (int i = 1; i <= K; i++) {
Zi tmp = work(n, i) * s1[K][i] * _fac[K];
if ((K - i) & 1) {
ans -= tmp;
} else {
ans += tmp;
}
}
return ans;
}
void __main__(ll l, ll r, int K) {
++l, ++r;
Zi ans = solve(r, K) - solve(l - 1, K);
ans = ans * ~Zi(r - l + 1);
printf("%d\n", ans.v);
}
}
namespace subtask2 {
typedef ComplexTemp<3> Complex;
const Zi inv2 ((Mod + 1) >> 1);
const Zi inv3 ((Mod + 1) / 3);
const Zi inv6 = inv2 * inv3;
const Complex c1 (inv2, -inv6), c2 (inv2, inv6);
const Complex q1 (2, Mod - 1), q2 (2, 1);
Complex pwq1[Kmx], pwq2[Kmx], pwc1[Kmx], pwc2[Kmx];
Complex pwqq[Kmx][Kmx];
Complex get_sum(int k1, int k2, ll n) {
Complex x = pwq1[k1] * pwq2[k2];
if (x == Zi(1))
return Zi(n);
return (pwqq[k1][k2] - 1) / (x - 1);
}
inline void init() {
pwq1[0] = pwq2[0] = pwc1[0] = pwc2[0] = Zi(1);
for (int k = 1; k < Kmx; k++) {
pwq1[k] = pwq1[k - 1] * q1;
pwq2[k] = pwq2[k - 1] * q2;
pwc1[k] = pwc1[k - 1] * c1;
pwc2[k] = pwc2[k - 1] * c2;
}
}
Zi work(ll n, int K) {
Complex ret (0, 0);
for (int k = 0; k <= K; k++) {
Complex tmp = get_sum(k, K - k, n) * comb[K][k];
Complex b = pwc1[k] * pwc2[K - k];
tmp = tmp * b;
ret = ret + tmp;
}
assert(ret.v.v == 0);
// cerr << n << " " << K << " " << ret.r.v << ‘\n‘;
return ret.r;
}
Zi solve(ll n, int K) {
n >>= 1;
Zi ans = 0;
Complex q1n = qpow(q1, n + 1);
Complex q2n = qpow(q2, n + 1);
pwqq[0][0] = Zi(1);
for (int i = 0; i <= K; i++) {
for (int j = (i == 0); i + j <= K; j++) {
pwqq[i][j] = ((!i) ? (pwqq[i][j - 1] * q2n) : (pwqq[i - 1][j] * q1n));
}
}
for (int i = 1; i <= K; i++) {
Zi tmp = work(n, i) * s1[K][i] * _fac[K];
if ((K - i) & 1) {
ans -= tmp;
} else {
ans += tmp;
}
}
return ans;
}
void __main__(ll l, ll r, int K) {
Zi ans = solve(r, K) - solve(l - 1, K);
ans = ans * ~Zi(r - l + 1);
printf("%d\n", ans.v);
}
}
int Case, ___;
ll l, r;
int K;
int main() {
prepare(502);
scanf("%d%d", &Case, &___);
if (___ == 3) {
subtask2::init();
} else {
subtask1::init();
}
while (Case--) {
scanf("%lld%lld%d", &l, &r, &K);
if (___ == 2) {
subtask1::__main__(l, r, K);
} else {
subtask2::__main__(l, r, K);
}
}
return 0;
}loj 3090 「BJOI2019」勘破神機 - 數學