According to test instructions: The last is to ask F (b) + f (k + B) + F (2 * k + b) + ... + f ((n-1) * k + b)
Obviously f (b) = A^b
Where a =
1 1
1 0
So sum (n-1) = a^b (E + a^ k + A ^ (2 * k) + ... + a ^ ((n-1) * k)
Set d = a^k
SUM (n-1) = A^b (E + D + D ^ 2 + ... + D ^ (n-1))
Part of the brackets can be divided into two-point recursion.
and a single matrix can be obtained with a matrix fast power.
/************************************************************************* > File Name:hdu1588.cpp > Auth Or:alex > Mail: [email protected] > Created time:2015 March 12 Thursday 18:25 07 seconds ******************************** ****************************************/#include <map>#include <set>#include <queue>#include <stack>#include <vector>#include <cmath>#include <cstdio>#include <cstdlib>#include <cstring>#include <iostream>#include <algorithm>using namespace STD;Const DoublePI =ACOs(-1.0);Const intINF =0x3f3f3f3f;Const DoubleEPS =1e-15;typedef Long LongLL;typedefPair <int,int> PLL; LL MoD, k, b;classmartix{ Public: LL mat[3][3]; Martix (); Martixoperator* (ConstMartix &b)Const; Martixoperator+ (ConstMartix &b)Const; martix&operator= (ConstMartix &b);} A, E, D; Martix:: Martix () {memset(Mat,0,sizeof(MAT));} Martix Martix::operator* (ConstMartix &b)Const{Martix ret; for(inti =0; I <2; ++i) { for(intj =0; J <2; ++J) { for(intK =0; K <2; ++K) {Ret.mat[i][j] + = ThisMAT[I][K] * B.mat[k][j]; RET.MAT[I][J]%= mod; } } }returnRET;} Martix Martix::operator+ (ConstMartix &b)Const{Martix ret; for(inti =0; I <2; ++i) { for(intj =0; J <2; ++J) {Ret.mat[i][j] = ThisMAT[I][J] + b.mat[i][j]; RET.MAT[I][J]%= mod; } }returnRET;} martix& Martix::operator= (ConstMartix &b) { for(inti =0; I <2; ++i) { for(intj =0; J <2; ++J) { ThisMAT[I][J] = B.mat[i][j]; } }return* This;} Martix Fastpow (Martix ret, LL N) {Martix ans; ans.mat[0][0] = ans.mat[1][1] =1; while(n) {if(N &1) {ans = ans * RET; } N >>=1; RET = ret * RET; }returnAns;}voidDebug (Martix A) { for(inti =0; I <2; ++i) { for(intj =0; J <2; ++J) {printf("%lld", A.mat[i][j]); }printf("\ n"); }}martix Binseach (LL N) {if(n = =1) {returnD } Martix NXT = Binseach (n >>1); Martix B = Fastpow (D, N/2); B = B + E; NXT = NXT * B;if(N &1) {Martix C = Fastpow (D, N); NXT = NXT + C; }returnNXT;}intMain () {LL n; e.mat[0][0] = e.mat[1][1] =1; a.mat[0][0] = a.mat[0][1] = a.mat[1][0] =1;//Debug (A); while(~scanf("%lld%lld%lld%lld", &k, &b, &n, &mod)) {if(n = =1) {Martix x = Fastpow (A, b);printf("%lld\n", x.mat[0][1]);Continue; } D = Fastpow (A, k); Martix ans = Binseach (n-1); Ans = ans + E; Martix base = Fastpow (A, b); Ans = base * ans;//Debug (ans); printf("%lld\n", ans.mat[0][1]); }return 0;}
hdu1588---Gauss Fibonacci (Matrix, linear recursion)