这道题目又跟fib有关,不过比一般的fib矩阵难度提高了,关键要有比较强的数学思想,
还要掌握牢固的矩阵公式。
/* *State: 1588 0MS 312K 4305 B *题目大意: * 有等差数列:g(i)=k*i+b; * fib数列: * f(0)=0 * f(1)=1 * f(n)=f(n-1)+f(n-2) (n>=2) * 有 k,b,n ,calculate the sum of every f(g(i)) for 0<=i<n *解题思路: * 如下。其实挺简单。 *解题感想: * 题目还有困惑,就是最后怎么是取ans(0, 1)? */ /* 把斐波那契数列转化为矩阵: A={1 1} {1,0}; {f[n+1],f[n]} {f[n],f[n-1]} = A^n ;最后输出右上角那项 或者用 {f[n+2],f[n+1]} {f[n+1], f[n] } = A^(n+1); 最后输出右下角那项 我们用第一个公式 所求即为A^b + A^(k+b) + A^(2*k+b) + ... + A^((n-1)*k+b) =A^b * ( A^0 + A^k + A^(2*k) + ... + A^((n-1)*k) ) =A^b * ( (A^k)^0 + (A^k)^1 + (A^k)^2 + ...+ (A^k)^(n-1) ); B=A^k; 上式=A^b * ( B^0 + B^1 + B^2 + ... + B^(n-1) ); B^0 + B^1 + B^2 + ... + B^(n-1)用上篇介绍到的递归二分 方法求解 最后输出矩阵的第二项(右上角)即可; 对于求解 B^0 + B^1 + B^2 + ... + B^(n-1) 我们也可以构造矩阵的方法。 我们来设置这样一个矩阵 B I O I 其中O是零矩阵,I是单位矩阵 将它乘方,得到 B^2 I+B O I 乘三方,得到 B^3 I+B+B^2 O I 乘四方,得到 B^4 I+B+B^2+B^3 O I 用快速幂求出n方,让第二项再与A^b相乘即可。 */
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include <stdio.h> 2 #include <string.h> 3 #include <iostream> 4 #define MAX_DIMENSION 5 5 6 using namespace std; 7 8 typedef int MATRIX_TYPE; 9 typedef __int64 MAX_INT_TYPE; //temporary var to avoid overflowing 10 typedef MATRIX_TYPE Matrix[MAX_DIMENSION][MAX_DIMENSION]; 11 int ndim=MAX_DIMENSION; 12 int mod=1;//mod 13 14 void m_zero(Matrix x) 15 { 16 memset(x, 0, sizeof(MATRIX_TYPE)*MAX_DIMENSION*MAX_DIMENSION); 17 } 18 19 void m_one(Matrix x) 20 { 21 int i; 22 m_zero(x); 23 for(i=0;i<ndim;++i)x[i][i]=1; 24 } 25 26 void m_copy(Matrix src,Matrix dest) 27 { 28 memcpy(dest,src, sizeof(MATRIX_TYPE)*MAX_DIMENSION*MAX_DIMENSION); 29 } 30 31 //z=x+y; 32 void m_add(Matrix x,Matrix y,Matrix z) 33 { 34 int i,j; 35 for(i=0;i<ndim;i++) 36 for(j=0;j<ndim;j++) 37 if(mod<=1)z[i][j]=x[i][j]+y[i][j]; 38 else z[i][j]=(x[i][j]+(MAX_INT_TYPE)y[i][j])%mod;//module 39 } 40 41 42 //c=a*b 43 void m_multiple(Matrix a,Matrix b,Matrix c) 44 { 45 int i,j,k; 46 MAX_INT_TYPE t; 47 48 for(i=0;i<ndim;i++) 49 for(j=0;j<ndim;j++) 50 { 51 t=0; 52 if(mod<=1) 53 for(k=0;k<ndim;k++) t+=a[i][k]*b[k][j];//module 54 else 55 for(k=0;k<ndim;k++){ 56 t+=(a[i][k]*(MAX_INT_TYPE)b[k][j])%mod; 57 t%=mod; 58 }//module 59 c[i][j]=t; 60 } 61 } 62 63 void m_pow(Matrix x, unsigned int n, Matrix y) 64 { 65 Matrix temp,temp_x; 66 m_one(y); 67 m_copy(x,temp_x); 68 for(;n;) 69 { 70 if ((n & 1) != 0) 71 { 72 m_multiple(temp_x,y,temp); 73 m_copy(temp,y); 74 } 75 if ((n >>= 1)) 76 { 77 m_multiple(temp_x,temp_x,temp); 78 m_copy(temp,temp_x); 79 } 80 } 81 } 82 83 84 void m_pow_with_saved(Matrix x_p[],unsigned int n, Matrix y) 85 { 86 Matrix temp; 87 m_one(y); 88 for(int k=1;n;++k,n>>=1) 89 { 90 if ((n & 1) != 0) 91 { 92 m_multiple(x_p[k],y,temp); 93 m_copy(temp,y); 94 } 95 } 96 } 97 98 99 void m_pow_sum1(Matrix x_p[],unsigned int n, Matrix y) 100 { 101 m_zero(y); 102 if(n==0)return; 103 if(n==1) m_copy(x_p[1],y); 104 else 105 { 106 Matrix temp; 107 //calculate x^1+...+^(n/2) 108 m_pow_sum1(x_p,n>>1,temp); 109 //add to y 110 m_add(temp,y,y); 111 //calculate x^(1+n/2)+...+x^n 112 Matrix temp2; 113 m_pow_with_saved(x_p,n>>1,temp2); 114 Matrix temp3; 115 m_multiple(temp,temp2,temp3); 116 //add to y 117 m_add(temp3,y,y); 118 if(n&1) 119 { 120 //calculate x^(n-1) 121 m_multiple(temp2,temp2,temp3); 122 //calculate x^n 123 m_multiple(temp3,x_p[1],temp2); 124 //add x^n 125 m_add(temp2,y,y); 126 } 127 } 128 129 } 130 131 void m_pow_sum(Matrix x, unsigned int n, Matrix y) 132 { 133 unsigned int i=0,logn,n2=n; 134 for(logn=0,n2=n;n2;logn++,n2 >>= 1); 135 Matrix *pow_arry=new Matrix[logn==0?2:(logn+1)]; 136 m_one(pow_arry[0]); 137 m_copy(x,pow_arry[1]); 138 for(i=1;i<logn;i++) 139 { 140 m_multiple(pow_arry[i],pow_arry[i],pow_arry[i+1]); 141 } 142 143 m_pow_sum1(pow_arry,n,y); 144 delete []pow_arry; 145 } 146 147 void fib(int m) 148 { 149 Matrix a={{1,1},{1,0}},y; 150 ndim=2; 151 mod=1024*1024; 152 153 m_pow(a,m,y); 154 //fib(m+1) 155 printf("fib(%d)=%d\n",m,y[0][0]); 156 157 m_pow_sum(a,m,y); 158 //fib(m+1) 159 printf("sum(fib(i))=%d,i=1,..,%d\n",y[0][0],m); 160 161 } 162 163 void view_arr(Matrix a, int n) 164 { 165 for(int i = 0; i < n; i++) 166 { 167 for(int j = 0; j < n;j++) 168 cout << a[i][j] << " "; 169 cout << endl; 170 } 171 } 172 173 int main() 174 { 175 #ifndef ONLINE_JUDGE 176 freopen("in.txt", "r", stdin); 177 #endif 178 Matrix a = {{1, 1}, {1, 0}}; 179 //view_arr(a, 2); 180 ndim = 2; 181 int k, b, n, M; 182 while(scanf("%d %d %d %d", &k, &b, &n, &M) == 4) 183 { 184 mod = M; 185 186 Matrix a_k, a_b, a_k_sum, res; 187 m_one(res); 188 m_pow(a, k, a_k); 189 m_pow(a, b, a_b); 190 m_pow_sum(a_k, n - 1, a_k_sum); 191 Matrix addans, ans; 192 m_add(res, a_k_sum, addans); 193 m_multiple(addans, a_b, ans); 194 printf("%d\n", ans[0][1]); 195 //view_arr(ans, 2); 196 //cout << endl; 197 } 198 return 0; 199 }