結果

問題 No.978 Fibonacci Convolution Easy
ユーザー tarattata1
提出日時 2019-07-14 00:04:44
言語 C++11（廃止可能性あり）
(gcc 13.3.0)
結果
AC  
実行時間 1 ms / 2,000 ms
コード長 3,830 bytes
コンパイル時間 438 ms
コンパイル使用メモリ 44,544 KB
実行使用メモリ 5,376 KB
最終ジャッジ日時 2024-09-18 20:50:28
合計ジャッジ時間 1,143 ms
ジャッジサーバーID
（参考情報） 		judge1 / judge2
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 21
権限があれば一括ダウンロードができます
コンパイルメッセージ
main.cpp: In function ‘int main(int, char**)’:
main.cpp:119:10: warning: ignoring return value of ‘int scanf(const char*, ...)’ declared with attribute ‘warn_unused_result’ [-Wunused-result]
  119 |     scanf("%d%d", &n, &p);
      |     ~~~~~^~~~~~~~~~~~~~~~

ソースコード

diff # 

#include <stdio.h>
#include <vector>
#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
typedef pair<ll, ll> QQ;   // {a+b√q)}

QQ ADD(QQ a, QQ b) {
    return make_pair( (a.first + b.first)%MOD, (a.second + b.second)%MOD );
}

QQ SUB(QQ a, QQ b) {
    return make_pair( (a.first + MOD - b.first)%MOD, (a.second + MOD - b.second)%MOD );
}

QQ PROD(QQ a, QQ b) {
    return make_pair((a.first * b.first %MOD + (a.second * b.second)%MOD * q %MOD )%MOD,
                     (a.first * b.second %MOD + a.second * b.first %MOD)%MOD );
}

QQ INV(QQ x) {
    ll det=(x.first*x.first%MOD+(MOD-(x.second*x.second%MOD)*q%MOD))%MOD;
    ll tmp=minv(det);
    return make_pair(x.first*tmp%MOD, (MOD-x.second)*tmp%MOD);
}

QQ POW(QQ x, ll n)
{
    QQ ans = make_pair(1, 0);
    while(n != 0){
        if(n&1) ans = PROD(ans, x);
        x = PROD(x,x);
        n = n >> 1;
    }
    return ans;
}

ll solve(int n, int p)
{
    q = ((ll)p*p + 4)%MOD;

    QQ val1 = make_pair(1,0);  // 1
    ll inv2 = minv(2);      // 1/2
    QQ a1=make_pair(0, 0);
    QQ a2=make_pair(1, 0);
    QQ aa=make_pair(p*inv2%MOD, 1*inv2%MOD);         // α
    QQ bb=make_pair(p*inv2%MOD, (MOD-1)*inv2%MOD);   // β
    QQ C0 = SUB(a2, PROD(bb,a1));         // a[2]-βa[1]
    QQ C  = PROD(C0, INV(SUB(aa, bb)));   // C
    QQ D0 = SUB(a2, PROD(aa,a1));         // (a[2]-αa[1])
    QQ D  = PROD(D0, INV(SUB(bb, aa)));   // D

    QQ suma = PROD(SUB(val1, POW(aa, n)), INV(SUB(val1, aa)));   // Σ(k=1..n-1)(α^k)
    QQ sumb = PROD(SUB(val1, POW(bb, n)), INV(SUB(val1, bb)));   // Σ(k=1..n-1)(β^k)
    QQ aa2=PROD(aa,aa);
    QQ suma2 = PROD(SUB(val1, POW(aa2, n)), INV(SUB(val1, aa2)));   // Σ(k=1..n-1)((α^2)^k)
    QQ bb2=PROD(bb,bb);
    QQ sumb2 = PROD(SUB(val1, POW(bb2, n)), INV(SUB(val1, bb2)));   // Σ(k=1..n-1)((β^2)^k))
    QQ ab2 = PROD(aa,bb);
    QQ sumab = PROD(SUB(val1, POW(ab2, n)), INV(SUB(val1, ab2)));    //Σ(k=1..n-1)((αβ)^k)
    QQ ans0 = ADD(PROD(C,suma), PROD(D,sumb));
    QQ ans1 = ADD(ADD(ADD(PROD(PROD(C,C),suma2), PROD(PROD(C,D),sumab)), PROD(PROD(C,D),sumab)), PROD(PROD(D,D),sumb2));
    QQ ans = ADD(PROD(ans0,ans0), ans1);
    ans = make_pair(ans.first*inv2%MOD, ans.second*inv2%MOD);

    return ans.first;
}

int main(int argc, char* argv[])
{
    int n;
    scanf("%d%d", &n, &p);
    printf("%lld\n", solve(n,p));
    return 0;
}
0