F. 雷方块¶

分数：200分

向着星辰...... 欢迎来到......

你向往常一样领取派遣探索和每日任务奖励，但今天凯瑟琳给了你一件特别的委托。

你看了一眼任务说明，只需要完成一个雷方块解密就能拿到1600原石，想着已经75发大保底了，做完就能看到奶奶在壶里脸红的样子，立刻按照指引传送到天云山上下顶端的锚点。原本雷鸟所在的平台已经密密麻麻铺上了雷方块，你扔了发天动万象，毫无头绪。

任务¶

给定一个由 0，1，2，3 这四个数构成的矩阵，代表由三种状态的方块构成的阵列，其中 0 表示空缺，1，2，3 分别为三种状态。敲击一个方块会使这个方块和相邻四个位置中的方块按照 (1, 2, 3) 的顺序轮回。求解每个方块需要敲击多少次能让全部方块状态变成3。

注意¶

注意到一个方块连续敲三次整体状态还原。题目保证在每个方块敲击次数不超过2的情况下具有唯一解。

输入¶

一个二进制文件 in.data ，文件内容为

2 个 32 位整数，n1 和 n2，保证 $n1 \leq n2$
一个 n2 行 n1 列的 32 位整数矩阵，按行优先方式存储，每个数表示这个位置是否有方块以及方块的初始状态。

样例输入¶

4, 4
0, 1, 1, 2
1, 3, 0, 3
0, 3, 1, 1
3, 1, 2, 2

二进制形态：

4 4 0 1 1 2 1 3 0 3 0 3 1 1 3 1 2 2

输出¶

一个二进制文件 out.data ，文件内容为一个 n2 行 n1 列的 32 位整数矩阵，按行优先方式存储，每个数表示这个位置的方块需要敲击的次数。

额外要求¶

为了保证输出唯一，额外提出以下两点要求。

敲击次数必须大于等于 0，小于等于 2。
如果一个位置没有方块，则矩阵中这个位置的整数必须为 0。

样例输出¶

0, 1, 1, 0
2, 0, 0, 0
0, 0, 0, 0
0, 0, 2, 2

附件¶

本题没有附件，如果需要测试数据可以使用以下代码生成测试数据，需要自行挑出有唯一解的情况

import numpy
n1=16
n2=16
x=numpy.zeros(2+n1*n2,dtype=numpy.int32)
x[0]=n1
x[1]=n2
r=numpy.random.rand(n1*n2)
x[2:][r>0.1]=1
x[2:][r>0.4]=2
x[2:][r>0.7]=3
x.tofile('in.data')

可以使用以下代码校验结果

import numpy
x=numpy.fromfile('in.data',dtype=numpy.int32)
n1=x[0]
n2=x[1]
m=x[2:].reshape((n2,n1))

m_=numpy.fromfile('out.data',dtype=numpy.int32).reshape((n2,n1))
assert((m_<=2).all()) # 敲3次等于没敲
assert((m_>=0).all()) # 怎么想都不能敲-1次
m__=numpy.array(m_)
m__[:,1:]+=m_[:,:-1]
m__[:,:-1]+=m_[:,1:]
m__[1:,:]+=m_[:-1,:]
m__[:-1,:]+=m_[1:,:]
assert(((m+m__)%3*m==0).all()) # 所有方块状态为3
assert(((m==0)*m_==0).all())   # 没有方块的位置必须是0

提示¶

本题写成模3下的线性方程组后可以使用直接求解的方式，包括高斯消元和LU分解，大致过程如下：

import numpy
x=numpy.fromfile('in.data',dtype=numpy.int32)
n1=x[0]
n2=x[1]
m=x[2:].reshape((n2,n1))
im=numpy.zeros((n2,n1),dtype=numpy.int32)

ci=0
for i in range(n2):
    for j in range(n1):
        if m[i,j]:
            im[i,j]=ci
            ci+=1
        else:
            im[i,j]=-1
nl=((0,1),(0,-1),(1,0),(-1,0),(0,0))

import sage.all
t=sage.all.Integers(3)

a=sage.all.matrix(t,ci,ci)
y=sage.all.vector(t,ci)
for i in range(n2):
    for j in range(n1):
        ci=im[i,j]
        if ci>=0:
            y[ci]=3-m[i,j]
            for nl_ in nl:
                i_=i+nl_[0]
                j_=j+nl_[1]
                if 0<=i_ and i_<n2 and 0<=j_ and j_<n1:
                    ci_=im[i_,j_]
                    if ci_>=0:
                        a[ci_,ci]=1

x_t=a.solve_right(y)

x=numpy.zeros((n2,n1),dtype=numpy.int32)
for i in range(n2):
    for j in range(n1):
        ci=im[i,j]
        if ci>=0:
            x[i,j]=x_t[ci]
x.tofile('out.data')

复杂度为 n1^3*n2 的朴素高斯消元算法确能拿到满分。

n方可能过百万，但n方过百万有点不可能

评分标准¶

提交一个 cpp 文件，程序使用 g++ -O3 -march=native -fopenmp 编译，在 2 个核心上运行，设置环境变量 OMP_NUM_THREADS=2，运行脚本：

g++ -O3 -march=native -fopenmp program.cpp -o program
time OMP_NUM_THREADS=2 ./program

共 3 个测试点，每个测试点得到正确结果获得基本分数，性能分数与运行时间倒数成正比。

final_score = correct ? base_score + perf_score * min(goal / time, 1) : 0

	n1	n2	empty rate	base score	perf score	limit	goal
1	4	4	18.75%	10	0	1 s	1 s
2	512	512	9.96%	5	90	2 min	1.35 s
3	512	1024	9.96%	5	90	5 min	2.5 s

评测环境¶

容器镜像：gcc 。已安装 g++ 和 openmp。