HDU 4914 Linear recursive sequence(矩阵乘法递推的优化)

最新推荐文章于 2019-02-26 23:49:50 发布

原创最新推荐文章于 2019-02-26 23:49:50 发布 · 2.4k 阅读

0 ·

CC 4.0 BY-SA版权

ACM_数论and数学类专栏收录该内容

55 篇文章

订阅专栏

本文深入探讨了X姐论文中关于矩阵乘法递推优化的方法，并提供了实际的C++代码实现，旨在提高计算效率。主要内容包括复数运算、快速傅里叶变换（FFT）算法以及模运算优化等关键技术。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题解见X姐的论文矩阵乘法递推的优化，只是mark一下。。

#include 
#include 
#include 
#include 
using namespace std;
typedef __int64 ll;
const double pi = acos(-1.0);
const int N = 10000+5;
const double eps = 1e-6;
const int MOD = 119;

struct Complex {
        double x,y;
        Complex(double _x=0,double _y=0)
                        :x(_x),y(_y) {}
        Complex operator + (const Complex& tt)const { return Complex(x+tt.x,y+tt.y); }
        Complex operator - (const Complex& tt)const { return Complex(x-tt.x,y-tt.y); }
        Complex operator * (const Complex& tt)const { return Complex(x*tt.x-y*tt.y,x*tt.y+y*tt.x); }
};
void FFT(Complex *a, int n, int rev) {
        // n是大于等于相乘的两个数组长度的2的幂次
        // 从0开始表示长度，对a进行操作
        // rev==1进行DFT，==-1进行IDFT
        for (int i = 1,j = 0; i < n; ++ i) {
                for (int k = n>>1; k > (j^=k); k >>= 1);
                if (i= MOD)    x -= MOD;
}

void mul(int *tmp1, int *tmp2, int maxn) {
    //for(int i = 0;i < q; i++)   printf("%d ", tmp1[i]); puts("");
    //for(int i = 0;i < q; i++)   printf("%d ", tmp2[i]); puts("");
    for(int i = 0;i < maxn; i++) {
        a1[i].x = a1[i].y = a2[i].x = a2[i].y = 0;
        if(i < q)   {
            a1[i].x = 1.0*tmp1[i];
            a2[i].x = 1.0*tmp2[i];
        }
    }
    FFT(a1, maxn, 1); FFT(a2, maxn, 1);
    for(int i = 0;i < maxn; i++)    a1[i] = a1[i]*a2[i];
    FFT(a1, maxn, -1);
    for(int i = 0;i < maxn; i++)    Fuck[i] = (ll)(a1[i].x+eps)%MOD;
    for(int i = q*2-2;i >= q; i--) {
        if(Fuck[i]) {
            Fuck[i-q] = (Fuck[i-q] + Fuck[i]*b)%MOD;
            Fuck[i-p] = (Fuck[i-p] + Fuck[i]*a)%MOD;
        }
    }
    //printf("q = %d\n", q);
    for(int i = 0;i < q; i++) {
        tmp1[i] = (int)Fuck[i];
    }
}

int ans[N], tmp[N], f[N];

int solve() {
    a %= MOD; b %= MOD;
    f[0] = 1;
    for(int i = 1;i < q; i++) {
        if(i-p < 0) f[i] = a+b;
        else if(i < q)  f[i] = (a*f[i-p]+b)%MOD;
    }
    if(n < q) {
        return f[n];
    }
    int maxn = 1;
    while(maxn <= (q-1)*2)  maxn <<= 1;
    for(int i = 0;i < q; i++)    ans[i] = tmp[i] = 0;
    ans[0] = 1; tmp[1] = 1;
    int fuck = n;
    while(n) {
        if(n&1) {
            mul(ans, tmp, maxn);
        }
        mul(tmp, tmp, maxn);
        n >>= 1;
    }
    int ret = 0;
    for(int i = 0;i < q; i++) {
        Add(ret, ans[i]*f[i]%MOD);
    }
    n = fuck;
    return ret;
}

int main() {
    while(input()) {
        printf("%d\n", solve());
    }
    return 0;
}