23 Julia编程示例–递归趣例

23.1 汉诺塔问题

设有三根柱子A, B, C, 有大小依次为\(1, 2, \dots, n\)\(n\)个空心圆盘, 在A柱子上依次从低向上穿了\(n, n-1, \dots, 1\)大小的圆盘。 任务是要把这\(n\)个圆盘移动到C柱子上, 仍按照从低向上越来越小的次序。 移动的要求为:

  • 每次仅移动一个圆盘到另一个柱子上;
  • 每次的移动,都不能使得增加一个圆盘的柱子上的大圆盘压在小圆盘上。

这个问题是典型的递归问题:

  • 如果仅有一个圆盘(\(n=1\)), 就直接将A柱子的圆盘放到C柱子,问题解决;
  • 如果有\(n>1\)个圆盘,分为三个步骤:
    • 将A柱子上面的\(n-1\)个圆盘以C柱为缓存转移到B柱子;
    • 将A柱子最下面的\(n\)号圆盘转移到C柱子;
    • 将B柱子的\(n-1\)个圆盘以A柱子为缓存转移到C柱子,问题解决。

其中,在将B柱子的\(n-1\)个圆盘转移到\(C\)时, 会将最上面的\(n-2\)个先转移到A柱子, 所以递归过程中三个柱子的地位是可以互相改变的。

用三个数组表示A, B, C三个柱子, 数组第一个元素记录柱子编号, 数组第二个元素表示放在最下面的圆盘编号(1号最小,\(n\)号最大), 后续元素依次表示上面的圆盘。

function demo_hanoi(n=3)
    # 用第一个元素保存柱子编号,
    # 后续元素保存圆盘大小
    a = [1; collect(n:-1:1)]
    b = Int[2]
    c = Int[3]
    labels = ("A", "B", "C")

    function solve!(a, b, c, n)
        if n == 1
            plate = pop!(a)
            push!(c, plate)
            println("$(plate): $(labels[a[1]]) ==> $(labels[c[1]])")
        else 
            # 将A柱上面的n-1个通过C移动到B
            solve!(a, c, b, n-1)
            # 将A柱最下面的n号圆盘移动到C
            plate = pop!(a)
            push!(c, plate)
            println("$(plate): $(labels[a[1]]) ==> $(labels[c[1]])")
            # 将B柱上面的n-1个通过A移动到C
            solve!(b, a, c, n-1)
        end 
    end 

    solve!(a, b, c, n)
end 

demo_hanoi()

结果如:

1: A ==> C
2: A ==> B
1: C ==> B
3: A ==> C
1: B ==> A
2: B ==> C
1: A ==> C

23.2 高斯八皇后问题

国际象棋的棋盘是 8×8 的格子。 国际象棋手马克斯·贝瑟尔于 1848 年提出如下问题: 在 8×8 格的国际象棋上摆放八个皇后, 使其不能互相攻击, 即任意两个皇后都不能处于同一行、同一列或同一斜线上, 问有多少种摆法。 高斯认为有 76 种方案。 1854 年在柏林的象棋杂志上不同的作者发表了 40 种不同的解, 后来有人用图论的方法解出 92 种结果。

23.2.1 穷举法

要用编程的方法解决, 首先想到的是, 每个皇后有64个位置可选, 穷举可以考虑\(64^8 \approx 2.8 \times 10^{14}\)种组合, 即百万亿级别的组合, 计算量太大。

实际上, 对于每一种可行的摆法, 因为每一行有且仅有一个皇后, 穷举时所以只需要考虑每一行的皇后在那一列, 有\(8^8 = 16,777,216\)种组合, 即约一千万种组合。 这个组合数目虽然很大, 但对于现代计算机来说可以很快地完成穷举。

这里记录一个Python编程实现。 将一种可能的摆法存入一个长度为8的列表, 每个元素代表对应的行的皇后位置。 虽然可以写一个8重循环, 但是这样的程序过于繁琐。 所以, 采用八进制的想法, 从[0,0,0,0,0,0,0,0]出发, 每次给个位加1, 遇到8就向下一位进1, 直到循环\(8^8\)次。 写一个函数判断新的一个位置是否合法。

import time

## 检查一种摆法是否合法
## x是长度为8的列表,代表从低位到高位的8进制数字,
## 用每个数字代表一行的棋子所在的列。
def is_solution(x):
    for i in range(7):
        for j in range(i+1, 8):
            di = x[j] - x[i]
            if di==0 or di==j-i or di==i-j :
                return False
    return True
            
