1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
#include <bits/stdc++.h>
using namespace std;

using U64 = unsigned long long;
template <typename T> using Mat = vector<vector<T>>;
template <typename T> T load() {
	T r;
	cin >> r;
	return r;
}
template <typename T> Mat<T> loadMatrix(int width, int height) {
	auto mat = Mat<T>(width, vector<T>(height));
	for (int y=0; y<height; ++y)
		for (int x=0; x<width; ++x)
			mat[x][y] = load<T>();
	return mat;
}

template <typename T> T qpow(T x, unsigned e) {
	if (e == 0) return T(1);
	if (e % 2 == 1) return qpow(x, e-1) * x;
	auto sub = qpow(x, e/2);
	return sub*sub;
}

template <typename T> struct V2 {
	T x, y;
};
template <U64 mod> class Mint {
	U64 val;
	static_assert(mod <= 4294967295, "Too big modulo value for U64-based Mint");
	public:
	explicit Mint(U64 raw):val(raw % mod){}
	Mint operator+(Mint other) const {
		return Mint(val + other.val);
	}
	Mint operator-(Mint other) const {
		return Mint(mod + val - other.val);
	}
	Mint operator*(Mint other) const {
		return Mint(val * other.val);
	}
	Mint operator/(Mint other) const {
		static_assert(mod >= 3, "Too small modulo value for Mint division");
		return (*this) * qpow(other, mod-2);
	}
	Mint& operator++() {
		return *this = *this + Mint(1);
	}
	unsigned long long value() const {
		return val;
	}
};

template <typename F> inline void vonNeumannNeighbourhood(int x, int y, int n, F f) {
	if (x > 0) f(x-1, y);
	if (x < n-1) f(x+1, y);
	if (y > 0) f(x, y-1);
	if (y < n-1) f(x, y+1);
}

vector<V2<int>> computer0Layer(int n, const Mat<char>& mat) {
	auto results = vector<V2<int>>();
	for (int x=0; x<n; ++x)
		for (int y=0; y<n; ++y)
			if (mat[x][y] == '#')
				results.push_back(V2<int>{ x, y });
	return results;
}
template <typename F> void layerBFS(int n, const Mat<char>& mat, F f) {
	auto visited = Mat<bool>(n, vector<bool>(n, false));
	auto curr = computer0Layer(n, mat);
	auto next = vector<V2<int>>();
	auto layer = 0;
	for (auto zero : curr)
		visited[zero.x][zero.y] = true;
	while (not curr.empty()) {
		for (auto vertex : curr) {
			f(vertex.x, vertex.y, layer);
			vonNeumannNeighbourhood(vertex.x, vertex.y, n, [&](int kidx, int kidy){
				if (not visited[kidx][kidy]) {
					visited[kidx][kidy] = true;
					next.push_back(V2<int>{ kidx, kidy });
				}
			});
		}
		curr = move(next);
		++layer;
	}
}
Mat<int> computeLayers(int n, const Mat<char>& mat) {
	auto layers = Mat<int>(n, vector<int>(n, -1));
	layerBFS(n, mat, [&](int x, int y, int layer){
		layers[x][y] = layer;
	});
	return layers;
}

using M = Mint<1000000007>;

template <typename T> inline int countMatrix(const Mat<T>& mat, const T& val) {
	auto result = 0;
	for (auto&& column : mat)
		for (auto&& cell : column)
			if (cell == val)
				++result;
	return result;
}
template <typename T, typename F> inline void anchorAt(const Mat<T>& mat, const T& val, F f) {
	auto n = static_cast<int>(mat.size());
	for (int x=0; x<n; ++x)
		for (int y=0; y<n; ++y)
			if (mat[x][y] == val)
				f(x, y);
}
template <typename T, typename F> inline void nearby(int x, int y, const Mat<T>& mat, const T& val, F f) {
	vonNeumannNeighbourhood(x, y, mat.size(), [&](int x1, int y1){
		if (mat[x1][y1] == val)
			f(x1, y1);
	});
}

M solveK1(int n, const Mat<char>& mat) {
	auto layers = computeLayers(n, mat);
	return M(countMatrix(layers, 1));
}
M solveK2(int n, const Mat<char>& mat) {
	auto layers = computeLayers(n, mat);
	auto onesCount = M(0);
	auto strain12Count = M(0);
	anchorAt(layers, 1, [&](int x, int y){
		++onesCount;
		nearby(x, y, layers, 2, [&](int, int){
			++strain12Count;
		});
	});
	return onesCount * (onesCount - M(1)) / M(2) + strain12Count;
}

M solve(int n, int k, const Mat<char>& mat) {
	if (k == 1) {
		return solveK1(n, mat);
	} else if (k == 2) {
		return solveK2(n, mat);
	} else {
		throw runtime_error("unsupported case :(");
	}
}

int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	auto n = load<int>();
	auto k = load<int>();
	auto mat = loadMatrix<char>(n, n);
	auto result = solve(n, k, mat);
	cout << result.value() << '\n';
}