POJ 3420 Quad Tiling

#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<iostream>
using namespace std;
typedef long long LL;
const double pi=acos(-1.0),eps=1e-8;
void File()
{
    freopen("D:\\in.txt","r",stdin);
    freopen("D:\\out.txt","w",stdout);
}
inline int read()
{
    char c = getchar();  while(!isdigit(c)) c = getchar();
    int x = 0;
    while(isdigit(c)) { x = x * 10 + c - ‘0‘; c = getchar(); }
    return x;
}

LL n,MOD;

struct Matrix
{
    long long A[18][18];
    int R, C;
    Matrix operator*(Matrix b);
};

Matrix X, Y, Z;

Matrix Matrix::operator*(Matrix b)
{
    Matrix c;
    memset(c.A, 0, sizeof(c.A));
    int i, j, k;
    for (i = 1; i <= R; i++)
        for (j = 1; j <= b.C; j++)
            for (k = 1; k <= C; k++)
                c.A[i][j] = (c.A[i][j] + (A[i][k] * b.A[k][j]) % MOD) % MOD;
    c.R = R; c.C = b.C;
    return c;
}

void dfs(int a,int b,int t,int c)
{
    if(t==4) { X.A[c+1][b+1]=1; return; }
    if(a&(1<<t)) dfs(a,b,t+1,c);
    else
    {
        dfs(a|(1<<t),b|(1<<t),t+1,c);
        if(t+1<4&&(a&(1<<(t+1)))==0) dfs(a+(1<<t)+(1<<(t+1)),b,t+1,c);
    }
}

void init()
{
    memset(X.A, 0, sizeof X.A);
    memset(Y.A, 0, sizeof Y.A);
    memset(Z.A, 0, sizeof Z.A);

    Y.R = 16; Y.C = 16;
    for (int i = 1; i <= 16; i++) Y.A[i][i] = 1;

    X.R = 16; X.C = 16;
    for(int i=0;i<=15;i++) dfs(i,0,0,i);

    Z.R = 1; Z.C = 16;
    Z.A[1][1]=1;
}

void work()
{
    while (n)
    {
        if (n % 2 == 1) Y = Y*X;
        n = n >> 1;
        X = X*X;
    }
    Z = Z*Y;
    printf("%lld\n",Z.A[1][16]);
}

int main()
{
    while(~scanf("%lld%lld",&n,&MOD))
    {
        if(n==0&&MOD==0) break;
        n++; init();
        work();
    }
    return 0;
}