hdu4914 Linear recursive sequence

用矩阵求解线性递推式通项

用fft优化矩阵乘法

首先把递推式求解转化为矩阵求幂，再利用特征多项式f(λ)满足f(A) = 0,将矩阵求幂转化为多项式相乘，

最后利用傅里叶变换的高效算法（迭代取代递归）（参见算法导论）解决。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <map>
#include <string>
#include <vector>
#include <set>
#include <cmath>
#include <ctime>
#pragma comment(linker, "/STACK:102400000,102400000")
using namespace std;
#define lson (u << 1)
#define rson (u << 1 | 1)
typedef long long ll;
const double eps = 1e-6;
const double pi = acos(-1.0);
const int maxn = 4e4 + 10;
const int maxm = 1050;
const int mod = 119;
const int inf = 0x3f3f3f3f;

int n, a, b, p, q;
int size;
int f[maxn], g[maxn];

struct Complex{
    double ii, ij;//ii::real, ij::image
    Complex(double ii = 0, double ij = 0) : ii(ii), ij(ij) {}
//    Complex clear() { this->ii = this->ij = 0; }
    Complex operator + (const Complex &rhs) const{
        return Complex(ii + rhs.ii, ij + rhs.ij);
    }
    Complex operator - (const Complex &rhs) const{
        return Complex(ii - rhs.ii, ij - rhs.ij);
    }
    Complex operator * (const Complex &rhs) const{
        return Complex(ii * rhs.ii - ij * rhs.ij, ii * rhs.ij + ij * rhs.ii);
    }
};

Complex a1[maxn], a2[maxn];

void fft(Complex *src, int len, int rev){
    //len is power of 2
    //rev == 1::dft rev == -1::idft
    for(int i = 1, j = 0; i < len; i++){
        for(int k = len >> 1; k > (j ^= k); k >>= 1) ;
        if(i < j) swap(src[i], src[j]);
    }
    for(int i = 2; i <= len; i <<= 1){
        Complex wi(cos(2 * pi * rev / i), sin(2 * pi * rev / i));
        //(wi)^i = 1
        for(int j = 0; j < len; j += i){
            //using iteration insetad of recursion
            Complex w(1.0, 0.0);
            //w = (wi)^0
            for(int k = j; k < j + i / 2; k++){
                Complex tem = w * src[k + i / 2];
                src[k + i / 2] = src[k] - tem;
                src[k] = src[k] + tem;
                w = w * wi;
            }
        }
    }
    if(rev == -1){
        for(int i = 0; i < len; i++) src[i].ii = (src[i].ii / len + eps);
    }
}

void multi(int *src1, int *src2, int len){
    for(int i = 0; i < len; i++){
        a1[i].ii = a1[i].ij = a2[i].ii = a2[i].ij = 0;
        if(i < q){
            a1[i].ii = (double)src1[i];
            a2[i].ii = (double)src2[i];
        }
    }
    fft(a1, len, 1), fft(a2, len, 1);
    for(int i = 0; i < len; i++) a1[i] = a1[i] * a2[i];
    fft(a1, len, -1);
    for(int i = 0; i < len; i++) g[i] = (int)((ll)(a1[i].ii + eps) % mod);
    for(int i = 2 * q - 2; i >= q; i--){
        //this is because for the fisrt row in matrix A,
        //which satisfies ths (f(n + q),...,f(n))T = A((f(n + q - 1),...,f(n - 1))T)
        //only two elements are nonzero integers
        g[i - q] = (g[i - q] + g[i] * b) % mod;
        g[i - p] = (g[i - p] + g[i] * a) % mod;
    }
    memcpy(src1, g, sizeof(int) * q);
}

int tmp[maxn], ans[maxn];

int main(){
    //freopen("in.txt", "r", stdin);
    while(~scanf("%d%d%d%d%d", &n, &a, &b, &p, &q)){
        a %= mod, b %= mod;
        f[0] = 1;
        for(int i = 1; i < q; i++){
            f[i] = i < p ? a + b : a * f[i - p] + b;
            f[i] %= mod;
        }
        if(n < q){
            printf("%d
", f[n]);
            continue;
        }
        size = 1;
        while(size <= (q - 1) * 2) size <<= 1;
        memset(tmp, 0, sizeof tmp);
        memset(ans, 0, sizeof ans);
        ans[0] = tmp[1] = 1;
        while(n){
            if(n & 1) multi(ans, tmp, size);
            multi(tmp, tmp, size);
            n >>= 1;
        }
        int res = 0;
        for(int i = 0; i < q; i++){
            res = (res + ans[i] * f[i]) % mod;
        }
        printf("%d
", res);
    }
    return 0;
}