标签：out dbi ems tip ota udt udp imu pxe
传送门：https://jzoj.net/senior/#main/show/5142
【题目大意】
给出n个数，a[1]...a[n]，称作集合S，求
其中f[i] = 2f[i-1] + 3f[i-2]，给出f[0],f[1]。mod 99991
n<=100000
【题解】
暴力dp，用矩阵作为存储值，复杂度O(n^2)
# include <ctype.h> # include <stdio.h> # include <assert.h> # include <iostream> # include <string.h> # include <algorithm> using namespace std; typedef long long ll; typedef unsigned long long ull; typedef long double ld; const int N = 1e5 + 10, M = 2e5 + 10; const int mod = 99991; inline int getint() { int x = 0, f = 1; char ch = getchar(); while(!isdigit(ch)) { if(ch == ‘-‘) f = 0; ch = getchar(); } while(isdigit(ch)) { x = (x<<3) + (x<<1) + ch - ‘0‘; ch = getchar(); } return f ? x : -x; } struct mat { int n, m, a[3][3]; inline void set(int _n, int _m) { n = _n, m = _m; memset(a, 0, sizeof a); } friend mat operator * (mat a, mat b) { mat c; c.set(a.n, b.m); assert(a.m == b.n); for (int i=1; i<=c.n; ++i) for (int j=1; j<=c.m; ++j) for (int k=1; k<=a.m; ++k) { c.a[i][j] += 1ll * a.a[i][k] * b.a[k][j] % mod; if(c.a[i][j] >= mod) c.a[i][j] -= mod; } return c; } friend mat operator ^ (mat a, ll b) { mat c; c.set(a.n, a.n); assert(a.n == a.m); for (int i=1; i<=c.n; ++i) c.a[i][i] = 1; while(b) { if(b&1) c = c * a; a = a * a; b >>= 1; } return c; } friend mat operator + (mat a, mat b) { assert(a.n == b.n && a.m == b.m); for (int i=1; i<=a.n; ++i) for (int j=1; j<=a.m; ++j) { a.a[i][j] += b.a[i][j]; if(a.a[i][j] >= mod) a.a[i][j] -= mod; } return a; } }A, T, U, g[2][5010], e[M]; inline int f(ll n) { if(n == 0) return A.a[1][1]; U = (T^(n-1)) * A; return U.a[2][1]; } int n, K, a[M]; inline void force() { int now = 1, pre = 0; for (int i=0; i<=1; ++i) for (int j=0; j<=n; ++j) g[i][j].set(2, 1); for (int i=1; i<=n; ++i) { g[now][0] = g[pre][0] = A; for (int j=1, jto = min(i, K); j<=jto; ++j) { g[now][j] = g[pre][j] + (e[i] * g[pre][j-1]); // printf("%d %d %d\n", i, j, g[i][j].a[1][1]); } swap(pre, now); } printf("%d\n", g[pre][K].a[1][1]); } int main() { // freopen("see.in", "r", stdin); // freopen("see.out", "w", stdout); n = getint(), K = getint(); for (int i=1; i<=n; ++i) a[i] = getint(); A.set(2, 1); T.set(2, 2); cin >> A.a[1][1] >> A.a[2][1]; T.a[2][1] = 3, T.a[2][2] = 2, T.a[1][1] = 0, T.a[1][2] = 1; for (int i=1; i<=n; ++i) e[i] = (T^a[i]); force(); return 0; } /* 4 2 1 2 1 3 1 1 20 10 125 3162 3261 152 376 238 462 4382 376 4972 16 1872 463 9878 688 308 125 236 3526 543 1223 3412 */
根据特征根那套理论，设f[i] = x^i,那么有x^2-2x-3=0，解出x=3和x=-1两个根。
那么f[i] = a*3^i + b*(-1)^i
带入f[0],f[1]可以解出a,b（取模意义下）
然后f[a1+a2+...+ak] = a * 3^a1 * 3^a2 * ... * 3^ak + b * (-1)^a1 * (-1)*a2 * ... * (-1)^ak。（取模意义下）
明显3和-1互不干扰，所以分开来算。
n个数选出k个，提示了我们生成函数。
对于每个a[i]，构造多项式(1+base^a[i] * x)，base为3/-1.
然后答案就是(1+base^a[1]) * (1+base^a[2]) * ... * (1+base^a[n])
这个东西可以分治FFT解决，O(nlogn^2n)
做两遍分治FFT即可解决问题。
# include <math.h> # include <vector> # include <stdio.h> # include <iostream> # include <string.h> # include <algorithm> using namespace std; typedef long long ll; typedef unsigned long long ull; typedef long double ld; # define RG register # define ST static const int N = 1e5 + 10, M = 2e5 + 10; const int mod = 99991; const double pi = acos(-1.0); inline int getint() { int x = 0, f = 1; char ch = getchar(); while(!isdigit(ch)) { if(ch == ‘-‘) f = 0; ch = getchar(); } while(isdigit(ch)) { x = (x<<3) + (x<<1) + ch - ‘0‘; ch = getchar(); } return f ? x : -x; } // x^2 - 2x - 3 = 0 // x1 = -1, x2 = 3 // f(x) = A * (-1)^x + B * 3^x // f(x) = A * (-1)^{a1+a2+...+ak} + B * 3^{a1+a2+...+ak} // f(x) = A * (-1)^{a1} * (-1)^{a2} * ... + B * 3^{a1} * 3^{a2} * ... // (1 + 3^a) ST int n, K, a[M], sr[2][M]; ST int A, B, F0, F1; struct cp { double x, y; cp() {} cp(double x, double y) : x(x), y(y) {} friend cp operator + (cp a, cp b) { return cp(a.x+b.x, a.y+b.y); } friend cp operator - (cp a, cp b) { return cp(a.x-b.x, a.y-b.y); } friend cp operator * (cp a, cp b) { return cp(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x); } }; ST cp u[M], v[M]; inline int pwr(int a, int b) { int ret = 1; while(b) { if(b&1) ret = 1ll * ret * a % mod; a = 1ll * a * a % mod; b >>= 1; } return ret; } ST int p[M]; namespace FFT { ST cp w[2][M]; ST int n, lst[M]; inline void init(int _n) { n = 1; while(n < _n) n <<= 1; for (RG int i=0; i<n; ++i) w[0][i] = cp(cos(pi * 2.0 / n * i), sin(pi * 2.0 / n * i)), w[1][i] = cp(w[0][i].x, -w[0][i].y); RG int len = 0; while((1 << len) < n) ++len; for (RG int i=0; i<n; ++i) { int t = 0; for (RG int j=0; j<len; ++j) if(i & (1 << j)) t |= (1 << (len - j - 1)); lst[i] = t; } } inline void DFT(cp *a, int op) { cp *o = w[op]; for (RG int i=0; i<n; ++i) if(i < lst[i]) swap(a[i], a[lst[i]]); for (RG int len=2; len<=n; len <<= 1) { int m = len>>1; for (RG cp *p = a; p != a+n; p += len) { for (RG int k=0; k<m; ++k) { cp t = o[n/len*k] * p[k+m]; p[k+m] = p[k] - t; p[k] = p[k] + t; } } } if(op == 1) for (RG int i=0; i<n; ++i) a[i].x /= (double)n; } inline void merge(int l, int mid, int r) { int m = r-l+2; init(m); // printf("%d\n", n); u[0] = cp(1.0, 0.0), v[0] = cp(1.0, 0.0); for (RG int i=1; i<=mid-l+1; ++i) u[i] = cp(p[i+l-1], 0.0); for (RG int i=1; i<=r-mid; ++i) v[i] = cp(p[i+mid], 0.0); for (RG int i=mid-l+2; i<n; ++i) u[i] = cp(0.0, 0.0); for (RG int i=r-mid+1; i<n; ++i) v[i] = cp(0.0, 0.0); DFT(u, 0); DFT(v, 0); for (RG int i=0; i<n; ++i) u[i] = u[i] * v[i]; DFT(u, 1); for (RG int i=1; i<m; ++i) p[l+i-1] = ((ll)(u[i].x + 0.5)) % mod; } } int id; inline void solve(int l, int r) { if(l == r) { p[l] = sr[id][l]; return ; } RG int mid = l+r>>1; solve(l, mid); solve(mid+1, r); FFT::merge(l, mid, r); // printf("l = %d, r = %d\n", l, r); // for (int i=0; i<p[x].size(); ++i) printf("%d ", p[x][i]); // puts(""); } int main() { // freopen("see.in", "r", stdin); // freopen("see.out", "w", stdout); n = getint(), K = getint(); for (RG int i=1; i<=n; ++i) { a[i] = getint(); sr[0][i] = pwr(3, a[i]); sr[1][i] = pwr(mod-1, a[i]); } F0 = getint(), F1 = getint(); B = 1ll * (F0 + F1) * pwr(4, mod-2) % mod; A = (F0 - B + mod) % mod; id = 0; solve(1, n); int t = 1ll * B * p[K] % mod; id = 1; solve(1, n); t += 1ll * A * p[K] % mod; if(t >= mod) t -= mod; printf("%d\n", t); return 0; }
顺便贴一个被卡常的分治FFT.......
# include <math.h> # include <vector> # include <stdio.h> # include <iostream> # include <string.h> # include <algorithm> using namespace std; typedef long long ll; typedef unsigned long long ull; typedef long double ld; # define RG register # define ST static const int N = 1e5 + 10, M = 2e5 + 10; const int mod = 99991; const double pi = acos(-1.0); inline int getint() { int x = 0, f = 1; char ch = getchar(); while(!isdigit(ch)) { if(ch == ‘-‘) f = 0; ch = getchar(); } while(isdigit(ch)) { x = (x<<3) + (x<<1) + ch - ‘0‘; ch = getchar(); } return f ? x : -x; } // x^2 - 2x - 3 = 0 // x1 = -1, x2 = 3 // f(x) = A * (-1)^x + B * 3^x // f(x) = A * (-1)^{a1+a2+...+ak} + B * 3^{a1+a2+...+ak} // f(x) = A * (-1)^{a1} * (-1)^{a2} * ... + B * 3^{a1} * 3^{a2} * ... // (1 + 3^a) ST int n, K, a[M], sr[2][M]; ST int A, B, F0, F1; struct cp { double x, y; cp() {} cp(double x, double y) : x(x), y(y) {} friend cp operator + (cp a, cp b) { return cp(a.x+b.x, a.y+b.y); } friend cp operator - (cp a, cp b) { return cp(a.x-b.x, a.y-b.y); } friend cp operator * (cp a, cp b) { return cp(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x); } }; ST cp u[M], v[M]; inline int pwr(int a, int b) { int ret = 1; while(b) { if(b&1) ret = 1ll * ret * a % mod; a = 1ll * a * a % mod; b >>= 1; } return ret; } ST vector<int> x[M]; # define ls (x<<1) # define rs (x<<1|1) namespace FFT { ST cp w[2][M]; ST int n, lst[M]; inline void init(int _n) { n = 1; while(n < _n) n <<= 1; for (RG int i=0; i<n; ++i) w[0][i] = cp(cos(pi * 2.0 / n * i), sin(pi * 2.0 / n * i)), w[1][i] = cp(w[0][i].x, -w[0][i].y); RG int len = 0; while((1 << len) < n) ++len; for (RG int i=0; i<n; ++i) { int t = 0; for (RG int j=0; j<len; ++j) if(i & (1 << j)) t |= (1 << (len - j - 1)); lst[i] = t; } } inline void DFT(cp *a, int op) { cp *o = w[op]; for (RG int i=0; i<n; ++i) if(i < lst[i]) swap(a[i], a[lst[i]]); for (RG int len=2; len<=n; len <<= 1) { int m = len>>1; for (RG cp *p = a; p != a+n; p += len) { for (RG int k=0; k<m; ++k) { cp t = o[n/len*k] * p[k+m]; p[k+m] = p[k] - t; p[k] = p[k] + t; } } } if(op == 1) for (RG int i=0; i<n; ++i) a[i].x /= (double)n; } inline vector<int> merge(vector<int> a, vector<int> b) { vector<int> ret; ret.clear(); int m = max(a.size(), b.size()) * 2; init(m); // printf("%d\n", n); for (RG int i=0; i<a.size(); ++i) u[i] = cp(a[i], 0.0); for (RG int i=0; i<b.size(); ++i) v[i] = cp(b[i], 0.0); for (RG int i=a.size(); i<n; ++i) u[i] = cp(0.0, 0.0); for (RG int i=b.size(); i<n; ++i) v[i] = cp(0.0, 0.0); DFT(u, 0); DFT(v, 0); for (RG int i=0; i<n; ++i) u[i] = u[i] * v[i]; DFT(u, 1); for (RG int i=0; i<K+1 && i<m; ++i) ret.push_back(((ll)(u[i].x + 0.5)) % mod); return ret; } } int id; inline void solve(int x, int l, int r) { if(l == r) { p[x].clear(); p[x].push_back(1), p[x].push_back(sr[id][l]); return ; } RG int mid = l+r>>1; solve(ls, l, mid); solve(rs, mid+1, r); p[x] = FFT::merge(p[ls], p[rs]); // printf("l = %d, r = %d\n", l, r); // for (int i=0; i<p[x].size(); ++i) printf("%d ", p[x][i]); // puts(""); } # include <time.h> int main() { // freopen("see.in", "r", stdin); // freopen("see.out", "w", stdout); int beg = clock(); n = getint(), K = getint(); for (RG int i=1; i<=n; ++i) { a[i] = getint(); sr[0][i] = pwr(3, a[i]); sr[1][i] = pwr(mod-1, a[i]); } F0 = getint(), F1 = getint(); B = 1ll * (F0 + F1) * pwr(4, mod-2) % mod; A = (F0 - B + mod) % mod; id = 0; solve(1, 1, n); int t = 1ll * B * p[1][K] % mod; id = 1; solve(1, 1, n); t += 1ll * A * p[1][K] % mod; if(t >= mod) t -= mod; printf("%d\n", t); int end = clock() - beg; cerr << end << " ms" << endl; return 0; }
标签：out dbi ems tip ota udt udp imu pxe
原文地址：http://www.cnblogs.com/galaxies/p/jzoj5142.html