## 穷举进入下一个摆法
def next_board(x):
    x[0] += 1
    for i in range(7):
        if(x[i]==8):
            x[i+1] += 1
            x[i] = 0
    return x
            
def queens_exaust():
    board = [0]*8
    nmax = 8**8
    ## boards用来存放所有的摆法,
    ## 每种摆法表示为8个元素的列表
    boards = []
    nb = 0
    for i in range(nmax):
        if is_solution(board):
            print(str(nb) + ',' + ','.join(map(str, board)))
            boards.append(board.copy())
            nb += 1
        
        board = next_board(board)
    return boards


time_start=time.time()
res1 = queens_exaust()
time_end=time.time()
print('穷举法用时', time_end-time_start)
print('解的个数:', len(res1))

用时50秒,得到了92个解。

下面是Julia版本的程序:

## 检查一种摆法是否合法
## x是长度为8的数组,第i个元素是第i行的棋子所在的列序号
function is_solution(x)
    for i in 1:7
        for j in (i+1):8
            di = x[j] - x[i]
            if di==0 || di==j-i || di==i-j
                return false
            end
        end
    end
    return true
end

## 穷举进入下一个摆法
## 把长度为8的数组x看成是有8位数的8进制数,但数字都加1,
## 次序是从低位到高位
function next_board(x)
    x[1] += 1
    for i in 1:7
        if x[i]==9
            ## 进位
            x[i+1] += 1
            x[i] = 1
        end
    end

    return x
end

function queens_exaust()
    board = ones(Int8, 8)
    nmax = 8^8
    nbmax = 1000
    nb = 0
    boards = zeros(Int8, nbmax, 8)
    for i in 1:nmax
        if is_solution(board)
            nb += 1
            println(board)
            boards[nb,:] = board
        end

        board = next_board(board)
    end
    boards = boards[1:nb,:]
    return boards
end
@time queens_exaust()

用时0.93秒,得到92个解。

23.2.2 递归算法

穷举法完全不考虑已经摆放的棋子, 所以另一种更好的做法是使用递归算法, 先在第一行摆下一颗棋子, 这有8种摆法; 然后,在第二行试图摆下第二颗棋子, 使得与第一颗棋子不冲突; 再考虑第三颗棋子, 如此一直到第八颗棋子。 这个过程中如果某颗棋子无法摆放, 就退回到上一颗棋子的下一个位置; 找到一种摆法后, 也考虑最后一颗棋子的下一个位置。

下面是从百度百科网站复制的Python版本的递归算法程序。

## 用长度为8的列表A表示一种摆法,
## 用cur表示当前要摆放的棋子序号,0 <= cur < 8
def queens(A, cur=0):
    if cur == len(A):
        ## 如果当前棋子序号为8,表示找到了一种摆法
        print(','.join(map(str, A)))
        boards.append(A.copy())
        return 0
    ## 对每颗棋子col遍历
    for col in range(len(A)):
        ## 将当前棋子cur摆放在col列,
        A[cur] = col
        ## 设置摆放是否合法的初值为真
        flag = True
        ## 判断cur号棋子摆放在col列,是否与0到cur-1号棋子有冲突
        ## row是用来遍历0到cur-1号棋子的变量
        for row in range(cur):
            ## 以下判断摆放在第cur行、第col列的棋子,
            ## 是否与第row行、A[row]列的棋子冲突
            if A[row] == col or abs(col - A[row]) == cur - row:
                flag = False
                break
        ## 如果摆放第cur号棋子在第col列成功,
        ## 就考虑下一颗棋子(递归调用);
        ## 如果摆放不成功,就摆放到col的下一列继续判断
        ## 这样会遍历每颗棋子与前面棋子的不冲突的位置
        if flag:
            queens(A, cur+1)

boards = []
queens([None]*8)

也是92个结果,用时仅30毫秒。

下面是Julia版本的程序:

