矩阵连乘（3）

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const long long mod=10007;

typedef struct
{
  long long m[4][4];
}mat;

mat I={1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1};

mat calc(mat a,mat b)
{
    int i,j,k;
    mat c;
    for(i=0;i<4;i++)
    for(j=0;j<4;j++)
    {
        c.m[i][j]=0;
       for(k=0;k<4;k++)
      {
        c.m[i][j]+=(a.m[i][k]*b.m[k][j])%mod;
      }
      c.m[i][j]=c.m[i][j]%mod;
    }
    return c;
}

mat matirx(mat P,long long n)
{
    mat m=P,b=I;
    while(n>=1)
    {
        if(n&1) b=calc(b,m);
        n>>=1;
        m=calc(m,m);
    }
    return b;
}

int main()
{
    long long n,x,y;
    while(scanf("%lld%lld%lld",&n,&x,&y)!=EOF)
    {
        long long sum=0;
        x=x%mod;
        y=y%mod;   //2*x*y可能会溢出
        mat P={x*x,2*x*y,y*y,0,x,y,0,0,1,0,0,0,1,0,0,1};
        mat a;
        a=matirx(P,n);
        sum=(a.m[3][0]+a.m[3][1]+a.m[3][2]+a.m[3][3])%mod;
        printf("%lld\n",sum);
    }
    return 0;
}