[Project Euler 520]Simbers 题解

传送门：https://projecteuler.net/problem=520

一个假数位dp。

矩阵乘法调了一年。。。

刚开始以为是个倍增fwt，发现怎么算复杂度都不对。

后面发现好像所有奇数和偶数本质是一样的，除了前导0要容斥一下。

似乎可以推式子了。

枚举$t$个奇数要用，$f_{i,j,k}$表示前i位，有$j$个偶数出现次数为偶数，有$k$个奇数出现次数为奇数
考虑$f_{i,j,k}$的贡献：
$f_{i,j,k}*(5-j) o f_{i+1,j+1,k} $
$f_{i,j,k}*(t-k) o f_{i+1,j,k+1}$
$f_{i,j,k}*j o f_{i+1,j-1,k}$
$f_{i,j,k}*k o f_{i+1,j,k-1}$
我们发现这是个线性递推，矩阵维护即可

调了一年。

//waz
#include <bits/stdc++.h>

using namespace std;

#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define ALL(x) (x).begin(), (x).end()
#define SZ(x) ((int)((x).size()))

typedef pair<int, int> PII;
typedef vector<int> VI;
typedef long long int64;
typedef unsigned int uint;
typedef unsigned long long uint64;

#define gi(x) ((x) = F())
#define gii(x, y) (gi(x), gi(y))
#define giii(x, y, z) (gii(x, y), gi(z))

int F()
{
    char ch;
    int x, a;
    while (ch = getchar(), (ch < '0' || ch > '9') && ch != '-');
    if (ch == '-') ch = getchar(), a = -1;
    else a = 1;
    x = ch - '0';
    while (ch = getchar(), ch >= '0' && ch <= '9')
        x = (x << 1) + (x << 3) + ch - '0';
    return a * x;
}

const int mod = 1000000123;

int inc(int a, int b) { a += b; return a >= mod ? a - mod : a; }

int dec(int a, int b) { a -= b; return a < 0 ? a + mod : a; }

struct mat
{
    int a[50][50];
    int n;
    void clear() { n = 0; memset(a, 0, sizeof a); }
    friend mat operator * (const mat &a, const mat &b)
    {
        mat c;
        c.clear();
        c.n = a.n;
        for (int i = 1; i <= c.n; ++i)
            for (int j = 1; j <= c.n; ++j)
                for (int k = 1; k <= c.n; ++k)
                    c.a[i][j] = inc(c.a[i][j], 1LL * a.a[i][k] * b.a[k][j] % mod);
        return c; 
    }
    friend mat operator + (mat a, mat b)
    {
        for (int i = 1; i <= a.n; ++i)
            for (int j = 1; j <= b.n; ++j)
                a.a[i][j] = inc(a.a[i][j], b.a[i][j]);
        return a;
    }
    friend mat operator - (mat a, mat b)
    {
        for (int i = 1; i <= a.n; ++i)
            for (int j = 1; j <= b.n; ++j)
                a.a[i][j] = dec(a.a[i][j], b.a[i][j]);
        return a;
    }
} g;

int work(int x, int p, int t, int gg)
{
    static int id[10][10]; int idt = 0;
    for (int i = 0; i <= p; ++i)
        for (int j = 0; j <= t; ++j)
            id[i][j] = ++idt;
    g.clear();
    g.n = idt;
    for (int j = 0; j <= p; ++j)
        for (int k = 0; k <= t; ++k)
        {
            if (j + 1 <= p)
            {
                int x = id[j][k], y = id[j + 1][k];
                g.a[x][y] = p - j;
            }
            if (k + 1 <= t)
            {
                int x = id[j][k], y = id[j][k + 1];
                g.a[x][y] = t - k;
            }
            if (j - 1 >= 0)
            {
                int x = id[j][k], y = id[j - 1][k];
                g.a[x][y] = j;
            }
            if (k - 1 >= 0)
            {
                int x = id[j][k], y = id[j][k - 1];
                g.a[x][y] = k;
            }
        }
    /*for (int i = 1; i <= idt; ++i)
        for (int j = 1; j <= idt; ++j)
            printf("%d%c", g.a[i][j], " 
"[j == idt]);
    cerr << id[p][0] << ", " << id[p][t] << endl;*/
    mat ret;
    ret.clear();
    ret.n = idt;
    ret.a[1][id[p][gg]] = 1;
    for (; x; x >>= 1)
    {
        if (x & 1) ret = ret * g;
        g = g * g;
    }
    return ret.a[1][id[p][t]];
}