## 递归算法
## 用长度为8的数组A表示一种摆法,
## 用cur表示当前要摆放的棋子序号,1 <= cur <= 8
function queens()
    nbmax = 1000
    nb = 0
    boards = zeros(Int8, nbmax, 8)

    function queens_rec(A, cur=1)
        if cur == 9
            ## 如果当前棋子序号为9,表示找到了一种摆法
            nb += 1
            boards[nb,:] = A
            ##println(join(string.(A), ","))
            return
        end # if cur

        ## 对当前要摆放的棋子遍历所有列位置试摆放
        for col in 1:8
            ## 将当前棋子cur摆放在col列,
            A[cur] = col
            ## 设置摆放是否合法的初值为真
            flag = true
            ## 判断cur号棋子摆放在col列,是否与1到cur-1号棋子有冲突
            ## row是1到cur-1号棋子的变量
            for row in 1:(cur-1)
                ## 以下判断摆放在第cur行、第col列的棋子,
                ## 是否与第row行、A[row]列的棋子冲突
                if A[row] == col || abs(col - A[row]) == cur - row
                    flag = false
                    break
                end
            end # for row

            ## 如果摆放第cur号棋子在第col列成功,
            ## 就考虑下一颗棋子(递归调用);
            ## 如果摆放不成功,就摆放到col的下一列继续判断
            ## 这样会遍历每颗棋子与前面棋子的不冲突的位置
            if flag
                queens_rec(A, cur+1)
            end
        end # for col
    end # function queens_rec

    queens_rec(zeros(Int8, 8), 1)
    boards = boards[1:nb,:]
    println("找到解的个数:", nb)

    return boards
end # function queens
@time queens()

找到92个解,用时0.1秒。

23.2.3 找到变换等价的摆法

还可以继续考虑这些解法中哪些是不能通过左右反射、上下反射, 旋转90、180、270度、沿主对角线反转、沿反对角线反转得到的, 按如果能通过变换得到就分为一组。 可以验证,这些变换构成一个变换群, 从而可以用来定义等价类, 可以通过变换变成相同的解法为同一类。 找到12个类。

import DataFrames: DataFrame
import CSV

function board_equals(x, y)
    ## 棋盘的8元素存储到8x8的0-1矩阵存储
    bs = length(x)
    function board_to_mat(z)
        mat = zeros(Int,bs,bs)
        for i in 1:bs
            mat[i,z[i]] = 1
        end # for i
        return mat
    end # board_to_mat

    xm = board_to_mat(x)
    ym = board_to_mat(y)

    # 左右反射
    function reflect_lr(z)
        return z[:, end:-1:begin]
    end # reflect_lr

    # 上下反射
    function reflect_ud(z)
        return z[end:-1:begin, :]
    end # reflect_ud

    function idtran(z)
        return z
    end
    
    # 逆时针旋转90度, rotl90
    # 旋转180度, rot180
    # 顺时针旋转90度,rotr90
    # 转置:transpose
    # 沿反对角线反转:reflect_lr ∘ rotl90

    funcs = [reflect_lr, reflect_ud, 
        rotl90, rotr90, rot180,
        transpose, reflect_lr  rotl90]

    for i in 1:length(funcs)
        f = funcs[i]
        ynew = f(ym)
        if all(ynew .== xm)
            return true
        end # if
    end # for f

    return false
end # board_equals

## 标记通过变换可以得到的摆法
function group_boards(boards)
    b = copy(boards)
    nb = size(b, 1)    # 摆法个数
    g = zeros(Int, nb) # 分组编码,每种摆法一个值
    ng = 1   # 分组编码初值
    g[1] = 1 # 第一中摆法分组为1 

    for j in 2:nb # 对从第二种摆法开始的每一种摆法
        newgroup = true
        for i in 1:(j-1) 
            ## 判断第j种摆法是否与前面的第i种摆法可转换,
            ## 如果可转换,将第j摆法标记为第i中摆法的分组编码
            if board_equals(b[i,:], b[j,:])
                g[j] = g[i]
                newgroup = false
                break
            end # if
        end # for i
        if newgroup
            ng += 1
            g[j] = ng
        end # if
    end # for j

    ## 保存为csv
    println("找到分组个数:", ng)
    grouped = [1:nb boards g]
    qn = ["solution"; "Q" .* string.(1:8); "group"]
    df = DataFrame(grouped, Symbol.(qn))
    sort!(df, :group)
    df[:, :solution] = 1:nb
    CSV.write("queens_group_julia.csv", df)
    
    return grouped
end # group_boards

boards = queens()
grouped = group_boards(boards);

结果的CSV文件:

