頭がゴチャゴチャしてきて、詰め切るのがとても大変だった。
でも整理し切った後に出てくる性質は、至極単純なものだった。面白い。

問題概要

$H \times W$ のグリッドが与えられ、各マスは白黒に塗られている。この盤面をレベル 1 のフラクタルとする。

レベル $0$ のフラクタルとは黒いマスひとつからなる $1 \times 1$ のグリッド
レベル $k+1$ のフラクタルとは、レベル $k$ のフラクタルを、グリッドの全ての黒いマスに相当する位置に並べ、白いマスに相当する位置は白いマスで埋めたもの

レベル $K$ のフラクタルにおいて黒マスからなる連結成分の個数が何個あるか、1000000007 で割ったあまりを求めよ。

制約

$1 \le H, W \le 1000$
$0 \le K \le 10^{18}$

考えたこと

たとえば

...#.
#.###
###..
.#...
.###.

とかだと、実は何回レベルを引き上げても全部繋がっている。そうなる理由は、

列をみたときに 2 列目の両端が繋がっている
行をみたときに 3 行目の両端が繋がっている

からだ。一般に、両端の繋がっている行があって、両端の繋がっている行もあるとき、全部繋がる。

それは手を動かしてるとなんとなくわかるのだけど、数学的帰納法で示せる。もしレベル $k$ のフラクタルの連結成分が 1 個だけならば、それを上の形につないだとき、横方向にも縦方向にも全部つながるので、レベル $k+1$ のフラクタルも結局全部繋がるのだ。

(ここまで至るのが結構大変だった。頭が混乱したポイントとして「横一列すべて黒マス」と「横一列の両端が黒マス」との区別だった。結局前者は全然特徴量として意味がないことと思い至るまでが長かった)

さらに、「どの行も両端の少なくとも一方は白マス」かつ「どの列も両端の少なくとも一方は白マス」となっている場合もわかりやすい。この場合もまた、レベル $k$ のフラクタルを成長させたとき、それらは横方向にも縦方向にも繋がらない。よって、黒マスの個数を $B$ として

$B^{K-1}$ 個

となる。

縦にはつながらず、横にはつながる場合

最後に、

どの列を見ても、両端のうちのどちらかは白マス (縦にはつながらない)
ある列が存在して、両端が黒マス (横にはつながる)

という場合を考える。結論からいうと、レベル 1 のフラクタル (与えられた入力) において

黒マスの個数 = $B$
両端が黒マスとなっている行の本数 = $D$
黒マスが横方向に隣接している箇所の個数 = $A$

がパラメータとなる。そしてレベル $k$ において

連結成分の個数を $x_{k}$
レベル $k$ のフラクタルを横につないだときに減少する連結成分の個数 $y_{k}$

としたとき、以下のような漸化式が立てられる:

$x_{k+1} = Bx_{k} - Ay_{k}$
$y_{k+1} = Dy_{k}$

よって行列累乗で解くことができる。なお、初期条件として、 $x_{1} = 1$ , $y_{1} = 1$ が成立する。計算量は $O(\log{K} + HW)$ 。

...という具合に、結論はめちゃくちゃシンプルなのだが、最初のうちは

両端が黒マスとなっている行自体が縦方向に隣接していたらヤバイんじゃないか...
両端が黒マスとなっている行につき、連結成分が 1 減りそうだけど、実は奥で繋がったりして、変なダブルカウントをしないだろうか...

とか考えて頭がパニックだった。落ち着いて整理すると、元々フラクタルを成長させるときに縦方向には繋がらないことを仮定しているので関係ない。

#include <iostream>
#include <vector>
using namespace std;


// matrix
template<class T> struct Matrix {
    vector<vector<T> > val;
    Matrix(int n = 1, int m = 1, T v = 0) : val(n, vector<T>(m, v)) {}
    void init(int n, int m, T v = 0) {val.assign(n, vector<T>(m, v));}
    void resize(int n, int m) {
        val.resize(n);
        for (int i = 0; i < n; ++i) val[i].resize(m);
    }
    Matrix<T>& operator = (const Matrix<T> &A) {
        val = A.val;
        return *this;
    }
    size_t size() const {return val.size();}
    vector<T>& operator [] (int i) {return val[i];}
    const vector<T>& operator [] (int i) const {return val[i];}
    friend ostream& operator << (ostream& s, const Matrix<T>& M) {
        s << endl;
        for (int i = 0; i < (int)M.size(); ++i) s << M[i] << endl;
        return s;
    }
};

template<class T> Matrix<T> operator * (const Matrix<T> &A, const Matrix<T> &B) {
    Matrix<T> R(A.size(), B[0].size());
    for (int i = 0; i < A.size(); ++i)
        for (int j = 0; j < B[0].size(); ++j)
            for (int k = 0; k < B.size(); ++k)
                R[i][j] += A[i][k] * B[k][j];
    return R;
}

