【洛谷】P1962
题目链接：https://www.luogu.org/problemnew/show/P1962
题意：求fib数列的第n项，很大。mod 1e9+7.
题解：BM直接推。
代码：
1 #include <cstdio> 2 #include <cstring> 3 #include <cmath> 4 #include <algorithm> 5 #include <vector> 6 #include <string> 7 #include <map> 8 #include <set> 9 #include <cassert> 10 using namespace std; 11 #define rep(i,a,n) for (int i=a;i<n;i++) 12 #define per(i,a,n) for (int i=n-1;i>=a;i--) 13 #define pb push_back 14 #define mp make_pair 15 #define all(x) (x).begin(),(x).end() 16 #define fi first 17 #define se second 18 #define SZ(x) ((int)(x).size()) 19 typedef vector<int> VI; 20 typedef long long ll; 21 typedef pair<int, int> PII; 22 const ll mod = 1000000007; 23 ll powmod(ll a, ll b) 24 { 25 ll res = 1; a %= mod; 26 assert(b >= 0); 27 for (; b; b >>= 1) 28 { 29 if (b & 1) 30 res = res * a%mod; 31 a = a * a%mod; 32 } 33 return res; 34 } 35 // head 36 37 int _, n; 38 namespace linear_seq 39 { 40 const int N = 10010; 41 ll res[N], base[N], _c[N], _md[N]; 42 vector<int> Md; 43 void mul(ll *a, ll *b, int k) 44 { 45 rep(i, 0, k + k) _c[i] = 0; 46 rep(i, 0, k) 47 if (a[i]) 48 rep(j, 0, k) 49 _c[i + j] = (_c[i + j] + a[i] * b[j]) % mod; 50 for (int i = k + k - 1; i >= k; i--) 51 if (_c[i]) 52 rep(j, 0, SZ(Md)) 53 _c[i - k + Md[j]] = (_c[i - k + Md[j]] - _c[i] * _md[Md[j]]) % mod; 54 rep(i, 0, k) a[i] = _c[i]; 55 } 56 int solve(ll n, VI a, VI b) { // a 系数 b 初值 b[n+1]=a[0]*b[n]+... 57 // printf("%d\n",SZ(b)); 58 ll ans = 0, pnt = 0; 59 int k = SZ(a); 60 assert(SZ(a) == SZ(b)); 61 rep(i, 0, k) 62 _md[k - 1 - i] = -a[i]; _md[k] = 1; 63 Md.clear(); 64 rep(i, 0, k) 65 if (_md[i] != 0) Md.push_back(i); 66 rep(i, 0, k) 67 res[i] = base[i] = 0; 68 res[0] = 1; 69 while ((1ll << pnt) <= n) pnt++; 70 for (int p = pnt; p >= 0; p--) 71 { 72 mul(res, res, k); 73 if ((n >> p) & 1) 74 { 75 for (int i = k - 1; i >= 0; i--) res[i + 1] = res[i]; res[0] = 0; 76 rep(j, 0, SZ(Md)) res[Md[j]] = (res[Md[j]] - res[k] * _md[Md[j]]) % mod; 77 } 78 } 79 rep(i, 0, k) ans = (ans + res[i] * b[i]) % mod; 80 if (ans < 0) ans += mod; 81 return ans; 82 } 83 VI BM(VI s) 84 { 85 VI C(1, 1), B(1, 1); 86 int L = 0, m = 1, b = 1; 87 rep(n, 0, SZ(s)) 88 { 89 ll d = 0; 90 rep(i, 0, L + 1) d = (d + (ll)C[i] * s[n - i]) % mod; 91 if (d == 0) ++m; 92 else if (2 * L <= n) 93 { 94 VI T = C; 95 ll c = mod - d * powmod(b, mod - 2) % mod; 96 while (SZ(C) < SZ(B) + m) C.pb(0); 97 rep(i, 0, SZ(B)) C[i + m] = (C[i + m] + c * B[i]) % mod; 98 L = n + 1 - L; B = T; b = d; m = 1; 99 } 100 else 101 { 102 ll c = mod - d * powmod(b, mod - 2) % mod; 103 while (SZ(C) < SZ(B) + m) C.pb(0); 104 rep(i, 0, SZ(B)) C[i + m] = (C[i + m] + c * B[i]) % mod; 105 ++m; 106 } 107 } 108 return C; 109 } 110 int gao(VI a, ll n) 111 { 112 VI c = BM(a); 113 c.erase(c.begin()); 114 rep(i, 0, SZ(c)) c[i] = (mod - c[i]) % mod; 115 return solve(n, c, VI(a.begin(), a.begin() + SZ(c))); 116 } 117 }; 118 119 int main() { 120 long long n; 121 vector<int>v; 122 v.push_back(1); 123 v.push_back(1); 124 v.push_back(2); 125 v.push_back(3); 126 v.push_back(5); 127 v.push_back(8); 128 v.push_back(13); 129 v.push_back(21); 130 v.push_back(34); 131 while(~scanf("%lld",&n)){ 132 printf("%lld\n", linear_seq::gao(v, n - 1)); 133 //VI{1,2,4,7,13,24} 134 //printf("%lld\n", linear_seq::gao(v, n - 1)); 135 //printf("%d\n",linear_seq::gao(VI{x1,x2,x3,x4},n-1)); 136 } 137 }
update:
矩阵快速幂
\begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix} *
\begin{pmatrix} f(n-1)\\ f(n-2) \end{pmatrix} =
\begin{pmatrix} f(n)\\ f(n-1) \end{pmatrix}
代码：
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 using namespace std; 5 #define ll long long 6 const int maxn = 3; 7 const ll mod = 1e9+7; 8 9 ll n,p; 10 11 //矩阵结构体 12 struct Matrix{ 13 ll a[maxn][maxn]; 14 void init(){ //初始化为单位矩阵 15 memset(a, 0, sizeof(a)); 16 for(int i = 1; i <= maxn;i++){ 17 a[i][i] = 1; 18 } 19 } 20 }; 21 22 //矩阵乘法 23 Matrix mul(Matrix a, Matrix b){ 24 Matrix ans; 25 for(int i = 1;i <= 2;++i){ 26 for(int j = 1;j <= 2;++j){ 27 ans.a[i][j] = 0; 28 for(int k = 1;k <= 2;++k){ 29 ans.a[i][j] = ans.a[i][j] % mod + a.a[i][k] * b.a[k][j] % mod; 30 } 31 } 32 } 33 return ans; 34 } 35 36 //矩阵快速幂 37 Matrix qpow(Matrix a,ll b){ 38 Matrix ans; 39 ans.init(); 40 while(b){ 41 if(b & 1) 42 ans = mul(ans,a); 43 a = mul(a,a); 44 b >>= 1; 45 } 46 return ans; 47 } 48 49 void print(Matrix a){ 50 for(int i = 1; i <= n;++i){ 51 for(int j = 1;j <= n;++j){ 52 cout << a.a[i][j]%mod<< " "; 53 } 54 cout << endl; 55 } 56 } 57 58 int main(){ 59 Matrix base; 60 Matrix ans; 61 ans.a[1][1] = 1; 62 ans.a[1][2] = 0; 63 base.a[1][1] = 1; 64 base.a[1][2] = 1; 65 base.a[2][1] = 1; 66 base.a[2][2] = 0; 67 cin>>n; 68 ans = mul(ans,qpow(base,n-1)); 69 cout<<ans.a[1][1]<<endl; 70 return 0; 71 }
