在每行和每列中具有两个非相邻非零的等概率随机二元矩阵的算法
Chr*_*ngs 18 random algorithm matrix
如果有人能指出我允许我的算法,那会很棒:
- 创建一个随机的方阵,条目为0和1,这样就可以了
- 每一行和每一列恰好包含两个非零项,
- 两个非零条目不能相邻,
- 所有可能的矩阵都是等概率的.
现在我设法实现以下第1点和第2点:这样的矩阵可以使用合适的行和列排列转换为具有块形式的对角块矩阵
1 1 0 0 ... 0
0 1 1 0 ... 0
0 0 1 1 ... 0
.............
1 0 0 0 ... 1
所以我使用[0,...,n-1]的分区从这样的矩阵开始,并通过随机排列行和列来加扰它.不幸的是,我找不到一种方法来整合邻接条件,我很确定我的算法不会平等对待所有矩阵.
更新
我已经设法达到了第3点.答案实际上就在我的鼻子底下:我正在创建的块矩阵包含考虑邻接条件所需的所有信息.首先是一些属性和定义:
- 一个合适的矩阵定义
[1, ..., n]了可以这样构建的排列:选择一行11.包含此条目的列恰好包含另一个条目,该条目在a不同于1的行上等于1.同样,行a包含列中的另一个条目1,其中包含行上的第二个条目1b,依此类推.这开始了一种排列1 -> a -> b ....
例如,使用以下矩阵,从标记的条目开始
v
1 0 1 0 0 0 | 1
0 1 0 0 0 1 | 2
1 0 0 1 0 0 | 3
0 0 1 0 1 0 | 4
0 0 0 1 0 1 | 5
0 1 0 0 1 0 | 6
------------+--
1 2 3 4 5 6 |
我们得到排列1 -> 3 -> 5 -> 2 -> 6 -> 4 -> 1.
- 这种置换的周期导致我前面提到的块矩阵.我还提到使用行和列上的任意排列来加扰块矩阵,以重建与要求兼容的矩阵.
但我正在使用任何排列,这导致一些相邻的非零条目.为了避免这种情况,我必须选择分隔块矩阵中相邻的行(和列)的排列.实际上,更准确地说,如果两行属于同一个块并且是循环连续的(块的第一行和最后一行也被认为是连续的),那么我想要应用的排列必须将这些行移动到非连续的最终矩阵的行(在这种情况下,我将调用两行不兼容).
所以问题就变成了:如何构建所有这些排列？
最简单的想法是通过随机添加与前一个行兼容的行来逐步构建置换.作为示例,考虑n = 6使用分区6 = 3 + 3和相应块矩阵的情况
1 1 0 0 0 0 | 1
0 1 1 0 0 0 | 2
1 0 1 0 0 0 | 3
0 0 0 1 1 0 | 4
0 0 0 0 1 1 | 5
0 0 0 1 0 1 | 6
------------+--
1 2 3 4 5 6 |
这里是行1,2并且3是相互不兼容的4,5和6.比如选择一个随机行3.
我们将编写一个置换作为数组:[2, 5, 6, 4, 3, 1]含义1 -> 2,2 -> 5,3 -> 6,...这意味着该行2的块矩阵将成为最终的矩阵的第一行,行5会不会成为第二行,依此类推.
现在让我们通过随机选择一行来构建一个合适的排列,比如说3:
p = [3, ...]
下一行随后将能与相容的剩余的行中随机选择的3:4,5和6.说我们选择4:
p = [3, 4, ...]
接下来的选择必须要中进行1和2,例如1:
p = [3, 4, 1, ...]
等等:p = [3, 4, 1, 5, 2, 6].
将此置换应用于块矩阵,我们得到:
1 0 1 0 0 0 | 3
0 0 0 1 1 0 | 4
1 1 0 0 0 0 | 1
0 0 0 0 1 1 | 5
0 1 1 0 0 0 | 2
0 0 0 1 0 1 | 6
------------+--
1 2 3 4 5 6 |
这样做,我们设法垂直隔离所有非零条目.必须对列进行相同的操作,例如通过使用排列p' = [6, 3, 5, 1, 4, 2]来最终获得
0 1 0 1 0 0 | 3
0 0 1 0 1 0 | 4
0 0 0 1 0 1 | 1
1 0 1 0 0 0 | 5
0 1 0 0 0 1 | 2
1 0 0 0 1 0 | 6
------------+--
6 3 5 1 4 2 |
所以这似乎非常有效,但构建这些排列需要谨慎,因为很容易被卡住:例如,使用n=6和分区6 = 2 + 2 + 2,遵循之前设置的构造规则可能会导致p = [1, 3, 2, 4, ...].不幸的是,5并且6不相容,所以选择一个或另一个使得最后的选择变得不可能.我想我发现所有情况都会导致死路一条.我会用r剩下的选择来表示:
p = [..., x, ?],r = {y}与x和y不兼容p = [..., x, ?, ?],r = {y, z}与y和z既是与不相容x(不可以作出选择)p = [..., ?, ?],r = {x, y}有x和y不相容(任何选择会导致情况1)p = [..., ?, ?, ?],r = {x, y, z}与x,y和z被循环连续(选择x或z将导致情况2,在选择y到情况3)p = [..., w, ?, ?, ?],r = {x, y, z}与xwy作为一个3个周期(均未x也不y可以被选择,在选择z将导致情况3)p = [..., ?, ?, ?, ?],r = {w, x, y, z}与wxyz作为一个4循环(任何选择将导致情况4)p = [..., ?, ?, ?, ?],r = {w, x, y, z}与xyz作为一个3个周期(选择w将导致情况4中,在选择任何其他会导致情况4)
现在似乎以下算法给出了所有合适的排列:
- 只要有超过5个数字可供选择,请在兼容的数字中随机选择.
- 如果剩下5个数字可供选择:如果剩余数字包含3个周期或4个周期,则中断该周期(即选择属于该周期的数字).
- 如果剩下4个数字可供选择:如果剩余数字包含三个循环连续数字,请选择其中一个.
- 如果剩下3个数字可供选择:如果剩余数字包含两个循环连续数字,请选择其中一个.
我很确定这允许我生成所有合适的排列,因此产生所有合适的矩阵.
不幸的是,每个矩阵都会被多次获得,具体取决于所选的分区.
sas*_*cha 10
介绍
这是一些原型方法,试图解决统一组合采样的更一般的任务 ,这对我们的方法来说意味着:我们可以将这种方法用于我们可以制定为SAT问题的一切.
它没有直接利用你的问题而且需要大量的绕道而行.绕过SAT问题可以帮助理论(更强大的一般理论结果)和效率(SAT求解器).
话虽如此,如果你想在几秒或更短的时间内(在我的实验中)进行采样,这不是一种方法,至少在关注均匀性时.
理论
该方法基于复杂性理论的结果,遵循以下工作:
GOMES,Carla P.; SABHARWAL,Ashish; 塞尔曼,巴特.使用XOR约束对组合空间进行近似均匀采样.在:神经信息处理系统的进展.2007. S. 481-488.
基本理念:
- 将问题表述为SAT问题
- 为问题添加随机生成的xors(仅对决策变量起作用!这在实践中很重要)
- 这将减少解决方案的数量(一些解决方案将变得不可能)
- 在一个循环(使用调整参数)中执行此操作,直到只剩下一个解决方案!
- SAT求解器或#SAT求解器(=模型计数)正在寻找一些解决方案
- 如果有多个解决方案:不会添加xors,但会完成重启:将random-xors添加到启动问题中!
保证:
- 在正确调整参数时,这种方法可以实现接近均匀的采样
- 这种调整成本很高,因为它基于近似可能的解决方案的数量
- 根据经验,这也可能是昂贵的!
- Ante的答案,提到数字序列A001499实际上给出了解决方案空间的一个很好的上限(因为它只是忽略了邻接约束!)
缺点:
- 对于大问题效率低下(一般情况下;不一定与MCMC和co等替代方案相比)
- 需要更改/减少参数以生成样本
- 那些减少的参数失去了理论保证
- 但凭经验:仍然可以取得好成绩!
参数:
在实践中,参数是:
- N:添加的xors数
- L:一个xor约束的最小变量数
- U:一个xor-constraint的最大变量数
N对于减少可能的解决方案的数量很重要.给定N常数,其他变量当然也会对此产生一些影响.
理论说(如果我正确解释),我们应该使用L = R = 0.5*#dec-vars.
这在实践中是不可能的,因为xor约束会伤害SAT解算器!
关于L和U的影响,这里有一些更科学的幻灯片.
它们称xor大小为8-20短XORS,而我们稍后需要使用更短的xors!
履行
最终版本
这是python中的一个非常hacky实现,使用此处的XorSample脚本.
正在使用的底层SAT解算器是Cryptominisat.
代码基本归结为:
- 将问题转换为联合正常形式
- 作为DIMACS-CNF
- 实施抽样方法:
- 调用XorSample(基于管道+基于文件)
- 调用SAT解算器(基于文件)
- 将样本添加到某个文件以供以后分析
代码:(我希望我已经警告过你关于代码质量的信息)
from itertools import count
from time import time
import subprocess
import numpy as np
import os
import shelve
import uuid
import pickle
from random import SystemRandom
cryptogen = SystemRandom()
""" Helper functions """
# K-ARY CONSTRAINT GENERATION
# ###########################
# SINZ, Carsten. Towards an optimal CNF encoding of boolean cardinality constraints.
# CP, 2005, 3709. Jg., S. 827-831.
def next_var_index(start):
next_var = start
while(True):
yield next_var
next_var += 1
class s_index():
def __init__(self, start_index):
self.firstEnvVar = start_index
def next(self,i,j,k):
return self.firstEnvVar + i*k +j
def gen_seq_circuit(k, input_indices, next_var_index_gen):
cnf_string = ''
s_index_gen = s_index(next_var_index_gen.next())
# write clauses of first partial sum (i.e. i=0)
cnf_string += (str(-input_indices[0]) + ' ' + str(s_index_gen.next(0,0,k)) + ' 0\n')
for i in range(1, k):
cnf_string += (str(-s_index_gen.next(0, i, k)) + ' 0\n')
# write clauses for general case (i.e. 0 < i < n-1)
for i in range(1, len(input_indices)-1):
cnf_string += (str(-input_indices[i]) + ' ' + str(s_index_gen.next(i, 0, k)) + ' 0\n')
cnf_string += (str(-s_index_gen.next(i-1, 0, k)) + ' ' + str(s_index_gen.next(i, 0, k)) + ' 0\n')
for u in range(1, k):
cnf_string += (str(-input_indices[i]) + ' ' + str(-s_index_gen.next(i-1, u-1, k)) + ' ' + str(s_index_gen.next(i, u, k)) + ' 0\n')
cnf_string += (str(-s_index_gen.next(i-1, u, k)) + ' ' + str(s_index_gen.next(i, u, k)) + ' 0\n')
cnf_string += (str(-input_indices[i]) + ' ' + str(-s_index_gen.next(i-1, k-1, k)) + ' 0\n')
# last clause for last variable
cnf_string += (str(-input_indices[-1]) + ' ' + str(-s_index_gen.next(len(input_indices)-2, k-1, k)) + ' 0\n')
return (cnf_string, (len(input_indices)-1)*k, 2*len(input_indices)*k + len(input_indices) - 3*k - 1)
# K=2 clause GENERATION
# #####################
def gen_at_most_2_constraints(vars, start_var):
constraint_string = ''
used_clauses = 0
used_vars = 0
index_gen = next_var_index(start_var)
circuit = gen_seq_circuit(2, vars, index_gen)
constraint_string += circuit[0]
used_clauses += circuit[2]
used_vars += circuit[1]
start_var += circuit[1]
return [constraint_string, used_clauses, used_vars, start_var]
def gen_at_least_2_constraints(vars, start_var):
k = len(vars) - 2
vars = [-var for var in vars]
constraint_string = ''
used_clauses = 0
used_vars = 0
index_gen = next_var_index(start_var)
circuit = gen_seq_circuit(k, vars, index_gen)
constraint_string += circuit[0]
used_clauses += circuit[2]
used_vars += circuit[1]
start_var += circuit[1]
return [constraint_string, used_clauses, used_vars, start_var]
# Adjacency conflicts
# ###################
def get_all_adjacency_conflicts_4_neighborhood(N, X):
conflicts = set()
for x in range(N):
for y in range(N):
if x < (N-1):
conflicts.add(((x,y),(x+1,y)))
if y < (N-1):
conflicts.add(((x,y),(x,y+1)))
cnf = '' # slow string appends
for (var_a, var_b) in conflicts:
var_a_ = X[var_a]
var_b_ = X[var_b]
cnf += '-' + var_a_ + ' ' + '-' + var_b_ + ' 0 \n'
return cnf, len(conflicts)
# Build SAT-CNF
#############
def build_cnf(N, verbose=False):
var_counter = count(1)
N_CLAUSES = 0
X = np.zeros((N, N), dtype=object)
for a in range(N):
for b in range(N):
X[a,b] = str(next(var_counter))
# Adjacency constraints
CNF, N_CLAUSES = get_all_adjacency_conflicts_4_neighborhood(N, X)
# k=2 constraints
NEXT_VAR = N*N+1
for row in range(N):
constraint_string, used_clauses, used_vars, NEXT_VAR = gen_at_most_2_constraints(X[row, :].astype(int).tolist(), NEXT_VAR)
N_CLAUSES += used_clauses
CNF += constraint_string
constraint_string, used_clauses, used_vars, NEXT_VAR = gen_at_least_2_constraints(X[row, :].astype(int).tolist(), NEXT_VAR)
N_CLAUSES += used_clauses
CNF += constraint_string
for col in range(N):
constraint_string, used_clauses, used_vars, NEXT_VAR = gen_at_most_2_constraints(X[:, col].astype(int).tolist(), NEXT_VAR)
N_CLAUSES += used_clauses
CNF += constraint_string
constraint_string, used_clauses, used_vars, NEXT_VAR = gen_at_least_2_constraints(X[:, col].astype(int).tolist(), NEXT_VAR)
N_CLAUSES += used_clauses
CNF += constraint_string
# build final cnf
CNF = 'p cnf ' + str(NEXT_VAR-1) + ' ' + str(N_CLAUSES) + '\n' + CNF
return X, CNF, NEXT_VAR-1
# External tools
# ##############
def get_random_xor_problem(CNF_IN_fp, N_DEC_VARS, N_ALL_VARS, s, min_l, max_l):
# .cnf not part of arg!
p = subprocess.Popen(['./gen-wff', CNF_IN_fp,
str(N_DEC_VARS), str(N_ALL_VARS),
str(s), str(min_l), str(max_l), 'xored'],
stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
result = p.communicate()
os.remove(CNF_IN_fp + '-str-xored.xor') # file not needed
return CNF_IN_fp + '-str-xored.cnf'
def solve(CNF_IN_fp, N_DEC_VARS):
seed = cryptogen.randint(0, 2147483647) # actually no reason to do it; but can't hurt either
p = subprocess.Popen(["./cryptominisat5", '-t', '4', '-r', str(seed), CNF_IN_fp], stdin=subprocess.PIPE, stdout=subprocess.PIPE)
result = p.communicate()[0]
sat_line = result.find('s SATISFIABLE')
if sat_line != -1:
# solution found!
vars = parse_solution(result)[:N_DEC_VARS]
# forbid solution (DeMorgan)
negated_vars = list(map(lambda x: x*(-1), vars))
with open(CNF_IN_fp, 'a') as f:
f.write( (str(negated_vars)[1:-1] + ' 0\n').replace(',', ''))
# assume solve is treating last constraint despite not changing header!
# solve again
seed = cryptogen.randint(0, 2147483647)
p = subprocess.Popen(["./cryptominisat5", '-t', '4', '-r', str(seed), CNF_IN_fp], stdin=subprocess.PIPE, stdout=subprocess.PIPE)
result = p.communicate()[0]
sat_line = result.find('s SATISFIABLE')
if sat_line != -1:
os.remove(CNF_IN_fp) # not needed anymore
return True, False, None
else:
return True, True, vars
else:
return False, False, None
def parse_solution(output):
# assumes there is one
vars = []
for line in output.split("\n"):
if line:
if line[0] == 'v':
line_vars = list(map(lambda x: int(x), line.split()[1:]))
vars.extend(line_vars)
return vars
# Core-algorithm
# ##############
def xorsample(X, CNF_IN_fp, N_DEC_VARS, N_VARS, s, min_l, max_l):
start_time = time()
while True:
# add s random XOR constraints to F
xored_cnf_fp = get_random_xor_problem(CNF_IN_fp, N_DEC_VARS, N_VARS, s, min_l, max_l)
state_lvl1, state_lvl2, var_sol = solve(xored_cnf_fp, N_DEC_VARS)
print('------------')
if state_lvl1 and state_lvl2:
print('FOUND')
d = shelve.open('N_15_70_4_6_TO_PLOT')
d[str(uuid.uuid4())] = (pickle.dumps(var_sol), time() - start_time)
d.close()
return True
else:
if state_lvl1:
print('sol not unique')
else:
print('no sol found')
print('------------')
""" Run """
N = 15
N_DEC_VARS = N*N
X, CNF, N_VARS = build_cnf(N)
with open('my_problem.cnf', 'w') as f:
f.write(CNF)
counter = 0
while True:
print('sample: ', counter)
xorsample(X, 'my_problem', N_DEC_VARS, N_VARS, 70, 4, 6)
counter += 1
输出看起来像(删除了一些警告):
------------
no sol found
------------
------------
no sol found
------------
------------
no sol found
------------
------------
sol not unique
------------
------------
FOUND
核心:CNF配方
我们为矩阵的每个单元引入一个变量.N = 20表示400个二进制变量.
Adjancency:
预先计算所有对称减少的冲突并添加冲突条款.
基础理论:
a -> !b
<->
!a v !b (propositional logic)
行/列基数:
这在CNF中难以表达,天真的方法需要指数数量的约束.
我们使用一些基于加法器电路的编码(SINZ,Carsten.对布尔基数约束的最佳CNF编码),它引入了新的辅助变量.
备注:
sum(var_set) <= k
<->
sum(negated(var_set)) >= len(var_set) - k
这些SAT编码可以放入精确的模型计数器中(对于小N;例如<9).解决方案的数量等于Ante的结果,这是正确转换的强烈迹象!
还有一些有趣的近似模型计数器(也非常基于xor约束),例如approxMC,它显示了我们可以用SAT配方做的另一件事.但实际上我还没能使用这些(approxMC = autoconf;没有评论).
其他实验
我还使用pblib构建了一个版本,为SAT-CNF配方使用更强大的基数配方.我没有尝试使用基于C++的API,而只使用简化的pbencoder,它自动选择一些最佳编码,这比我上面使用的编码更糟糕(这最好仍然是一个研究问题;通常甚至是冗余约束可以帮助).
实证分析
为了获得一些样本大小(考虑到我的耐心),我只计算N = 15的样本.在这种情况下,我们使用:
- N = 70 xors
- L,U = 4,6
我还用(100,3,6)计算了N = 20的一些样本,但这需要几分钟,我们减少了下限!
可视化
这里有一些动画(加强了我与matplotlib的爱恨交织):
编辑:与N = 5(NXOR,L,U = 4,10,30)的强力均匀采样(减少)比较:
(我还没有决定添加绘图代码.它和上面的一样丑陋,人们可能会对我的统计错误看起来太过分了;规范化和合作.)
理论
统计分析可能很难做到,因为潜在的问题是这种组合性质.甚至不完全清楚最终的单元格PDF应该是什么样子.在N =奇数的情况下,它可能是非均匀的并且看起来像棋盘(我做蛮力检查N = 5来观察这个).
我们可以肯定的一件事(imho):对称性!
给定一个cell-PDF矩阵,我们应该期望矩阵是对称的(A = AT).这在可视化中被检查,并且绘制了随时间的差异的欧几里德范数.
我们可以在其他一些观察中做同样的事情:观察到配对.
对于N = 3,我们可以观察以下对:
- 0,1
- 0,2
- 1,2
现在我们可以按行和每列执行此操作,并且应该期望对称性!
可悲的是,对于方差的说法可能并不容易,因此需要的样本来谈论信心!
意见
根据我的简化感知,虽然尚未实现收敛(或者我们远离均匀性),但是当前样本和单元格PDF看起来很好.
更重要的方面可能是两个规范,很好地逐渐减少到0.(是的;可以通过使用prob = 0.5进行转置来调整一些算法;但这不是在这里完成的,因为它会破坏它的目的).
潜在的后续步骤
- 调整参数
- 使用#SAT-solvers/Model-counters而不是SAT-solvers查看方法
- 尝试不同的CNF配方,特别是在基数编码和xor编码方面
- 默认情况下,XorSample使用类似tseitin的编码来以指数方式增长
- 对于较小的xors(使用时),使用naive编码(传播速度更快)可能是个好主意
- XorSample在理论上支持它; 但剧本在实践中的工作方式不同
- Cryptominisat以专用的XOR处理而闻名(因为它是用于分析加密包括许多xors的构建),并且可能通过天真编码获得一些东西(因为推断CNF爆炸的xors要困难得多)
- 对于较小的xors(使用时),使用naive编码(传播速度更快)可能是个好主意
- 默认情况下,XorSample使用类似tseitin的编码来以指数方式增长
- 更多的统计分析
- 摆脱XorSample脚本(shell + perl ......)
摘要
- 这种方法很一般
- 该代码产生可行的样本
- 应该不难证明,每个可行的解决方案都可以被采样
- 其他人已证明某些参数的均匀性的理论保证
- 不适合我们的参数
- 其他人根据经验/理论分析较小的参数(在这里使用)
- 哦不,你有照片!你本可以告诉我们你正在研究这个问题,然后我们就不会打扰了:-) (2认同)
(更新了测试结果,示例贯穿期和下面的代码片段.)
您可以使用动态编程来计算每个状态产生的解决方案的数量(以比蛮力算法更有效的方式),并使用这些(预先计算的)值来创建等概率随机解.
考虑一个7x7矩阵的例子; 在一开始,国家是:
0,0,0,0,0,0,0
意味着有七个相邻的未使用的列.在向第一行添加两个之后,状态可以是例如:
0,1,0,0,1,0,0
有两列现在有一列.在将第二行添加到第二行之后,状态可以是例如:
0,1,1,0,1,0,1
在填充三行之后,列可能最多有两行; 这有效地将矩阵分成两个独立的区域:
1,1,1,0,2,0,1 -> 1,1,1,0 + 0,1
这些区域是独立的,因为当将不相邻的规则添加到不同区域时,不相邻规则没有效果,并且区域的顺序对解决方案的数量没有影响.
为了将这些状态用作解决方案类型的签名,我们必须将它们转换为规范符号.首先,我们必须考虑这样一个事实:在它们中只有1个的列可能在下一行中不可用,因为它们在当前行中包含一个.因此,我们不得使用二进制表示法,而是使用三元符号,例如:
2,1,1,0 + 0,1
其中2表示此列在当前行中使用(而不是列中有2列).在下一步,我们应该将两个转换回一个.
此外,我们还可以镜像单独的组,将它们放入按字典顺序排列的最小符号:
2,1,1,0 + 0,1 -> 0,1,1,2 + 0,1
最后,我们将单独的组从小到大排序,然后按字典顺序排列,以便更大矩阵中的状态可以是例如:
0,0 + 0,1 + 0,0,2 + 0,1,0 + 0,1,0,1
然后,当计算每个状态产生的解的数量时,我们可以使用每个状态的标准符号作为关键字进行存储.
创建状态字典和每个状态的解决方案的数量只需要进行一次,而较大矩阵的表也可以用于较小的矩阵.
实际上,您将生成一个介于0和解决方案总数之间的随机数,然后对于每一行,您将查看可以从当前状态创建的不同状态,查看每个状态的唯一解决方案的数量生成,并查看哪个选项导致与您随机生成的数字对应的解决方案.
请注意,每个状态和相应的键只能出现在特定的行中,因此您可以将键存储在每行的单独词典中.
检测结果
使用未经优化的JavaScript的第一个测试给出了非常有希望的结果.使用动态编程,计算10x10矩阵的解决方案数量现在需要一秒钟,其中蛮力算法需要几个小时(这是算法的一部分,只需要执行一次).具有签名和解决方案数量的字典的大小随着每个步骤的大小逐渐减小2.5; 生成它的时间增长了大约3倍.
这些是创建的解决方案,状态,签名(字典的总大小)和每行(每行最大字典)的最大签名数:
size unique solutions states signatures max/row
4x4 2 9 6 2
5x5 16 73 26 8
6x6 722 514 107 40
7x7 33,988 2,870 411 152
8x8 2,215,764 13,485 1,411 596
9x9 179,431,924 56,375 4,510 1,983
10x10 17,849,077,140 218,038 13,453 5,672
11x11 2,138,979,146,276 801,266 38,314 14,491
12x12 304,243,884,374,412 2,847,885 104,764 35,803
13x13 50,702,643,217,809,908 9,901,431 278,561 96,414
14x14 9,789,567,606,147,948,364 33,911,578 723,306 238,359
15x15 2,168,538,331,223,656,364,084 114,897,838 1,845,861 548,409
16x16 546,386,962,452,256,865,969,596 ... 4,952,501 1,444,487
17x17 155,420,047,516,794,379,573,558,433 12,837,870 3,754,040
18x18 48,614,566,676,379,251,956,711,945,475 31,452,747 8,992,972
19x19 17,139,174,923,928,277,182,879,888,254,495 74,818,773 20,929,008
20x20 6,688,262,914,418,168,812,086,412,204,858,650 175,678,000 50,094,203
(使用简单的128位整数实现,使用C++获得的其他结果.)
例
5x5矩阵的字典如下所示:
row 0: 00000 -> 16 row 3: 101 -> 0
1112 -> 1
row 1: 20002 -> 2 1121 -> 1
00202 -> 4 1+01 -> 0
02002 -> 2 11+12 -> 2
02020 -> 2 1+121 -> 1
0+1+1 -> 0
row 2: 10212 -> 1 1+112 -> 1
12012 -> 1
12021 -> 2 row 4: 0 -> 0
12102 -> 1 11 -> 0
21012 -> 0 12 -> 0
02121 -> 3 1+1 -> 1
01212 -> 1 1+2 -> 0
解决方案总数为16; 如果我们随机选择一个从0到15的数字,例如13,我们可以找到相应的(即第14个)解决方案,如下所示:
state: 00000
options: 10100 10010 10001 01010 01001 00101
signature: 00202 02002 20002 02020 02002 00202
solutions: 4 2 2 2 2 4
这告诉我们第14个解决方案是选项00101的第二个解决方案.下一步是:
state: 00101
options: 10010 01010
signature: 12102 02121
solutions: 1 3
这告诉我们第二个解决方案是选项01010的第一个解决方案.下一步是:
state: 01111
options: 10100 10001 00101
signature: 11+12 1112 1+01
solutions: 2 1 0
这告诉我们第一个解决方案是选项10100的第一个解决方案.下一步是:
state: 11211
options: 01010 01001
signature: 1+1 1+1
solutions: 1 1
这告诉我们第一个解决方案是选项01010的第一个解决方案.最后一步是:
state: 12221
options: 10001
对应于随机选择的数字13的5x5矩阵是:
0 0 1 0 1
0 1 0 1 0
1 0 1 0 0
0 1 0 1 0
1 0 0 0 1
这是一个快速的代码示例; 运行代码片段以生成签名和解决方案计数字典,并生成一个随机的10x10矩阵(生成字典需要一秒钟;一旦完成,它会在半毫秒内生成随机解决方案):
function signature(state, prev) {
var zones = [], zone = [];
for (var i = 0; i < state.length; i++) {
if (state[i] == 2) {
if (zone.length) zones.push(mirror(zone));
zone = [];
}
else if (prev[i]) zone.push(3);
else zone.push(state[i]);
}
if (zone.length) zones.push(mirror(zone));
zones.sort(function(a,b) {return a.length - b.length || a - b;});
return zones.length ? zones.join("2") : "2";
function mirror(zone) {
var ltr = zone.join('');
zone.reverse();
var rtl = zone.join('');
return (ltr < rtl) ? ltr : rtl;
}
}
function memoize(n) {
var memo = [], empty = [];
for (var i = 0; i <= n; i++) memo[i] = [];
for (var i = 0; i < n; i++) empty[i] = 0;
memo[0][signature(empty, empty)] = next_row(empty, empty, 1);
return memo;
function next_row(state, prev, row) {
if (row > n) return 1;
var solutions = 0;
for (var i = 0; i < n - 2; i++) {
if (state[i] == 2 || prev[i] == 1) continue;
for (var j = i + 2; j < n; j++) {
if (state[j] == 2 || prev[j] == 1) continue;
var s = state.slice(), p = empty.slice();
++s[i]; ++s[j]; ++p[i]; ++p[j];
var sig = signature(s, p);
var sol = memo[row][sig];
if (sol == undefined)
memo[row][sig] = sol = next_row(s, p, row + 1);
solutions += sol;
}
}
return solutions;
}
}
function random_matrix(n, memo) {
var matrix = [], empty = [], state = [], prev = [];
for (var i = 0; i < n; i++) empty[i] = state[i] = prev[i] = 0;
var total = memo[0][signature(empty, empty)];
var pick = Math.floor(Math.random() * total);
document.write("solution " + pick.toLocaleString('en-US') +
" from a total of " + total.toLocaleString('en-US') + "<br>");
for (var row = 1; row <= n; row++) {
var options = find_options(state, prev);
for (var i in options) {
var state_copy = state.slice();
for (var j in state_copy) state_copy[j] += options[i][j];
var sig = signature(state_copy, options[i]);
var solutions = memo[row][sig];
if (pick < solutions) {
matrix.push(options[i].slice());
prev = options[i].slice();
state = state_copy.slice();
break;
}
else pick -= solutions;
}
}
return matrix;
function find_options(state, prev) {
var options = [];
for (var i = 0; i < n - 2; i++) {
if (state[i] == 2 || prev[i] == 1) continue;
for (var j = i + 2; j < n; j++) {
if (state[j] == 2 || prev[j] == 1) continue;
var option = empty.slice();
++option[i]; ++option[j];
options.push(option);
}
}
return options;
}
}
var size = 10;
var memo = memoize(size);
var matrix = random_matrix(size, memo);
for (var row in matrix) document.write(matrix[row] + "<br>");下面的代码片段显示了大小为10x10的矩阵的签名和解决方案计数字典.我使用了与上面的解释略有不同的签名格式:区域由'2'而不是加号分隔,而在前一行中有一个的列标有'3'而不是'2'.这显示了如何将键存储在文件中作为具有2×N位的整数(用2填充).
function signature(state, prev) {
var zones = [], zone = [];
for (var i = 0; i < state.length; i++) {
if (state[i] == 2) {
if (zone.length) zones.push(mirror(zone));
zone = [];
}
else if (prev[i]) zone.push(3);
else zone.push(state[i]);
}
if (zone.length) zones.push(mirror(zone));
zones.sort(function(a,b) {return a.length - b.length || a - b;});
return zones.length ? zones.join("2") : "2";
function mirror(zone) {
var ltr = zone.join('');
zone.reverse();
var rtl = zone.join('');
return (ltr < rtl) ? ltr : rtl;
}
}
function memoize(n) {
var memo = [], empty = [];
for (var i = 0; i <= n; i++) memo[i] = [];
for (var i = 0; i < n; i++) empty[i] = 0;
memo[0][signature(empty, empty)] = next_row(empty, empty, 1);
return memo;
function next_row(state, prev, row) {
if (row > n) return 1;
var solutions = 0;
for (var i = 0; i < n - 2; i++) {
if (state[i] == 2 || prev[i] == 1) continue;
for (var j = i + 2; j < n; j++) {
if (state[j] == 2 || prev[j] == 1) continue;
var s = state.slice(), p = empty.slice();
++s[i]; ++s[j]; ++p[i]; ++p[j];
var sig = signature(s, p);
var sol = memo[row][sig];
if (sol == undefined)
memo[row][sig] = sol = next_row(s, p, row + 1);
solutions += sol;
}
}
return solutions;
}
}
var memo = memoize(10);
for (var i in memo) {
document.write("row " + i + ":<br>");
for (var j in memo[i]) {
document.write(""" + j + "": " + memo[i][j] + "<br>");
}
}