PE 440 【二阶递推】

被拖进的坑,555
题目大意：有10种长度为1的木块，1种长度为2的木块，设T(n)表示长度为n的木板拼凑方案数，求∑gcd(T(ca),T(cb))mod987898789(1<=a,b,c<=2000)

容易得到递推式 T(n)=10T(n−1)+T(n−2),首项 T(1)=10,T(2)=101
和Fibonacci数列相似，我们可以推出它一些相关的性质

T(n)的通项公式：T(n)=(5+26√)n+1−(5−26√)n+1226√
可得 T(a+b)=T(a)T(b)+T(a−1)T(b−1)

1) gcd(Tn,Tn−1)==1
证明 :　gcd(Tn,Tn−1)=gcd(10Tn−1+Tn−2,Tn−1)=gcd(Tn−2,Tn−1)⋯=gcd(T2,T1)=1

2) gcd(Ta,Tb)==Tgcd(a+1,b+1)−1
证明：gcd(Ta,Tb)=gcd(Ta−bTb+Ta−b−1Tb−1,Tb)=gcd(Ta−b−1Tb−1,Tb)=gcd(T(a+1)−(b+1)−1,T(b+1)−1)=Tgcd(a+1,b+1)−1

所以可得gcd(T(ca),T(cb))=T(gcd(ca+1,cb+1)−1)
设 g=gcd(a,b)
gcd(ca+1,cb+1)=⎧⎩⎨⎪⎪⎪⎪cg+1c(mod2)+1agandbgareoddotherwise

证明：
gcd(ca+1,cb+1)=gcd(ca−cb,cb+1)=gcd(ca−b−1,cb+1)=gcd(ca−b+cb,cb+1)=gcd(c|a−2b|+1,cb+1)

i) 最终 a = b = g,值为 cg+1
    若ag或bg为偶数那么一定可以再分
ii) 出现 a 或 b == 0，值为 c mod 2 + 1

现在就可以考虑枚举 a,b，这样需要先预处理出sum[c]=∑gcd(T(ca),T(cb))

但是 ca 很大，所以还要找到Tmod987898789的循环节
循环节长度为987898788

//PS:这里并没有形成 ρ 环，蒟蒻表示这里还搞得不是很清楚。。。

所以求T(ca)的时候还要加上矩乘。。。

复杂度O(N2logN),(其实还有个23,就算做常数啦。。。

【答案】970746056

#include<iostream>#include<algorithm>#include<cstdlib>#include<cstdio>#include<string>#include<cstring>#include<ctime>#include<cmath>#define N 2000#define M 987898789#define LL long longusing namespace std;LL ans;LL sum[N+5];LL A[2][2],B[2][2],t[2][2];LL ksm(LL a,LL b,LL c){    LL ret = 1;    for (;b;b >>= 1,a = a * a % c)        if (b & 1) ret = ret * a % c;    return ret;}LL mul(LL a[][2],LL b[][2]){    for (int i = 0;i < 2;i ++)        for (int j = 0;j < 2;j ++)        {            t[i][j] = 0;            for (int k = 0;k < 2;k ++)                (t[i][j] += a[i][k] * b[k][j]) %= M;        }    for (int i = 0;i < 2;i ++)        for (int j = 0;j < 2;j ++)            a[i][j] = t[i][j];}LL matrix_ksm(LL b){    A[0][0] = 10,A[0][1] = A[1][0] = 1,A[1][1] = 0;    B[0][0] = B[1][1] = 1,B[0][1] = B[1][0] = 0;    for (;b;b >>= 1,mul(A,A))        if (b & 1) mul(B,A);    return B[0][0];}int gcd(int a,int b){    return b ? gcd(b,a % b) : a;}int main(){    for (int i = 1;i <= N;i ++)        for (int j = 1;j <= N;j ++)            (sum[i] += matrix_ksm(ksm(j,i,M - 1))) %= M;    for (int a = 1;a <= N;a ++)        for (int b = 1;b <= N;b ++)        {            int g = gcd(a,b);            if (((a / g) & 1) && ((b / g) & 1)) (ans += sum[g]) %= M;                else (ans += 11000) %= M;        }    cout << ans << endl;    return 0;}