论MC中的数学

塞巴斯特

嘛，可能我觉得6R人才众多又比较友好所以发这里吧。
一道困扰我很久的问题。
一道数学题，或者说，应用题。
问：在Minecraft里面9x9范围的平面内，最多能种多少格甘蔗

故事背景是这样的，大概我高一还是高二的时候，和朋友玩MC，里面装了各种MOD
“在很久很久以前，玩家A发现了镜中世界”
“于是决定把家搬进去、但是遇到了一个问题、镜中世界只有9X9X9的空间、种一般的田是可以的、但是甘蔗这种必须和水相邻的种进去就显得占地了。”
“于是就引伸出了甘蔗难题”
“至今未能解决”

我那时候自己随便试了下（后来上网查有没有别的方法的时候，好像大家的方法也都是这样），用了下面这个排列

这样的结果是，17格水，60格甘蔗，4格实在碰不到空着。
当时我们都以为这是最优解，甚至开玩笑说什么MC里的数学之类的

可能我觉得这件事有那么点意思吧，和我姐闲聊的时候说到了
然后我姐听了之后，问了下这个问题的详情，甘蔗要种有什么条件之类的
就抱着玩玩的心态，用几格水放在一起的组合方式，折腾出了61格

然后我就不淡定了。

而她老人家做出这个比我的最优解更好的最优解之后就没继续做下去了。

那时候挺震惊的，感受到了，数学的魅力？？？
反正一下子就对这问题特别感兴趣
现在借@px.凤翔九天 写的程序在试着穷举出最优解
不过我还是挺好奇，如果现实中遇到这样的题，有没有穷举以外的，类似我们平常做题的方法，能做出最优解。
数学发展到现在我感觉应该是有的。然而凭我的数学水平与搜索能力并不知道这种题属于什么范畴，更不用提自己搜到答案或搜到相关知识去学习着尝试自己解决。


那么，6R的各位对这道题有什么看法，思路，以及，有没有人能提供穷举以外的算出最优解的方法呢。
触手们大胆地来卖触吧（

colorlemon

大大们玩游戏，都在用数学考虑最优解问题。
我这种笨蛋玩游戏，只会考虑怎么种林子会比较好看，想了半天弄不出好看的样子还要去网上抄别人；然而第二天全死光像秃头一样（x

taroxd

设第 i 行有 xi 片水，则可以得到如下优化问题（不等式解释为，水能提供的甘蔗数量 >= 实际的甘蔗数量）

1. 
   
2. minimize x1 + x2 + ... + x9
   
3. subject to  
   
4.     2x1      + x2 >= 9 - x1
   
5.     2x2 + x1 + x3 >= 9 - x2
   
6.     ...
   
7.     2x9 + x8      >= 9 - x9
   
8.     (x1, x2, ..., x9) ∈ {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}^9
   

复制代码

解这个整数线性规划问题，得 (3, 1, 3, 1, 3, 1, 3, 1, 3) 是唯一最优解，达到最优值 19.

因此能种甘蔗的位置不超过 81-19=62 块.

对每一列也同理。因此，验证 62 块做不到，只需要考虑满足每行分别有  (3, 1, 3, 1, 3, 1, 3, 1, 3)  块水，每列分别有  (3, 1, 3, 1, 3, 1, 3, 1, 3)  块水的种法就可以了。我相信这是做不到的，但我懒得枚举或证明了。

因此答案是 61.

下图的例子由触瞎的 @RyanBem @RyanBern 以 GPL 许可提供：

感谢楼主提供有趣的题目。希望更多人能感受到数学的魅力，并报考北京大学数学系。

塞巴斯特