solution,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,group
1,1,5,8,6,3,7,2,4,1
2,1,7,5,8,2,4,6,3,1
3,3,6,4,2,8,5,7,1,1
4,4,2,7,3,6,8,5,1,1
5,5,7,2,6,3,1,4,8,1
6,6,3,5,7,1,4,2,8,1
7,8,2,4,1,7,5,3,6,1
8,8,4,1,3,6,2,7,5,1
9,1,6,8,3,7,4,2,5,2
10,1,7,4,6,8,2,5,3,2
11,3,5,2,8,6,4,7,1,2
12,4,7,5,2,6,1,3,8,2
13,5,2,4,7,3,8,6,1,2
14,6,4,7,1,3,5,2,8,2
15,8,2,5,3,1,7,4,6,2
16,8,3,1,6,2,5,7,4,2
17,2,4,6,8,3,1,7,5,3
18,3,8,4,7,1,6,2,5,3
19,4,2,8,6,1,3,5,7,3
20,4,7,3,8,2,5,1,6,3
21,5,2,6,1,7,4,8,3,3
22,5,7,1,3,8,6,4,2,3
23,6,1,5,2,8,3,7,4,3
24,7,5,3,1,6,8,2,4,3
25,2,5,7,1,3,8,6,4,4
26,3,6,2,7,1,4,8,5,4
27,4,1,5,8,2,7,3,6,4
28,4,6,8,3,1,7,5,2,4
29,5,3,1,6,8,2,4,7,4
30,5,8,4,1,7,2,6,3,4
31,6,3,7,2,8,5,1,4,4
32,7,4,2,8,6,1,3,5,4
33,2,5,7,4,1,8,6,3,5
34,3,6,2,7,5,1,8,4,5
35,3,6,8,1,4,7,5,2,5
36,4,8,1,5,7,2,6,3,5
37,5,1,8,4,2,7,3,6,5
38,6,3,1,8,5,2,4,7,5
39,6,3,7,2,4,8,1,5,5
40,7,4,2,5,8,1,3,6,5
41,2,6,1,7,4,8,3,5,6
42,3,1,7,5,8,2,4,6,6
43,3,5,7,1,4,2,8,6,6
44,4,6,1,5,2,8,3,7,6
45,5,3,8,4,7,1,6,2,6
46,6,4,2,8,5,7,1,3,6
47,6,8,2,4,1,7,5,3,6
48,7,3,8,2,5,1,6,4,6
49,2,6,8,3,1,4,7,5,7
50,3,7,2,8,6,4,1,5,7
51,4,2,5,8,6,1,3,7,7
52,4,8,5,3,1,7,2,6,7
53,5,1,4,6,8,2,7,3,7
54,5,7,4,1,3,8,6,2,7
55,6,2,7,1,3,5,8,4,7
56,7,3,1,6,8,5,2,4,7
57,2,7,3,6,8,5,1,4,8
58,2,8,6,1,3,5,7,4,8
59,4,1,5,8,6,3,7,2,8
60,4,7,5,3,1,6,8,2,8
61,5,2,4,6,8,3,1,7,8
62,5,8,4,1,3,6,2,7,8
63,7,1,3,8,6,4,2,5,8
64,7,2,6,3,1,4,8,5,8
65,2,7,5,8,1,4,6,3,9
66,3,6,4,1,8,5,7,2,9
67,4,2,7,3,6,8,1,5,9
68,4,8,1,3,6,2,7,5,9
69,5,1,8,6,3,7,2,4,9
70,5,7,2,6,3,1,8,4,9
71,6,3,5,8,1,4,2,7,9
72,7,2,4,1,8,5,3,6,9
73,3,5,2,8,1,7,4,6,10
74,4,6,8,2,7,1,3,5,10
75,5,3,1,7,2,8,6,4,10
76,6,4,7,1,8,2,5,3,10
77,3,5,8,4,1,7,2,6,11
78,3,6,8,2,4,1,7,5,11
79,3,7,2,8,5,1,4,6,11
80,4,2,8,5,7,1,3,6,11
81,5,7,1,4,2,8,6,3,11
82,6,2,7,1,4,8,5,3,11
83,6,3,1,7,5,8,2,4,11
84,6,4,1,5,8,2,7,3,11
85,3,6,2,5,8,1,7,4,12
86,3,6,8,1,5,7,2,4,12
87,4,2,7,5,1,8,6,3,12
88,4,7,1,8,5,2,6,3,12
89,5,2,8,1,4,7,3,6,12
90,5,7,2,4,8,1,3,6,12
91,6,3,1,8,4,2,7,5,12
92,6,3,7,4,1,8,2,5,12