洛谷T54037 【最开始】

发表于 | 更新于 | 分类于 | 阅读次数：

题目

给定两个正整数$n,m$，表示问题为$a+\frac{1}{a}=n,ask~~ a^m+\frac{1}{a^m}=?$

思路

  1. 这题一看就是数学题啊！而且好像不难，我想练练！

抱着这样的心情，我运用了$(x^n+y^n)=(x+y)(x^{n-1}-x^{n-2}*y+…)$这个公式，想要递归求出每一项，结果完美的全WA了。。

翻了翻题解，发现标答是用了递推的思想：

$$(a+\frac{1}{a})(a^{n-1}+\frac{1}{a^{n-1}})=(a^n+\frac{1}{a^n})+(a^{n-2}+\frac{1}{a^{n-2}})$$

这样就搞出了$f[x]$与$f[x-1],f[x-2]$的关系，如果直接滚动数组递推的话可得$40~point$。

  1. 如果想进一步优化就得用到计算机知识了。

抱着这样的想法，我把题解往下翻了翻，发现还有一个公式：
$(a^n+\frac{1}{a^n})^2=a^{2n}+\frac{1}{a^{2n}}+2$ ，
$(a^n+\frac{1}{a^n})(a^{n+1}+\frac{1}{a^{n+1}})=a^{2n+1}+\frac{1}{a^{2n+1}}+(a+\frac{1}{a})$

所以本质上应该是一个倍增的过程，代码十行都不到，估计这就是数学的魅力吧。

但是这种公式只是适用于此类式子，并不具有通用性。

既然得到了递推式，我们便可以利用矩阵加速达到同样的目的；

我们设$(a+\frac{1}{a})=k$，便可得到：

$$f[n]+f[n-2]=k*f[n-1]$$

容易构造以下递推式：
$\begin{bmatrix}0&1&0\\0&0&1\\0&{-1}&k\\\end{bmatrix}*\begin{bmatrix}F_{n-2}\\F_{n-1}\\F_n\\\end{bmatrix}=\begin{bmatrix}F_{n-1}\\F_n\\k*F_n-F_{n-1}\end{bmatrix}$

搞个矩阵加速玩玩即可，但是被细节坑的很惨：

  1. 可以注意到构造矩阵中有负数，若取模时候不进行处理可能会爆负，直接在转移时多加一个mo即可。

  2. 这玩意构造基本矩阵时候居然会爆LL！还要搞一个慢速乘处理一下。

此外就没了，上丑陋的代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
//1. n=2^k ---> f[k] get
//2. n%2==0 --> (n/2)%2==1 
//3. n%2==0 --> 公式     
//以上放屁
//以下为沙雕代码。。

#include <iostream>
#include <cstring>
#include <cmath>
#define ll long long
#define mo 1000000007
using namespace std;

int n,k;
ll b[4][4],p[4][4],e[4][4],kp[4][4],g[4][4],mul[4][4];

void ans()
{
    memset(kp,0,sizeof(kp));
    for (int i=1;i<=3;i++)
    for (int j=1;j<=3;j++)
    for (int k=1;k<=3;k++)
      kp[i][j]=(kp[i][j]%mo+(g[i][k]%mo*p[k][j]%mo+mo)%mo+mo)%mo;
    for (int i=1;i<=3;i++)
    for (int j=1;j<=3;j++) g[i][j]=kp[i][j]%mo;
}

void mi()
{
    memset(kp,0,sizeof(kp));
    for (int i=1;i<=3;i++)
    for (int j=1;j<=3;j++)
    for (int k=1;k<=3;k++)
      kp[i][j]=(kp[i][j]%mo+(p[i][k]%mo*p[k][j]%mo+mo)%mo+mo)%mo;
    for (int i=1;i<=3;i++)
    for (int j=1;j<=3;j++) p[i][j]=kp[i][j]%mo;
}

void sol()
{
    memset(kp,0,sizeof(kp));
    for (int i=1;i<=3;i++)
    for (int j=1;j<=3;j++)
    for (int k=1;k<=3;k++)
      kp[i][j]=(kp[i][j]%mo+(g[i][k]%mo*b[k][j]%mo+mo)%mo+mo)%mo;
    for (int i=1;i<=3;i++)
    for (int j=1;j<=3;j++) e[i][j]=kp[i][j]%mo;
}

ll qmul(int x,int y)
{
    int ans=0;
    while (y>0)
    {
        if (y%2) ans=(ans%mo+x%mo)%mo;
        x=(x+x)%mo; y/=2;
    }
    return ans;
}

int main()
{
    cin>>k>>n;
    //begin:
    b[1][1]=2; b[2][1]=k; b[3][1]=(qmul(k,k)-2)%mo;
    //jiasu:
    p[1][1]=p[2][1]=p[2][2]=p[1][3]=p[3][1]=0;
    p[1][2]=p[2][3]=1; p[3][2]=-1; p[3][3]=k;
    //jiasu begin
    for (int i=1;i<=3;i++) g[i][i]=1;
    n-=2;
    while (n>0)  //处理构造矩阵的快速幂 
    {
        if (n%2) ans();
        mi(); n/=2;
    }
    sol();
    cout<<e[3][1]<<endl;
}