Problem Description

As we all known , the Fibonacci series : F(0) = 1, F(1) = 1, F(N) = F(N - 1) + F(N - 2) (N >= 2).Now we define another kind of Fibonacci : A(0) = 1 , A(1) = 1 , A(N) = X * A(N - 1) + Y * A(N - 2) (N >= 2).And we want to Calculate S(N) , S(N) = A(0)² +A(1)²+……+A(n)².

 

Input

There are several test cases.

Each test case will contain three integers , N, X , Y .

N : 2<= N <= 2³¹ – 1

X : 2<= X <= 2³¹– 1

Y : 2<= Y <= 2³¹ – 1

 

Output

For each test case , output the answer of S(n).If the answer is too big , divide it by 10007 and give me the reminder.

 

Sample Input

2 1 1

3 2 3 

 

Sample Output

6

196

 

#include <iostream>

#include <cstdio>

#include <cstring>

#include <algorithm>

using namespace std;

const int mod=10007;

struct matrix

{

    long long ma[5][5];

};

matrix multi(matrix x,matrix y)//矩阵相乘

{

    matrix ans;

    memset(ans.ma,0,sizeof(ans.ma));

    for(int i=1;i<=4;i++)

    {

        for(int j=1;j<=4;j++)

        {

            if(x.ma[i][j])//稀疏矩阵优化

            for(int k=1;k<=4;k++)

            {

                ans.ma[i][k]=(ans.ma[i][k]+(x.ma[i][j]*y.ma[j][k])%mod)%mod;

            }

        }

    }

    return ans;

}

matrix pow(matrix a,long long m)

{

    matrix ans;

    for(int i=1;i<=4;i++)

    {

        for(int j=1;j<=4;j++)

        {

            if(i==j)

            ans.ma[i][j]=1;

            else

            ans.ma[i][j]=0;

        }

    }

    while(m)

    {

        if(m&1)

        ans=multi(ans,a);

        a=multi(a,a);

        m=m>>1;

    }

    return ans;

}

int main()

{

    long long x,y,n;

    while(~scanf("%I64d%I64d%I64d",&n,&x,&y))

    {

        matrix a,b;

        memset(a.ma,0,sizeof(a.ma));

        memset(b.ma,0,sizeof(b.ma));

        a.ma[1][1]=1;

        a.ma[1][2]=1;

        a.ma[2][2]=(x*x)%mod;

        a.ma[2][3]=(y*y)%mod;

        a.ma[2][4]=(2*x*y)%mod;

        a.ma[3][2]=1;

        a.ma[4][2]=x;

        a.ma[4][4]=y;

        b.ma[1][1]=1;

        b.ma[2][1]=1;

        b.ma[3][1]=1;

        b.ma[4][1]=1;

        a=pow(a,n);

        a=multi(a,b);

        printf("%I64d\n",a.ma[1][1]);

    }

    return 0;

}