int solve1(int x, int p, int t, int gg)
{
    static int id[10][10]; int idt = 0;
    for (int i = 0; i <= p; ++i)
        for (int j = 0; j <= t; ++j)
            id[i][j] = ++idt;
    g.clear();
    g.n = idt;
    for (int j = 0; j <= p; ++j)
        for (int k = 0; k <= t; ++k)
        {
            if (j + 1 <= p)
            {
                int x = id[j][k], y = id[j + 1][k];
                g.a[x][y] = p - j;
            }
            if (k + 1 <= t)
            {
                int x = id[j][k], y = id[j][k + 1];
                g.a[x][y] = t - k;
            }
            if (j - 1 >= 0)
            {
                int x = id[j][k], y = id[j - 1][k];
                g.a[x][y] = j;
            }
            if (k - 1 >= 0)
            {
                int x = id[j][k], y = id[j][k - 1];
                g.a[x][y] = k;
            }
        }
    mat ret, c;
    ret.clear();
    c.clear();
    c.n = idt;
    ret.n = idt;
    ret.a[1][id[p][gg]] = 1;
    static mat f[50], tt[50];
    f[0] = g;
    tt[0] = g;
    for (int i = 1; i <= x; ++i)
    {
        f[i].clear();
        f[i] = f[i - 1] * f[i - 1];
        tt[i] = f[i - 1] * tt[i - 1] + tt[i - 1];
    }
    for (int i = 0; i <= x; ++i) 
        c = c + tt[i];
    ret = ret * c;
    return ret.a[1][id[p][t]];
}

int solve2(int x, int p, int t, int gg)
{
    static int id[10][10]; int idt = 0;
    for (int i = 0; i <= p; ++i)
        for (int j = 0; j <= t; ++j)
            id[i][j] = ++idt;
    g.clear();
    g.n = idt;
    for (int j = 0; j <= p; ++j)
        for (int k = 0; k <= t; ++k)
        {
            if (j + 1 <= p)
            {
                int x = id[j][k], y = id[j + 1][k];
                g.a[x][y] = p - j;
            }
            if (k + 1 <= t)
            {
                int x = id[j][k], y = id[j][k + 1];
                g.a[x][y] = t - k;
            }
            if (j - 1 >= 0)
            {
                int x = id[j][k], y = id[j - 1][k];
                g.a[x][y] = j;
            }
            if (k - 1 >= 0)
            {
                int x = id[j][k], y = id[j][k - 1];
                g.a[x][y] = k;
            }
        }
    mat ret, c;
    ret.clear();
    c.clear();
    c.n = idt;
    ret.n = idt;
    ret.a[1][id[p][gg]] = 1;
    static mat f[50], tt[50], pp[50];
    f[0] = g;
    tt[0] = g;
    pp[0].clear(); pp[0].n = idt;
    for (int i = 1; i <= x; ++i)
    {
        f[i].clear();
        f[i] = f[i - 1] * f[i - 1];
        tt[i] = f[i - 1] * tt[i - 1] + tt[i - 1];
        pp[i] = tt[i] - f[i];
    }
    for (int i = 0; i <= x; ++i) 
        c = c + pp[i];
    ret = ret * c;
    return ret.a[1][id[p][t]];
}

int C[110][110];

int work(int m)
{
    int ans = 0;
    for (int n = 1; n <= m; ++n)
    for (int t = 0; t <= 5; ++t)
    {
        int x = dec(work(n, 5, t, 0), work(n - 1, 4, t + 1, 0));
        ans = (ans + 1LL * x * C[5][t]) % mod;
    }
    return ans;
}

int main()
{
    for (int i = *C[0] = 1; i <= 100; ++i)
        for (int j = C[i][0] = 1; j <= i; ++j)
            C[i][j] = inc(C[i - 1][j - 1], C[i - 1][j]);
    int m;
    int c = 39;
    /*int ret = 0;
    for (int i = 0; i <= 3; ++i)
    {
        ret = inc(ret, work(1 << i));
    }
    cerr << ret << endl;*/
    int ans = 0;
    for (int t = 0; t <= 5; ++t) ans = (ans + 1LL * solve1(c, 5, t, 0) * C[5][t]) % mod;
    for (int t = 0; t <= 5; ++t) ans = (ans - 1LL * solve2(c, 4, t + 1, 0) * C[5][t]) % mod;
    ans = dec(ans, work(1));
    cout << ans << endl;
    return 0;
}