HDU 3483 A Very Simple Problem

/*
*  [题意]
*   输入n, x, m
*   求(1^x)*(x^1)+(2^x)*(x^2)+(3^x)*(x^3)+...+(n^x)*(x^n)
*  [解题方法]
*   设f[n] = [x^n, n*(x^n), (n^2)*(x^n),..., (n^x)*(x^n)]
*   则f[n][k] = (n^k)*(x^n)
*   问题转化为求：( g[n] = f[1][x]+f[2][x]+...+f[n][x] )
*   设C(i,j)为组合数，即i种元素取j种的方法数
*   所以有：f[n+1][k] = ((n+1)^k)*(x^(n+1)) （二次多项式展开）
*                     = x*( C(k,0)*(x^n)  +C(k,1)*n*(x^n)+...+C(k,k)*(n^k)*(x^n) )
*                     = x*( C(k,0)*f[n][0]+C(k,1)*f[n][1]+...+C(k,k)*f[n][k] )
*   所以得：
*   |x*1 0................................0|        |f[n][0]|       |f[n+1][0]|
*   |x*1 x*1 0............................0|        |f[n][1]|       |f[n+1][1]|
*   |x*1 x*2 x*1 0........................0|    *   |f[n][2]|   =   |f[n+1][2]|
*   |......................................|        |.......|       |.........|
*   |x*1 x*C(k,1) x*C(k,2)...x*C(k,x) 0...0|        |f[n][k]|       |f[n+1][k]|
*   |......................................|        |.......|       |.........|
*   |x*1 x*C(x,1) x*C(x,2).......x*C(x,x) 0|        |f[n][x]|       |f[n+1][x]|
*   |0................................0 1 1|        |g[n-1] |       | g[ n ]  |
*/
#include <iostream>
#include <cstring>
using namespace std;
#define M 55
#define LL long long
#define FF(i, n) for(int i = 0; i < n; i++)

int ans[M], mod, C[M][M];
int ret[M][M], init[M][M];

void ini(int n, int x)
{
    memset(init, 0, sizeof(init));
    FF(i, n) FF(j, i+1)
        init[i][j] = (LL)x*C[i][j] % mod;
    FF(i, n) ans[i] = x;
    ans[n] = 0;
    init[n][n-1] = init[n][n] = 1;
}

void matmul(int a[][M], int b[][M], int n)
{
    int tp[M][M] = {0};
    FF(i, n) FF(k, n) if(a[i][k]) FF(j, n) if(b[k][j])
        tp[i][j] = (tp[i][j] + (LL)a[i][k]*b[k][j]) % mod;
    FF(i, n) FF(j, n) a[i][j] = tp[i][j];
}

void matmul(int a[], int b[][M], int n)
{
    int tp[M] = {0};
    FF(j, n) if(a[j]) FF(i, n) if(b[i][j])
        tp[i] = (tp[i] + (LL)a[j]*b[i][j]) % mod;
    FF(i, n) a[i] = tp[i];
}

void qmod(int n, int b)
{
    FF(i, n) FF(j, n) ret[i][j] = (i==j);
    for( ; b; b >>= 1)
    {
        if (b & 1) matmul(ret, init, n);
        matmul(init, init, n);
    }
}

int main()
{
    int i, j, n, x;
    while (cin >> n >> x >> mod, n >= 1)
    {
        for(i = 0; i <= x; i++)
            C[i][0] = C[i][i] = 1;
        for(i = 2; i <= x; i++)
            for(j = 1; j < i; j++)
                C[i][j] = ((LL)C[i-1][j-1]+C[i-1][j]) % mod;
        ini(x+1, x);
        qmod(x+2, n);
        matmul(ans, ret, x+2);
        cout << ans[x+1] << endl;
    }
    return 0;
}