template<class T> Matrix<T> pow(const Matrix<T> &A, long long n) {
    Matrix<T> R(A.size(), A.size());
    auto B = A;
    for (int i = 0; i < A.size(); ++i) R[i][i] = 1;
    while (n > 0) {
        if (n & 1) R = R * B;
        B = B * B;
        n >>= 1;
    }
    return R;
}

template<class T> vector<T> operator * (const Matrix<T> &A, const vector<T> &B) {
    vector<T> v(A.size());
    for (int i = 0; i < A.size(); ++i)
        for (int k = 0; k < B.size(); ++k)
            v[i] += A[i][k] * B[k];
    return v;
}

template<class T> Matrix<T> operator + (const Matrix<T> &A, const Matrix<T> &B) {
    Matrix<T> R(A.size(), A[0].size());
    for (int i = 0; i < A.size(); ++i)
        for (int j = 0; j < A[0].size(); ++j)
            R[i][j] = A[i][j] + B[i][j];
    return R;
}

template<class T> Matrix<T> operator - (const Matrix<T> &A, const Matrix<T> &B) {
    Matrix<T> R(A.size(), A[0].size());
    for (int i = 0; i < A.size(); ++i)
        for (int j = 0; j < A[0].size(); ++j)
            R[i][j] = A[i][j] - B[i][j];
    return R;
}

// modint
template<int MOD> struct Fp {
    long long val;
    constexpr Fp(long long v = 0) noexcept : val(v % MOD) {
        if (val < 0) val += MOD;
    }
    constexpr int getmod() { return MOD; }
    constexpr Fp operator - () const noexcept {
        return val ? MOD - val : 0;
    }
    constexpr Fp operator + (const Fp& r) const noexcept { return Fp(*this) += r; }
    constexpr Fp operator - (const Fp& r) const noexcept { return Fp(*this) -= r; }
    constexpr Fp operator * (const Fp& r) const noexcept { return Fp(*this) *= r; }
    constexpr Fp operator / (const Fp& r) const noexcept { return Fp(*this) /= r; }
    constexpr Fp& operator += (const Fp& r) noexcept {
        val += r.val;
        if (val >= MOD) val -= MOD;
        return *this;
    }
    constexpr Fp& operator -= (const Fp& r) noexcept {
        val -= r.val;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr Fp& operator *= (const Fp& r) noexcept {
        val = val * r.val % MOD;
        return *this;
    }
    constexpr Fp& operator /= (const Fp& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b; swap(a, b);
            u -= t * v; swap(u, v);
        }
        val = val * u % MOD;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr bool operator == (const Fp& r) const noexcept {
        return this->val == r.val;
    }
    constexpr bool operator != (const Fp& r) const noexcept {
        return this->val != r.val;
    }
    friend constexpr ostream& operator << (ostream &os, const Fp<MOD>& x) noexcept {
        return os << x.val;
    }
    friend constexpr Fp<MOD> modpow(const Fp<MOD> &a, long long n) noexcept {
        if (n == 0) return 1;
        auto t = modpow(a, n / 2);
        t = t * t;
        if (n & 1) t = t * a;
        return t;
    }
};

const int MOD = 1000000007;
using mint = Fp<MOD>;
long long H, W, K;
vector<string> fi;

mint solve() {
    if (K <= 1) return 1;
    long long N = 0, yoko = 0, tate = 0, adj = 0;
    for (int i = 0; i < H; ++i)
        for (int j = 0; j < W; ++j)
            if (fi[i][j] == '#') ++N;
    for (int i = 0; i < H; ++i)
        if (fi[i][0] == '#' && fi[i].back() == '#') ++yoko;
    for (int j = 0; j < W; ++j)
        if (fi[0][j] == '#' && fi[H-1][j] == '#') ++tate;
    if (yoko && tate) return 1;
    if (!yoko && !tate) return modpow(mint(N), K-1);

    Matrix<mint> A(2, 2, 0);
    A[0][0] = N;
    if (yoko) {
        for (int i = 0; i < H; ++i)
            for (int j = 0; j + 1 < W; ++j)
                if (fi[i][j] == '#' && fi[i][j+1] == '#') ++adj;
        A[0][1] = -adj;
        A[1][1] = yoko;
    }
    else {
        for (int i = 0; i + 1 < H; ++i)
            for (int j = 0; j < W; ++j)
                if (fi[i][j] == '#' && fi[i+1][j] == '#') ++adj;
        A[0][1] = -adj;
        A[1][1] = tate;
    }
    auto res = pow(A, K-1);
    return res[0][0] + res[0][1];
}

int main() {
    cin >> H >> W >> K;
    fi.resize(H);
    for (int i = 0; i < H; ++i) cin >> fi[i];
    cout << solve() << endl;
}