【BZOJ4002】【JLOI2015】有意义的字符串（数学，矩阵快速幂）

看来我数学还是太弱了，这道题也没有想到。
我们发现$(\frac{b+\sqrt d}{2})^n$是个无理数，而且这个数非常大，所以无法直接计算。我们又发现$(\frac{b-\sqrt d}{2})^n$的值虽然是个无理数，但是他的值向下取整之后非常好计算。然后我们又发现，$(\frac{b+\sqrt d}{2})^n+(\frac{b-\sqrt d}{2})^n$是个有理数，又因为$b\space mod\space2=1,d\space mod\space4=1$，所以不难证明这个数是一个整数。
所以，我们可以将问题转变为求$\lfloor[(\frac{b+\sqrt d}{2})^n+(\frac{b-\sqrt d}{2})^n]-(\frac{b-\sqrt d}{2})^n\rfloor$
然后我们就可以把问题拆成两部分。
首先我们考虑第一部分。我们可以构造一个一元二次方程：$x^2-b+\frac{b^2-d}{4}=0$，那么它的解$x_1=\frac{b+\sqrt d}{2},x_2=\frac{b-\sqrt d}{2}$。
那么第一部分就等于$x_1^n+x_2^n$
然后我们把这个式子变形一下就可以得到：$x_1^n+x_2^n=(x_1+x_2)(x_1^{n-1}+x_2^{n-1})-x_1x_2(x_1^{n-2}+x_2^{n-2})$
然后令$F(n)=x_1^n+x_2^n$，
可以得到$F(n)=(x_1+x_2)F(n-1)-x_1x_2F(n-2)$
然后我们构造的一元二次方程就派上用场了！
根据韦达定理，对于一个一元二次方程$ax^2+bx+c=0$的两个解$x_{1,2}$满足：$x_1+x_2=-\frac{b}{a},x_1x_2=\frac{c}{a}$
所以$F(n)=bF(n-1)+\frac{d-b^2}{2}F(n-2)$
其中$F(0)=2,F(1)=b$
推到这里就是非常经典的的矩乘快速递推啦！
对于第二部分，
$\because b^2\leq d<(b+1)^2$
又$\because b,d\in N$
$\therefore b\leq\sqrt d< b+1$
$\therefore \frac{b-\sqrt d}{2}\in(-1,0]$
不难发现，当且仅当$n$为奇数且$\sqrt d!=b$时，这一部分对答案有$-1$的贡献。所以我们判一下是否减$1$就行了。
有一个坑点，这个模数非常大，加起来可能会爆$long\space long$，所以我们加的时候要判一下它们的和是否小于$0$。
注意两个数相乘会爆$long\space long$，所以要用快速乘。
快速乘其实主要就是利用$long\space double$的较高的精度，再利用我们小学学的一个东西（被除数减除数乘商等于余数）,就可以在$O(1)$的时间内完成较大的数的乘法并取模。
如果有错在评论区吼一声哦！
代码：

#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll mod=7528443412579576937ll;
ll b,d,n;
ll Add(ll x,ll y){
    return x+y>=mod||x+y<0?x+y-mod:x+y;
}
ll Minus(ll x,ll y){
    return x-y<0?x-y+mod:x-y;
}
ll Mul(ll x,ll y){
    ll d=(ll)(x*(long double)y/mod+0.5);
    return Minus(x*y,d*mod);
}
struct Matrix{
    ll a[2][2];
    ll *operator[](int b){
        return a[b];
    }
    friend Matrix operator*(Matrix A,Matrix B){
        Matrix ret;
        for(int i=0;i<2;i++)
            for(int j=0;j<2;j++){
                ret[i][j]=0;
                for(int k=0;k<2;k++)
                    ret[i][j]=Add(ret[i][j],Mul(A[i][k],B[k][j]));
            }
        return ret;
    }
    Matrix operator*=(Matrix B){
        *this=*this*B;
        return *this;
    }
}A,B;
Matrix qpow(Matrix x,ll n){
    Matrix ret;
    ret[0][0]=ret[1][1]=1,ret[0][1]=ret[1][0]=0;
    while(n){
        if(n&1ll)
            ret*=x;
        x*=x;
        n>>=1ll;
    }
    return ret;
}
int main(){
    scanf("%lld%lld%lld",&b,&d,&n);
    if(!n){
        puts("1");
        return 0;
    }
    A[0][0]=b;
    A[0][1]=2;
    B[0][0]=b;
    B[0][1]=1;
    B[1][0]=(d-b*b)/4;
    A[1][0]=A[1][1]=B[1][1]=0;
    A*=qpow(B,n-1);
    printf("%lld",Minus(A[0][0],!(n&1ll)&&(d!=b*b)));
    return 0;
}