一开始眼瞎理解错了状元的意思以为是能62懒得举例，愣是画了2个小时
回来才发现原来是懒得证明62无法成立233
感谢解答_(:3
（可惜已经高考完了_(:3

（大概画了半小时左右发现如果要办到62得先解决个问题，偶数行（1格的那行）相交的那16个点


怎么说，得在满足行列313131313的前提下都种上甘蔗（即临近至少有1格水，如果临近一格放水都办不到的话那格本身也是放不了水的，所以直接假设种上了）
（像个大傻瓜一样试了一个半小时死活试不到解决方案也没法绕开，怕是理解错了回来看看状元怎么说，结果原来真的理解错了_(:3
（大概证明一下313131313的前提下没法让这16点旁边都有至少1格水（应该也有更简单的问题归纳）就能证明没办法了_(3
（再次感谢下帮忙想方案的天使姐姐和RB前辈

yangff

随便写写懒的优化了

1. #include <cstdio>
   
2. #include <algorithm>
   
3. #include <queue>
   
4. using namespace std;
   
5. // 0 water 1 veg
   
6. int dp[9][(1<<9)][1<<9];
   
7. 
   
8. bool checkok(int now, int upper, int upper1) {
   
9.         if (upper & 1){
   
10.                 if ((now & 1) && (upper1 & 1) && ((upper >> 1) & 1))
    
11.                         return false;
    
12.         }
    
13.         if ((upper >> 8) & 1){
    
14.                 if (((now>>8) & 1) && ((upper1>>8) & 1) && ((upper >> 7) & 1))
    
15.                         return false;
    
16.         }
    
17.         for (int i = 1; i < 8; i++) {
    
18.                 if ((upper >> i) & 1) {
    
19.                         if ( ((upper >> (i - 1))&1) && ((upper >> (i + 1))&1) && ((now >> (i))&1) && ((upper1 >> (i))&1)) {
    
20.                                 return false;
    
21.                         }
    
22.                 }
    
23.         }
    
24.         return true;
    
25. }
    
26. 
    
27. 
    
28. int count1(int x){
    
29.         int cnt = 0;
    
30.         while (x){
    
31.                 cnt+=x&1;
    
32.                 x >>= 1;
    
33.         }
    
34.         return cnt;
    
35. }
    
36. 
    
37. int main(){
    
38.         for (int i = 0; i < (1<<9); i++){
    
39.                 dp[0][i][(1 << 9) - 1] = count1(i);
    
40.         }
    
41. 
    
42.         for (int i = 1; i < 9; i++){
    
43.                 for (int now = 0; now < (1<<9); now++){
    
44.                         for (int upper = 0; upper < (1<<9); upper++)
    
45.                                 for (int upper1 = 0; upper1 < (1<<9); upper1++){
    
46.                                         if (checkok(now, upper, upper1)) {
    
47.                                                 dp[i][now][upper] = max(dp[i][now][upper], dp[i - 1][upper][upper1] + count1(now));
    
48.                                         }
    
49.                         }
    
50.                 }
    
51.         }
    
52. 
    
53.         int ans = 0;
    
54.         for (int upper = 0; upper < (1<<9); upper++) {
    
55.                 for (int upper1 = 0; upper1 < (1<<9); upper1 ++){
    
56.                         if (checkok((1<<9) - 1, upper, upper1)) {
    
57.                                 ans = max(ans, dp[8][upper][upper1]);
    
58.                         }
    
59.                 }
    
60.         }
    
61.         printf("%d\n", ans);
    
62. 
    
63. }
    

复制代码

RyanBern

状元说我的图是 GPL 的，因此我有义务把这张图怎么来的给写下来。

方法可能有点让楼主失望，仍然需要借助计算机完成，不过这里面的道理可不是一个简单的穷举法能说的清的。认识楼上 taroxd 和我的人可能知道，我们专业是一样的，因此脑回路可能非常像（我就不说状元的脑回路有 5% 是我调教出来的哼

作为一个数学系科学计算专业的学生，拿到这样的问题首先考虑的就是建立数学模型。楼上的状元已经给了一个很好的例子，但是那个建模非常简单，并且和原问题不等价。用数学语言来说，状元提供的解法只是原问题解的一个上界，从约束上来看是将约束集扩大来寻找一个必要条件。什么？说人话？那我这样讲：状元的做法只能告诉你就算你累秃了也不会找到一个比他所说的更好的解。我们对此自然是不太满意的。所以我这篇回贴实际就是告诉大家：上面问题的最优解就是 61，方案如图所示。

接下来就是我的计算过程。如果觉得没信心听懂的可以不用往下看了。拿到这个问题，我们的首要任务就是将它变成一个优化问题进行求解。对于优化问题，最重要的有两个要素：
- 目标函数：你的优化目标是什么-
- 自变量：目标函数所依赖的变量，以及对这些变量有无限制。

在我们这个问题里面，我们要优化的无非就是田地的最大块数，而变量就是我们在何处放置“水”。这个关系非常明确，一旦我们确定了所有“水”的位置，那我们就自动确定了该区域所能种植的田地数量。

下面就是用数学的语言描述这个问题。田地的最大快数非常好量化，但是在何处放置“水”这个变量本身不是数学对象。因此我们可以引入一个 81 个分量的向量 x 来表示各个位置上是否是水。如果是则对应位置为 1，如果不是则对应位置是 0。在这里我们将一个二维的表变成了一个一个长向量，看起来有些不直观。不过我们总有办法将这些对应的位置串起来，例如我们先将其拆成 9 行，然后首尾相接就可以做成一个长向量。

现在的问题来了，如何利用给定的一个向量来计算此方案下田地的块数？这就需要我们建立一个函数，将这个向量变成一个数字。这个函数应该有简单的结构。例如，我们现在考虑在 (i, j) 位置设置“水”，假设它在区域内部，那么它四周的块就可以安全地放置田地。从映射角度看来，(i, j) 位置对应着向量 x 的第 (9*i+j) 个分量（为了方便，i 和 j 从 0 开始计数），而这样放置，它四周的点 x(9*i+j-1)，x(9*i+j+1)，x(9*(i-1)+j)，x(9*(i+1)+j) 都会在这个点的影响内。如果用向量 x 表示某个地方是否是水，用向量 y 表示某个地方被水或者田地占用情况，那么我们可以构造这样的线性映射


如果直观点表示的话，A 的形状大概是这样，蓝点表示 1 其他地方表示 0


但是我们可能注意到了，映射 A 的作用就是表示某块地被占用的情况，大致可以分为两种
- 该片地被“水”占用，这从 A 的对角元可以看出。
- 该片地被“预计的田地”占用，这从 A 的非对角线可以看出。注意，如果一片田左右逢源，或者是和“水”重合了，那么这个位置要计算多次。

这样的 y 是不满足我们的要求的，我们想法是一片田地只算一次占用就好，并且去掉“水”的面积（因为显然水上不能种田）。这等价于我们将 y (=Ax) 中的所有比 1 大的元素都变成 1 就好。然后因为水的位置我们是知道的，变成 1 之后再减去 x，剩下的向量就是纯田地的信息了，最后我们把他们加起来，得到我们的目标函数值。

总之，我们写下我们的优化问题（这里 e 是一个全是 1 的向量，等价于对后面的分量求和）：


它等价于下面的问题，原因比较显然，这里就不证明了（方法是引入辅助变量 z）


这是一个整数规划问题，可以调用比较成熟的求解器进行求解。在这里我使用了商业软件 Gurobi 在 MATLAB 环境下进行求解，代码我贴在最后。得到的解大致就是这个样子了，蓝色是水，绿色是田。

MATLAB 代码复制

function [X, Out] = MC_gurobi(A, opts)
    %% Importing options
    init;
    %% Importing Gurobi API
    addpath(opts.path);
    %% Define model
    n = size(A, 1);
    model.obj = ones(2*n, 1);
    model.A = sparse([A, -eye(n)]);
    model.sense = '<';
    model.rhs = zeros(n, 1);
    model.lb = [zeros(n,1); -ones(n,1)];
    model.rb = [ones(n,1); ones(n,1)];
    model.vtype = 'I';
    model.modelname = 'MC-minimum';
    %% Set parameters
    param.method = -1;
    param.outputflag = 0;
    %% call solver
    tt = tic;
    Res = gurobi(model, param);
    tt = toc(tt);
    %% Exporting solutions
    x = Res.x;
    X = x(1:n);
    Out.tt = tt;
    Out.flag = Res.status;
    Out.fval = Res.objval;
    Out.itr = Res.baritercount;
    %% nested function
    function init
        if ~isfield(opts, 'path'); opts.path = '/opt/gurobi701/linux64/matlab'; end;
    end
end

function [X, Out] = MC_gurobi(A, opts)


    %% Importing options


    init;


    %% Importing Gurobi API


    addpath(opts.path);


    %% Define model


    n = size(A, 1);


    model.obj = ones(2*n, 1);


    model.A = sparse([A, -eye(n)]);


    model.sense = '<';


    model.rhs = zeros(n, 1);


    model.lb = [zeros(n,1); -ones(n,1)];


    model.rb = [ones(n,1); ones(n,1)];


    model.vtype = 'I';


    model.modelname = 'MC-minimum';


    %% Set parameters


    param.method = -1;


    param.outputflag = 0;


    %% call solver


    tt = tic;


    Res = gurobi(model, param);


    tt = toc(tt);


    %% Exporting solutions


    x = Res.x;


    X = x(1:n);


    Out.tt = tt;


    Out.flag = Res.status;


    Out.fval = Res.objval;


    Out.itr = Res.baritercount;


    %% nested function


    function init


        if ~isfield(opts, 'path'); opts.path = '/opt/gurobi701/linux64/matlab'; end;


    end


end

MATLAB 代码复制

dim = 9;
n = dim ^ 2;
A = gallery('poisson', dim);
for i = 1:n; A(i, i) = -1; end;
 
[x, Out] = MC_gurobi(A, []);
y = A * x;
X = reshape(x, 9, 9);
Y = reshape(y, 9, 9);
 
axis off;
%rectangle('Position', [1, 1, 9, 9], 'FaceColor', [0, .7, 0]);
for i=1:dim
    for j=1:dim
        pos = [i, j, 1, 1];
        if X(i, j) == 1
            rectangle('Position', pos, 'FaceColor', [0, 0, .7]);
        elseif Y(i, j) < 0
            rectangle('Position', pos, 'FaceColor', [0, .7, 0]);
        else
            rectangle('Position', pos);
        end
    end
end

dim = 9;


n = dim ^ 2;


A = gallery('poisson', dim);


for i = 1:n; A(i, i) = -1; end;


 


[x, Out] = MC_gurobi(A, []);


y = A * x;


X = reshape(x, 9, 9);


Y = reshape(y, 9, 9);


 


axis off;


%rectangle('Position', [1, 1, 9, 9], 'FaceColor', [0, .7, 0]);


for i=1:dim


    for j=1:dim


        pos = [i, j, 1, 1];


        if X(i, j) == 1


            rectangle('Position', pos, 'FaceColor', [0, 0, .7]);


        elseif Y(i, j) < 0


            rectangle('Position', pos, 'FaceColor', [0, .7, 0]);


        else


            rectangle('Position', pos);


        end


    end


end

tjjlb

【与题无关】
以前我闲的没事的时候出了一些MC的卷子哈哈哈哈哈
https://pan.baidu.com/s/1gfoy8dd
可以有更无聊的人尝试做做看
甚至出了一套略显正（zhong）式（er）试题
https://pan.baidu.com/s/1pKFCzWn