#include #include #pragma warning(disable:4996) typedef long long ll; const long long MOD = 1000000007; using namespace std; ll mpow(ll x, ll n){ //x^n(mod M) ll ans = 1; while(n != 0){ if(n&1) ans = ans*x % MOD; x = x*x % MOD; n = n >> 1; } return ans; } ll minv(ll x){ return mpow( x, MOD-2 ); } // <解法> // 求める答えは、(Σa[n])^2 + (Σ(a[n]^2))/2 // この計算は、O(n)ならば簡単に求めることができるが、 // nが大きい場合について解くには何か工夫が必要。 // 今回は、a[n]の一般項を求めて計算することにした。 // // x^2-px-1=0 の解は x=(p±√(p^2+4))/2 // これをα,βとおく。 // // a[n]-αa[n-1]=β(a[n-1]-αa[n-2]) // から // a[n]-αa[n-1]=β^(n-2) (a[2]-αa[1]) // 同様に // a[n]-βa[n-1]=α^(n-2) (a[2]-βa[1]) // 引き算して // (β-α)a[n-1]=β^(n-2) (a[2]-αa[1])-α^(n-1) (a[2]-βa[1]) // n-1 を n に書き換えると // a[n]=(β^(n-1))(a[2]-αa[1])-(α^(n-1))(a[2]-βa[1]))/(β-α) // これで // a[n]=C(α^(n-1))+D(β^(n-1)) // という形の一般項になったので、 // Σa[n] = CΣ(α^(n-1))+DΣ(β^(n-1)) // Σ(a[n]^2) = (C^2)Σ((α^2)^(n-1))+(2CD)Σ((αβ)^(n-1))+(D^2)Σ((β^2)^(n-1)) // ともに等比級数の和の公式から求めることができる // // α,βを使って実際に計算する際には、今回は実数による計算ではなく、 // {a+b√(p^2+4)}(pは固定、a,bは整数mod1000000007)という体の上で計算する int p; int q; // q=p^2+4 struct Q // {f+s√q} { ll f; // first ll s; // second public: Q(ll f0, ll s0): f(f0),s(s0){} Q(): f(0),s(0){} Q(ll f0): f(f0),s(0){} Q operator + (Q r) { return Q( (this->f + r.f)%MOD, (this->s + r.s)%MOD ); } Q operator - (Q r) { return Q( (this->f + MOD - r.f)%MOD, (this->s + MOD - r.s)%MOD ); } Q operator * (Q r) { return Q((this->f * r.f + (this->s * r.s)%MOD * q)%MOD, (this->f * r.s + this->s * r.f )%MOD ); } Q operator / (Q r) { ll det=(r.f*r.f+(MOD-(r.s*r.s%MOD)*q%MOD))%MOD; ll tmp=minv(det); return ((*this) * Q(r.f*tmp%MOD, (MOD-r.s)*tmp%MOD)); } }; Q POW (Q x, ll n) { Q ans = 1; while(n != 0){ if(n&1) ans = ans * x; x = x * x; n = n >> 1; } return ans; } int main(int argc, char* argv[]) { int n; scanf("%d%d", &n, &p); q = ((ll)p*p + 4)%MOD; ll inv2 = minv(2); // 1/2 Q a1=0; Q a2=1; Q aa=Q(p*inv2%MOD, 1*inv2%MOD); // α Q bb=Q(p*inv2%MOD, (MOD-1)*inv2%MOD); // β Q C = (a2 - bb*a1) / (aa - bb); // C Q D = (a2 - aa*a1) / (bb - aa); // D Q suma = (Q(1) - POW(aa,n)) / (Q(1) - aa); // Σ(k=0..n-1)(α^k) Q sumb = (Q(1) - POW(bb,n)) / (Q(1) - bb); // Σ(k=0..n-1)(β^k) Q suma2= (Q(1) - POW(aa*aa,n)) / (Q(1) - aa*aa); // Σ(k=0..n-1)((α^2)^k) Q sumb2= (Q(1) - POW(bb*bb,n)) / (Q(1) - bb*bb); // Σ(k=0..n-1)((β^2)^k)) Q sumab= (Q(1) - POW(aa*bb,n)) / (Q(1) - aa*bb); // Σ(k=0..n-1)((αβ)^k) Q ans0 = C * suma + D * sumb; Q ans1 = C * C * suma2 + C * D * 2 * sumab + D * D* sumb2; Q ans = ans0 * ans0 + ans1; ans = ans * inv2; printf("%lld\n", ans); return 0; }