魔方还原算法(二) 科先巴的二阶段算法

科先巴的二阶段算法

本文来具体介绍一种具体的魔方还原算法——科先巴的二阶段算法，有一部分相关内容在前篇讲述，主要是方向定义那一块儿，没有看的建议先看一下：

二阶段，顾名思义，解决问题分为两步，先完成一个目标，再最终复原。对于二阶段算法有一个生动的比喻，复原魔方就像是一条小船要在汪洋大海上行驶到一个固定的目的地。二阶段算法就是先让小船行驶到一个固定的特殊水域，再驶向最终的目的地，这显然比直接寻找目的地要容易简单的多。

比喻归比喻，终归抽象，我们来具体看看。

一、总述

为方便叙述，也避免搞混淆，我们定义任意打乱的魔方称为随机状态或者初始状态，处于特殊水域的那些状态称为目标状态，目的地为还原状态。

初始状态可以看作是由 U，R，F，D，L，B 这 6 个基本转动复合而成的，由这 6 个转动生成的群记为 $G=⟨U,R,F,D,L,B⟩$。

而科先巴定义的目标状态是只由 U，D，L2，R2，F2，B2 这 6 种转动复合生成的，同样的记为 $G_1=⟨U,D,L2,R2,F2,B2⟩$。

Ⅰ 目标状态性质

还记得上篇文章计算状态数时是怎么定义方向的么？再来复习复习重新看看：

1. 棱边有 2 种状态，未翻转(flip)用 0 表示，翻转用 1 表示， 只有 F，B 层的转动能够改变方向。

2. 角块有 3 种状态，未扭转(twist)用 0 表示，顺时针扭转用 1 表示，逆时针扭转用 2 表示，只有 U，D 层转动不改变角块方向。

目标状态由 U，D，L2，R2，F2，B2 这 6 种转动生成，而这 6 种转动不会引起角块和棱块方向的变化。因为不论怎么转动，角块拥有的 U 层面或 D 层面始终处于 U 层或者 D 层，所以方向不变。也不会引起棱块方向变化，因为只有 F，B 的转动会引起棱块方向的变化，而像 F2 这样的连续转动两次，翻转再反转棱块相当于未翻转，所以棱块的方向也是不会变化的。这 6 种转动也不会将中层块带到 U 层或者 D 层，所以处于中层的块只能处于中层。

因此目标状态的性质总结如下：

1. 棱块角块方向不变，即各自的方向值总和都为 0。
2. 还原状态下的中层块处于中层

Ⅱ 搜索算法

科先巴的二阶段算法使用的搜索算法是 IDA* 算法，也就是迭代加深搜索算法，而且两个阶段使用的都是 IDA* 算法。

关于IDA*，还有个 A*算法大家可能不太熟悉，也没什么神秘的，简单说说。A* 就是具有启发函数的广搜。IDA*就是有搜索深度限制，有启发函数的深搜。点到为止，详细的解释请百度 CSDN。

一般的，广搜是能够找到最短路径的，深搜不能，IDA*虽然使用的是深搜，但是它有深度限制，也是能够找到最短路径。比如说我搜索深度定为 1 ，没搜到，再定为 2 重新搜索，依次下去，就能够得到最短路径。

所以两个阶段都是能够得到相应的最短路径，能够使用最少的步数达到每个阶段的目标，但合起来会是最优解吗？

答案是否定的，因为阶段二得到的最优解是相对于阶段一来说的，举个例子

两个阶段都相应最优就好比上述方法 1 得到的距离，但两个阶段合起来明显不是最短距离，两点之间线段最短应该是方法 4。

那怎么解决呢，如何得到更优甚至最优的解法？也如上图所示，我们加深第一阶段的搜索深度，虽然第一阶段的路径变长了，但是第二阶段可能找到更短路径，这样它们合起来边可能更优。这样一步步下去，是总能找到最优解的，也就是说科先巴的二阶段算法一次搜索可能得不到最优解，但是只要给它时间，是一定能够找到最优解的。但是对于科先巴算法并不是拿来寻找最优解的，目的是短时间内找到一组比较优的解。

对于科先巴的二阶段算法应该有了很直观的基本认识了吧，如果还不清楚，再针对魔友举一个很形象的例子。咱们一般使用的 CFOP 还原公式，F2L 之后就应该是 OLL，PLL 的对吧，如果我们正常的按照 CFOP 的公式来，这几个阶段都得到的是最优解，合起来的话也是比较优的。但是如果我们在 F2L 时多花几步，是可能跳 O 甚至跳 P 的，这就会节省很多后续转动步骤，得到一个更优的解法。

上述有关两点之间线段最短的例子只是比喻，实际情况稍微有点不一样，前面说过二阶段算法就是小船先驶向一片特殊水域，再驶向目的地，但是某些情况的最短路径可能根本就不经过那片水域，如果还是先经过特殊水域的话，自然也就得不到最优解。那这是不是就与上述所讲的冲突呢？也不然，还原状态其实也是目标状态之一是吧，所以如果解决阶段一之后发现已经还原了，阶段二所需要的步骤降为 0 ，那这就是最优解。就好比完成 F2L 之后惊奇的发现跳 O 又 条 P 已经还原了。

好了，上述就是对科先巴的二阶段算法的一个总述，应该是很容易理解的。下面是算法细节，挺多，感兴趣的朋友可以看看，对魔方没什么兴趣的了解到这一步大概就够了。

不过学计算机的，关于魔方的还原算法还是很值得一看的。不知道大家有没有做过八数码的题，解魔方与解八数码本质上是一样的，英文也都用 puzzle 这个词来表示，只不过魔方的状态数比八数码要大得多，而且魔方方方正正的，有许多特殊的性质可以利用来优化算法，下面来具体看看：

二、定义魔方

Ⅰ 块层次

定义魔方，首先就得对角块，棱块这些基本元素标号，可按如下方式这样标号：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-O3jSeDeq-1627530609241)(https://cdn.jsdelivr.net/gh/Rand312/Rand_Picture_01@main/Magic_Cube_Pic/块的标号.6knwqhq4v2w0.png)]

代码表示如下：

typedef enum {URF,UFL,ULB,UBR,DFR,DLF,DBL,DRB} Corner;
typedef enum {UR,UF,UL,UB,DR,DF,DL,DB,FR,FL,BL,BR} Edge;

对角边进行标号了之后就可以从块这个角度来定义一个魔方：

struct corner_o
{	Corner c;    //角块
	char o;      //方向
};
struct edge_o
{	Edge e;      //棱块，也就是边
	char o;      //棱块方向
};
struct CubieCube   //一个魔方状态
{	
	struct corner_o co[8];   //8个棱块
	struct edge_o eo[12];    //12个角块
};

怎么存储的魔方状态？举个例子 c o [ U R F ] = { U R F , 0 } co[URF] = \{URF,0\} co[URF]={URF,0} 表示 URF 这个位置存储的 URF 这个块，方向为 0。 c o [ U R F ] = { U F L , 1 } co[URF] = \{UFL,1\} co[URF]={UFL,1} 表示 URF 这个位置存储的 UFL 这个块，方向为0 。

棱块同理，而在篇一就说过，我们不用考虑中心块，中心块保持不动，当做参考系。如此就可以用上述的 CubieCube 结构体来表示一个魔方状态，定义一个魔方。

Ⅱ 坐标层次(索引)

这名字是从英文名中直译过来的，这一层次简单来说就是我们要对魔方状态进行编码/哈希，给不同状态分配一个唯一数来标识，主要用来作为变换表，剪枝表的索引。当然不可能真的对魔方不同的完整状态进行编码，这状态数太多了，不现实。而是对魔方的部分状态编码，一个魔方的完整状态有 8 个角块，12 个棱块以及它们的方向组成，而算法就是将它们分割，比如某些棱块的位置，角块位置，角块方向，棱块方向等等，对它们进行编码哈希，这样数量就会小上很多。

因此主要是对于位置，方向两个指标编码，而编码方法主要用到以下三种：

① 方向

方向分为角块方向和棱块的方向，两者类似，在这儿以角块方向举例，假如有如下所示的一个带方向的角块排列：

第一行表示位置，第二行表示该位置上存储的角块和其方向，使用下面的方法进行编码：

$$2\times 3^6+0\times 3^5+0\times 3^4+1\times 3^3+1\times 3^2+0\times 3^1+0\times 3^0$$

应该很直观也很好理解吧，因为当 7 个角块方向确定之后，最后一个角块的方向也就确定了，所以不用编码最后一个角块的方向。代码如下：

short int idx_co = 0;    //corner orientation index，角块方向编码
for (Corner c = URF; c < DRB; c++) 
    idx_co = 3 * idx_co + co[c].o;

② 排列

对于位置变化常使用排列数，期间我们可能会用到多种排列，比如 8 个角块的排列，6 条边的排列等等，这些都是魔方部分状态的表示，我们需要进行编码，有关排列的编码方法有个名字叫做康托展开，同样直接来看一个例子来说明康托展开这种编码方式，假如 8 个角块有如下排列：

第一行表示位置，第二行表示对应位置上放的角块，由标号知道咱们的块是有大小之分的：$$URF<UFL<ULB<UBR<DFR<DLF<DBL<DRB$$ 所以第三行就表示位于左边的块比自己大的块有几个，比如说 URF 块（看第二行，第一行表示位置，第二行才是该位置上存放的块）左边有 3 个块，都比自己大，所以第三行值为 3；再比如 UBR 块，左边有 7 个块，比自己大的有 4 个，所以值为 4 。DFR 块处于第一个位置，左边没有块了，所以编码时第一个块不用考虑进去。

具体的编码：

$$1\times 1!+1\times 2!+3\times 3!+0\times 4!+1\times 5!+1\times 6!+4\times 7!$$

这个式子与上面第三行一一对应起来应该很好理解就不解释了，代码描述如下：

int idx_cp = 0;   //corner position index 角块位置编码
for (int i = DRB; i > URF; i--){
    int s = 0;
	for (int j = i - 1; i >= DRF; j--)
        if(co[j] > co[i]) s++;
    idx_cp = (idx_cp + s) * i;
}

③ 组合

对于位置变化有时也会用到组合数，需要对组合数进行编码，在阶段一我们会用到关于中间层 4 个棱块的一个组合数编码，阶段一的目标状态要求本来在中间层的棱块要处在中间层，所以算阶段一的棱块位置状态我们只用考虑中间层的 4 个棱块，有 $12\times 11\times 10\times 9=11880$ 种状态，但是中间层的 4 个棱块之间并无位置要求，所以还要除以它们自身的位置排列 $4!=24$ ，总共 495 种状态，而这就是一个组合数 $C_{12}^4$ ，我们要对这 495 种状态编码：

第一行表示位置，第二行表示存储的块，只用考虑中间层 4 个块的位置也就是 $FRFLBLBR$ 4 个棱块， 因为中间层的 4 个棱块无位置要求，所以用 X 代替，表明它们是相同的。第三行表示 X 对应的组合数是多少，以第二个举例说明怎么编码：

$$C_{11}^4+C_{10}^3+C_9^2+C_8^1$$

上述编码方式应该很容易看懂吧，按顺序来就行了，来看码：

int idx_slice = 0, x = 0   //中间层 4 个棱块状态的编码
for(int i = BR; i >= UR; i--)   //从后往前
    if(edge[i] >= FR && edge[i] <= BR)   //判断是否是中间层 4 条棱边   
        idx_slice += c_nk(11 - j, x + 1) //c_nk为，n选k的组合数值
        x += 1     //x加1，x就是上图括号里面第二个值

直接用到了组合数的值，所以要提前写一个求组合数的函数，n 选 k，如果 $k>n$ 规定等于 0，写函数时判断以下即可。上述编码也可以对另外的 8 个组合数和来表示，感兴趣的可以自己试试。

上述是编码过程，魔方的坐标层次表示，也就是将魔方的部分状态转变成一个特定的数来表示，编码方式不唯一，合理即可。我们还应需要一个逆过程，从编码到实际的魔方状态，这样才能实现上述两种层次的定义相互转化，这里不再赘述了，基本上就是完全反着来一遍。

三、定义转动

Ⅰ 位置变化

先来看篇一提到的 { 0 , 1 , 2 , 3 } \{0,1,2,3\} {0,1,2,3} 的排列问题，假如有两个排列，在程序中用数组如下表示：

int a[4] = {0, 1, 2, 3};
int b[4] = {1, 2, 3, 0};
int c[4] = {};

数组这个数据结构有两要素，一个是下标表示位置，一个是元素，都是同等重要的，但往往我们可能会忽略掉下标的作用。

b 这个数组表示什么？当然可以表示一个排列，0号位存的 1，1号位存的2 ，2号位存的3，3号位存的0。

这只是常规认识，还表示什么？本身就可以表示一个变换，它表示 0 号位的元素被 1 号位元素替代，1 号位的元素被 2 号位元素替代，2 号位的元素被 3 号位元素替代，3 号位的元素被 0 号位元素替代。使用这种替代的角度来描述变换有很直观的等式关系，比如，0 号位的元素被 1 号位的元素替代，可以写成 $a[0]=a[1]$ 。

所以给定一个排列，它不仅是排列，还可以看作是一个变换，这是因为排列本身就存在次序，从数组的角度来就是有个隐含的下标来表示元素所处的位置，通过一个元素前后次序的变化就可以说明一个变换操作。

特别的如果元素个数也就是下标取值范围等于元素的取值范围，我们就可以很方便的在其上定义封闭的变换。比如说上述的排列元素个数为 4，所以下标的取值范围为 0~3，而元素本身的取值也是 0~3 这 4 个数。这种情况就很方便定义变化操作，比如在 a 这个排列上做 b变换，得到的结果存在排列 c 里面，进行如下操作：

for(int i = 0; i < 4; i++)
    c[i] = a[b[i]];

总结下来就是如果给定一个排列，如果要把它当作是一个变换的话，那么不仅下标表示的是位置，它的元素也表示的是位置。如 b[2] = 3就表示位置 2 的元素被位置 3 的元素替代。

理解了上面之后，魔方转动的定义也就迎刃而解了，并不需要什么新的数据结构，前一节定义的魔方状态就可以一个变换，只不过在不同需求时我们对其解读不同罢了。

Ⅱ 方向变化

一个位置存储的不仅有块本身，还有块的方向，块可以按照上述替代的角度来看，但方向不行。

举个例子，在 a 状态下做 F 操作，前面说过转动操作也可以看作一个状态，我们把 F 操作记为 b ，得到的结果记为 ab，以 URF 这个位置为例来说明。

F 转动使得 URF 槽的元素被 UFL 槽的元素替代，可以写成

$$ab.co[URF].c=a.co[UFL].c=a.co[F(URF)].c=a.co[b.co[URF]].c$$

但是方向的取值不能用替代来表示，也就是说 $ab.co[URF].o\ne a.co[UFL].o$ ，F 转动下，UFL 位置上的元素移动到 URF 位置上，并顺时针扭转，所以对于 F 转动，方向应有以下变化：

$$ab.co[URF].o=(a.co[UFL].o+1)mod3=(a.co[UFL].o+b.co[URF].o)mod3=(a.co[b.co[URF]].o+b.co[URF].o)mod3$$

方向这里有点绕，可以实际找个魔方转一转模拟一下，物理实体上应该是很好理解的，主要是要用程序语言来表述的确是有一点绕。

Ⅲ 基本转动

因此我们应该能很容易的写出基本转动的定义，来看个例子 F ：

const cubieCube basicMoveF = {{{UFL,1},   //URF槽的元素被UFL槽的元素替代，且顺时针扭转，方向加1
    {DLF,2},{ULB,0},{UBR,0},{URF,2},{DFR,1},{DBL,0},{DRB,0}},
	  {{UR,0},{FL,1},{UL,0},{UB,0},{DR,0},{FR,1},
	  {DL,0},{DB,0},{UF,1},{DF,1},{BL,0},{BR,0}}};

Ⅳ 转动函数

转动分为角块转动和棱块转动，所以代码如下：

CubieCube cubeMove(CubieCube cc, int m){ //cc 上应用转动m
	  CubieCube ccRet;
	  cornerMultiply(&cc, &basicMove[m], &ccRet); //basicMove是基本转动数组，元素如上小节的基本转动
	  edgeMultiply(&cc, &basicMove[m], &ccRet);
	  return ccRet;
};

这个函数就是角块转动和棱块转动函数的封装，以角块的转动函数为例说明：

void cornerMultiply(const CubieCube *a, const CubieCube *b, CubieCube *ab){
	  for (Corner crn = URF; crn <= DRB; crn++){
		    ab->co[crn].c = a->co[b->co[crn].c].c;       //角块
		    ab->co[crn].o = (a->co[b->co[crn].c].o + b->co[crn].o) % 3   //角块方向
        /*********************/
    }
}

与科先巴的源程序不太一样，原程序还处理了镜像等情况，这里先不管，后面再说。上面讲述的角度是将 a 看作魔方的一个状态，b 看作是一个变换，但同样也可以将 a 看作是一个变换，那么 ab 就相当于一个复合操作，一样的道理，只是理解角度不同而已。

Ⅴ 转动表

转动在程序里面会经常用，每次都去做一遍变换操作太费时间，我们要建立转动表，只要在第一次运行程序时建立好这个表，以后的转动操作直接去查表就行了。

转动表有很多，取决于使用了哪些坐标来表示魔方状态，以角块方向为例：

unsigned short int twistMove[Ntwist][Nmove];   //Ntwist=3^7-1, Nmove=18

这里就用上了魔方状态定义的坐标层次，把它当作索引。举个例子说明上面二维数组的意思，假如有 $twistMove[a][U]=b$，说明状态 a 下，进行 U 操作，得到状态 b。a, b 是两个数，分别唯一标识着魔方的角块方向状态。

这个二维数组初始化也很容易，直接看码：

void initTwistMoveTable(){    //初始化角块方向转动表
    CubieCube a, b;
    int i,j,m;
    for (int i = 0;i < NTWIST; i++){     //枚举方向状态
        CubieCube a = twist2cube(i);     //坐标层次->块层次
        for (int j = U; j <= B; j++){    //6种基本转动
            for(int k = 0; k < 3; k++){  //基本转动的3种方向，如U,U2,U'，分别对应90°，180°，270°
                cornerMultiply(&a, &basicCubeMove[j], &b);
                twistMoveTable[i][j * 3 + k] = get_twist(b);   //变化后的状态从块层次->坐标层次
            }
            cornerMultiply(&a, &basicCubeMove[j], &b); //再转一次，返回原状态
        }
    }
}

四、对称

魔方是一个高度对称的物体，我们可以利用对称性来减少时间空间上的开销。

对称的形式有 4 种基本情况，原本状态，只颜色变化，只位置变化，镜像。

它们都是等价的，如果我们能解决其中一个，就能解决其他的，只要在解法上加上相应的对称变化即可。

来看一个同位置不同颜色的例子，看看原状态经过怎样的变化后得到的一种等价状态：

假如对还原状态的魔方进行 A(一个转动序列) 操作得到如下状态：

我们要得到它的一个等价状态，进行如下操作，还原状态先绕 UD 轴逆时针转动 90°，再做 A 操作，最后绕 UD 轴顺时针转动 90°。

看到这大家有没有联想到什么，是不是很像矩阵的相似变换，有一些还原算法就是把魔方存成矩阵的形式，魔方变换就是相应的矩阵变换，道理都是相通的。

Ⅰ 对称变换

① 基本变换

上述的一些等价形式都是通过对称变换得来的，对于每个状态加上自己有 48 种变换形式，一般来说每个状态也就有 48 种等价形式，它们都是由 4 种“基本对称”组合而成：

1. S_U4：围绕 UD 这个轴转动 90°，4 种
2. S_F2：围绕 FD 这个轴转动 180°，2 种
3. S_URF3：围绕 URF 这个角块转动120°，3 种
4. S_LR2：关于 LR 层的镜像，2 种

合起来 $4\times 2\times 3\times 2=48$ 种

具体怎么转动的我就不具体画图表示了(唉主要是变换的 3D 制图太难了，尝试了一下，效果不好，见谅)，既然是变换操作，所以对称变换还是可以用状态的块层次来表示，具体来看 4 种基本对称变换的定义，从块的替代关系中应该就能理解具体是怎么变换的了。

const CubieCube basicSymCube[4] =
{
	{{{URF,1},{DFR,2},{DLF,1},{UFL,2},{UBR,2},{DRB,1},{DBL,2},{ULB,1}},
	{{UF,1},{FR,0},{DF,1},{FL,0},{UB,1},{BR,0},
	{DB,1},{BL,0},{UR,1},{DR,1},{DL,1},{UL,1}}},//S_URF3
	{{{DLF,0},{DFR,0},{DRB,0},{DBL,0},{UFL,0},{URF,0},{UBR,0},{ULB,0}},
	{{DL,0},{DF,0},{DR,0},{DB,0},{UL,0},{UF,0},
	{UR,0},{UB,0},{FL,0},{FR,0},{BR,0},{BL,0}}},//S_F2
	{{{UBR,0},{URF,0},{UFL,0},{ULB,0},{DRB,0},{DFR,0},{DLF,0},{DBL,0}},
	{{UB,0},{UR,0},{UF,0},{UL,0},{DB,0},{DR,0},
	{DF,0},{DL,0},{BR,1},{FR,1},{FL,1},{BL,1}}},//S_U4
	{{{UFL,3},{URF,3},{UBR,3},{ULB,3},{DLF,3},{DFR,3},{DRB,3},{DBL,3}},
	{{UL,0},{UF,0},{UR,0},{UB,0},{DL,0},{DF,0},
	{DR,0},{DB,0},{FL,0},{FR,0},{BR,0},{BL,0}}} //S_LR2
};

48 种对称变换操作初始化也很容易，4 个 for 循环，调用 multiply 函数进行复合操作，得到的结果填到 symCube[48] 数组里面去，这里就不赘述重复操作了。

注意上面的镜像变换，定义角块方向变化全部加 3 ，也就是说镜像的角块方向取值为 $[3,5]$ 。

② 逆变换

上述为了得到等价状态做了两次对称变换，两种变换互为逆过程，群论里面叫做互为逆元，我们这需要为每一个对称变换找到它们对应的逆变换。

那怎么得到相应的逆变换？在群论里面互逆的两个元素复合操作等于幺元，也就是相当于乘法里面的互为倒数的两个数相乘为 1，而在魔方里面还原状态就是幺元就是那个 1 。所以我们只要对 48 个变换进行两两复合操作，看结果是否达到还原状态即可。

init inv_idx[48] -1;  //初始化为-1，表示每种对称变换的逆变换
for(i = 0; i < 48; i++){
    for(j = 0; j < 48; j++){
        CubieCube cc = multiply(symCube[i], symCube[j]);  //两种变换相乘
        if (cc == identity)    //如果结果等于幺元还原状态
            inv_idx[i] = j;    //填数组，对称变换i的逆变换是j
    }
}

Ⅱ 等价

对于上述的 48 种对称变换操作我们记为 $S[48]$ ，也就是数组 symCube[48]，假如有 A，B 两个状态，如果存在 $A=S[i]BS[i{]}^{-1}$ ,也就是 $S[i{]}^{-1}AS[i]=B$ ，我们则认为 A，B 是等价的。这种写法是不是就跟矩阵的相似变换一模一样了，本质上都是相通的，所以真的要学好数学啊。

既然是等价的，那么我们就可以把这些等价的状态集合在一起，也就是说我们对于状态进行了类似下面的分类：

KaTeX parse error: Unknown column alignment: 1 at position 15: \begin{array}{1̲} 0&A0&A1&A2&..…

因此对于每种状态我们只需要记录属于哪一个集合(行号)，所使用的对称变换(列号) ，这么一个集合我们称为等价类。如此一来就可以省下很多空间和时间，因为存的信息少了，搜索的状态空间也少了。在科先巴的二阶段算法里，保持了 UD 轴不动，所以只用到了 16 种对称变换。

加入一个状态属于等价类 y ，使用的对称变换为 i，那么我么就可以用新坐标 $y\times 16+i$ 来描述原状态。但是并不是每个等价类里面都有不同的 16 种状态，也就是说不是每个状态都有 16 种不同对称状态。比如说有等价类 y 和对称变换 i，j，$y\times 16+i$ 与 $y\times 16+j$ 是有可能描述的同一个状态的。

关键问题在于我们如何得知一个坐标它属于哪一个等价类，使用的是哪一个对称变换？直接来看伪码分析：

init idx2class[Nindex] -1   //初始化为-1，idx属于哪一个等价类
init idx2sym[Nindex] -1     //初始化为-1，idx使用的哪一种对称变换
init class2rep[NCLASS] -1   //初始化为-1，该等价类的代表
    
classidx = 0;               //等价类标号
CubieCube cc = get_cube();  //获取魔方初始状态的块层次表示

for(idx = 0; idx < Nindex; i++){    //魔方部分状态的坐标层次表示，比如说角块方向Nindex = Ntwist = 2187
    if(idx2class[idx] == -1){       //还没填充
		idx2class[idx] = classidx;   //填
        idx2sym[idx] = 0;            //等价类里面的第一个元素，使用的对称变换为第一种0
        class2rep[classidx] = idx;   //等价类里面的代表，idx最小的那个作为代表
    }
    else 
        continue;     //已填，跳过        
    
    /**上面是处理等价类中的第一个元素，这里直接应用16种对称变换来处理其他元素**/
    for(sym = 0; sym < 16; sym++){  
        idx_new = idxSymTable[idx][sym];    //对称变换表，在idx上应用对称变换sym得到新坐标idx_new
        if(idx2class[idx_new] == -1){
            idx2class[idx_new] = classidx;
            idx2sym[idx_new] = sym;
        }
    }
    classidx++;   //填完当前等价类，加一准备填下一个
}

Ⅲ 对称变换表

上述的代码里面已经使用过了，类似前面的转动表，对称变换，转动都是变换操作，转动有转动表，对称变换也得有相应的变换表。比如定义以下的表：

int moveSymTable[NMOVE][16];       //对每一种转动m做 S[i]*m*S[i]^-1
int twistSymTable[NTWIST][16];     //对角块方向twist做 S[i]*twist*S[i]^-1

具体的初始化操作同转动表，都是建立在状态的复合操作之上，以角块方向的对称变换表举例说明：

void initTwistConjugate(){
    CubieCube cc,ccSym,ccTmp;
    int i,j;
    for (int i = 0;i < NTWIST; i++){
        cc = twist2cube(i);     //坐标层次转块边层次
        for (int j = 0;j < 16; j++){
            cornerMultiply(&symCube[i], &cc, &ccTmp);    //S[i]*cc
			cornerMultiply(&ccTmp, &symCube[inv_idx[i]], &ccSym);	//S[i]*cc*S[i]^-1
            twistSymTable[i][j] = get_twist(ccSym);   //结果转为坐标层次
        }	
    }	
}

源程序里面命名用的是 conjugate 这个词，共轭，我这儿直接用的 Sym 来说明，私以为可能更直观一些。

五、剪枝

剪枝，应该是科先巴的二阶段算法的核心了，综合运用了前面定义的所有数据结构来构建一个新的查找表——剪枝表，表里面存放的信息是距离目标状态(阶段一)或者还原状态(阶段二)最少还需要多少步。如果这个步数大于了还能走的步数，剪掉回溯。

两个阶段的剪枝表可能有多张，每个阶段都要选取最大的那个数作为到目标/还原状态的最少需要步数，比如说表一告诉我当前状态至少还需要 5 步到达目标状态，表二告诉我至少还需要 6 步，则将 6 作为至少需要步数。还是很好理解的是吧，表一表二分别描述了魔方的一部分状态，每部分状态都达到目标状态才算达到目标状态，所以要选取最大的。这也是 IDA 算法里面的 $h(n)$ 这个函数*。

Ⅰ 建立剪枝表

剪枝表竟然能够出到目标状态的最少步数？是不是很神奇，其实没什么神奇之处，就是一个暴力广搜得到的结果。设定目标状态的深度为 0 ，然后对目标状态进行 18 种转动操作，得到第一层的状态结点深度为 1，再对第一层的每种状态转动，得到第二层状态结点，深度为2…具体看下面伪码：

init pruneTable1[Nindex] -1  //初始化所有元素为 -1
depth = 0                    //初始深度为 0
pruneTable[0] = 0;           //本就是目标状态，所以深度/所需步数为0
done = 1;                    //已填了一个初始状态
while done < Nindex:         //表还没填完
    for i = 0 to Nindex - 1
        if pruneTable[i] == depth   
            for each m in Nmove    //对每个状态进行18种转动操作
                index = moveTable[i][m]  //根据转动表得出转动后的状态坐标值
                if pruneTable[index] != -1    //如果还没填表，填
                    pruneTable[index] = depth + 1  //深度加 1，表下一层
                    done++          
    depth++    //填完一层结点，深度加 1 
    

本质上就是广搜，只是形式上与我们平常看到的框架可能不太一样，这种形式费点时但省空间不用像常见广搜框架那样储存节点，魔方的状态数很多，我么要采用更省空间的形式。

Ⅱ 存储技巧

如果知道初始状态距离目标状态有多远，而后每执行一步，距离目标状态要么要么加 1，要么不变，要么减 1，就这三种情况，那么我们就能够计算出转动后的状态距离目标状态多远。所以其实剪枝表里面我们只需要 2个 bit 来存储 实际距离 mod 3 的值。根据当前状态的距离，再根据剪枝表里面查到的 mod3 值就能够算出转动后的状态到目标状态的距离的。如此就大大减少了空间的使用。

但关键在于初始状态到目标状态距离有多远？你可能会说，查表嘛，查表没错，但是要记住，现在剪枝表里面存放的不是实际距离了，存放的是只用了 2 bits 的 mod3 距离。那怎么办呢？我们直接看伪码分析解决办法：

int get_depth_phase1(CubieCube cc){
    index = get_index(cc);      //得到初始状态的坐标值
    depth_mod3 = pruneTable[index];  //初始状态到目标状态的mod3距离
    depth = 0;                  //初始化0
    while index != SOLVED:      //若没到目标状态，循环
        if(depth_mod3 == 0) depth_mod3 = 3;   //mod3为0，表3倍数，设为3
        for each m in Nmove             //对于每个转动
            index_new = moveTable[index][m];
        	if (pruneTable[index1] == depth_mod3 - 1) //如果转动后距离更近
                depth++;          //实际距离加 1 
    			index = index1;   //更新坐标索引值
                depth_mod3--;     //mod3距离减 1
    			break;            //找到一个举例目标转动更近的转动，跳出循环
    return depth;
}

根据伪码应该能看出大概意思吧，因为我们从剪枝表里面查出来的是 mod3 距离，所以 mod3 距离从 3 到 0 一直循环的递减，实际距离从 0 开始一直递增，直到到达目标状态。因为我们每次更新两个距离值都是要转动后距离目标状态更近，所以最终得出的实际距离 depth 一定是初始状态到目标状态需要的最少步骤。这似乎都已经把魔方解了一遍了，但其实并不是，因为 index 只是魔方部分状态的坐标索引，所以这个函数算出来的值只是当作一个初始指标。

六、阶段一

Ⅰ 所用坐标

阶段一用了三个坐标：

1. twist，角块方向，取值范围 $[0,2186]$
2. flip，棱边方向，取值范围 $[0,2047]$
3. slice，中间层的 4 条棱边的选取状态，取值范围 $[0,494]$

前两个不说了，很熟了，第三个在坐标层次那一块儿也说过了。想想第一阶段目标状态的性质，前两个坐标就是来解决方向问题的。第一阶段的目标状态对块边的位置要求不高，只要求属于中层的块要在中层呆着，而且都不要求要在家里呆着，只要在中层就好，第三个坐标就是来解决位置问题。

Ⅱ 坐标优化

但其实科先巴做了优化，并未实际定义 slice 而是 slice_sorted 这个坐标，这个坐标就是 slice 的加强版，slice 不考虑中间层 4 条棱边的位置吗，假如要考虑，那就是 slice_sorted，取值范围 $[0,11879]$ 。但当然实际用的时候用的还是用的 slice，slice 才是我们阶段一的指标，也就是说把 slice_sorted 除以 24 再用。为啥这么麻烦，不直接定义 slice ？因为阶段二也要用到这个 slice_sorted 这个坐标，为了方便统一所以这样处理吧，但当然，只要你自己清楚不混淆，不管使用那种，甚至使用其他坐标都行。

这个坐标既不是组合数又不是排列数，怎么编码？虽然这个坐标两不是却两相像，$12\times 11\times 10\times 9=\frac{A_{12}^{12}}{C_{12}^4}$ ，所以这个坐标也可以由排列组合的编码方式来表示。

Ⅲ 坐标合并/索引计算

flip 和 slice 组合成了一个坐标，flipslice，$flipslice=2048\times slice+flip$ ，取值范围 $[0,(2048\times 495-1)]$ 。

但是由于对称性质，flipslice 可以由 classidx 和 sym 替代，取值范围分别为 $[0,64429]$ ，$[0,15]$。 $64430\times 16>2048\times 495$ ，这儿就体现了不是每个等价类里面都有 16 种不同的对称状态存在。

classidx, sym 再和 twist 组合成一个坐标，暂且叫做 sym_flipslice_twist，取值范围 $[0,64430\times 2187-1]$ 这个坐标就是阶段一剪枝表的索引。

$$sym\_flipslice\_twist=2187\times classidx+twistSymTable[twist][sym]$$
为什么要这么计算？没有为什么，科先巴这么算的，也可以不这么算，它就是一个索引值，哈希值，只要能够一一区分开就行。

Ⅳ 搜索算法

前面已经说过搜索算法用的是 IDA*，不多说，直接看伪码：

main:       //主线程
	CubieCube cc = get_cube();   //得到输入的魔方
	dist = get_depth_phase1(cc); //获取阶段1需要的最少步数
	twist = get_twist(cc);       //获取坐标twist
	flip = get_flip(cc);         //获取坐标flip
	slice_sorted = get_slice_sorted(cc);    //获取坐标slice_sorted
	for(togo1 = dist; togo1 <= 20; i++)    //IDA* 限制搜索深度体现点，深度不能超过togo1
        search1(twist, flip, slice_sorted, dist, togo1); //阶段一搜索

void search1(twist, flip, slice_sorted, dist, togo1){  //深搜
	if(togo1 == 0)   //阶段一已解决
    /**省略后面讲,大概就是满足一定条件后search2()**/
    else{
        for each permitted m in MOVE{    //对每个允许的转动
        	flip_new = flipMoveTable[flip][m];     //获取新flip
            twist_new = twistMoveTable[twist][m];  //获取新twist
            slice_sorted_new = sliceSortedMoveTable[slice_sorted][m];  //获取新slice_sorted
            
            flipslice = 2048 * slice_sorted_new / 24 + flip_new;  //计算flipslice，slice_sorted除以24再用
            classidx = flipslice2class[flipslice]; //flipslice属于哪个等价类
            sym = flipslice2sym[flipslice];   //记录使用的哪种对称变换
            
            sym_flipslice_twist = 2187 * classidx + twistSymTable[twist][sym];  //计算索引
            
            dist_new_mod3 = get_flipslice_twist_depth3(sym_flipslice_twist);  //查阶段一剪枝表，获取mod3距离
            dist_new = get_real_dist(dist_new_mod3);   //mod3距离转实际距离
            
            if (dist_new > togo1)   //最少需要步数大于允许步数，剪枝
                continue;
            
            maneuver1.append(m);   //做选择，m加进解法的转动序列
            search1(twist_new, flip_new, slice_new, dist_new, togo1 - 1);  //下一层搜索
            maneuver1.pop();       //撤销选择，尝试兄弟结点或回溯
        }
    } 	
}

七、阶段二

Ⅰ 所用坐标

1. corners，角块排列，8个，取值范围 $[0,40319]$
2. ud_edges，U层D层的棱边排列，8个，取值范围 $[0,40319]$
3. slice_sorted，中层的棱边排列，4个，取值范围 $[0,23]$

现在已经处于阶段一的目标状态，所有棱块棱边的方向都是好的，我们只需要解决它们的位置问题，所以这 3 个坐标都是排列数，按照排列数编码建立坐标即可。

Ⅱ 坐标合并/计算索引

coners，取值范围 $[0,40319]$ ，由于对称性，可以由 classidx 和 sym 替代，取值范围分别为 $[0,2767]$，$[0,15]$ ，2768 个等价类，同样的因为不是每个等价类里面都有 16 种不同的状态，所以 $2768\times 16>40320$

classidx, sym, ud_edges 组合成一个坐标，暂且称为 sym_corners_ud_edges，取值范围 $[0,40320\times 2768-1]$这个坐标就是阶段二剪枝表的索引。
$$sym\_corners\_ud\_edges=40320\times classidx+ud\_edgesSymTable[ud\_edges][sym]$$

corners，slice_sorted 组合成另一个坐标，暂且称为 corners_slice_sorted，取值范围 $[0,40320\times 24-1]$ ，这个坐标是阶段二另一个剪枝表的索引。
$$corners\_slice\_sorted=24\times corners+slice\_sorted$$

Ⅲ 搜索算法

void search2(corners, ud_edges, slice_sorted, dist, togo2){  //深搜框架
    if(togo2 == 0)     //阶段二已解决
        maneuver = maneuver1 + maneuver2;   //合并两阶段的解法
    	if(len(maneuver) < shortest)
            shortest = len(maneuver)      //记录最短解法的长度
    else{
        for each permitted m in MOVE{     //对每个允许的转动
            corners_new = cornersMoveTable[corners][m];     //获取转动后的corners
            ud_edges_new = ud_edgesMoveTable[ud_edges][m];  //获取转动后的ud_edges
            slice_sorted_new = sliceSortedMoveTable[slice_sorted][m]; //获取转动后的slice_sorted
            
            classidx = corner2classidx[corners_new];    //获取corners的代表
            sym = corner2sym[corners_new];              //获取使用的对称变换
            
            sym_corners_ud_edges = 40320 * classidx + ud_edgesSymTable[ud_edges][sym]; //计算索引
            corners_slice_sorted = 24 * corners + slice_sorted;   //计算索引
            
            dist_new_mod3 = get_corners_ud_edges_depth3(sym_flipslice_twist);  //查剪枝表，获取转动后的mod3距离
            dist_new = get_real_dist(dist_new_mod3);    //mod3距离转成实际距离
            another_dist = corners_slice_sorted_pruneTable[corners_slice_sorted]  //查另一个剪枝表，查表，获取距离
                
            if(max(dist_new, another_dist) >= togo2) //两个距离指标，如果最大的那个大于允许的步数，剪枝
                continue;
            
            maneuver.append(m);      //做选择，m加进阶段二的解法
            seach2(corners_new, ud_edges_new, slice2_new, dist_new, togo2 - 1);  //搜索下层结点
            maneuver.pop();          //撤销选择，尝试兄弟结点或回溯
        }
    }
}

八、完善补充

科先巴的二阶段算法已经从头至尾的梳理了一遍，现在来看看其中省略或简单略过的部分

Ⅰ 完善阶段一

void search1(twist, flip, slice, dist, togo1){  //深搜
	if(togo1 == 0){   //阶段一已解决
        /***计算阶段二需要的参数***/
        corners = get_corners(cc);    //获取初始状态的corners坐标
        for m in maneuver1            //根据阶段一的解法，转动更新corners
            corners = cornersMoveTable[corners][m];
        
        //计算阶段二最多允许的步数，目前的最优解法的长度减去当前阶段一已花的步数，且阶段二用的步数不能超过11-1=10步
        //阶段二的步数下降的很快，通常不超过9步
        togo2_limit = min(shortest - len(maneuver1), 10);   
        corners_slice_sorted = 24 * corners + slice_sorted;   //计算索引
        if(corner_slice_sorted_pruneTable[corners_slice_sorted] >= togo2_limit) //查表剪枝
            return
        
        ud_edges = get_ud_edges(cc);  //获取初始状态的ud_edges坐标
        for m in maneuver1            //根据阶段一的解法，转动更新ud_edges
            ud_edges = ud_edgesMoveTable[ud_edges][m];
        
        dist2 = get_depth_phase2(corners, ud_edges);  
        for(togo2 = dist2; togo2 <= togo2_limit; togo2++)      //限制步数的深搜，IDA*
            search2(corners, ud_edges, slice, dist2, togo2)
    }

Ⅱ 允许的转动

前面很多代码中都写道所有 18 中转动操作中允许的转动，什么意思呢？还是用例子来说明，比如我这一步采用 L ，则下一次转动就不该使用 { L , L 2 , L 3 } \{L, L2,L3\} {L,L2,L3} ，因为采用的话，总有更短的转动序列可以替代，比如说 $LL=L2,LL2=L3,LL3=0,L2L3=L$ 。所以 { L , L 2 , L 3 } \{L, L2,L3\} {L,L2,L3} 后面不应使用 { L , L 2 , L 3 } \{L, L2,L3\} {L,L2,L3} 。

魔方转动一般步满足交换律的，比如 $FR\ne RF$ ，但是相对的两层转动是满足交换律的，比如 $LR=RL$ ，所以如果 L 类转动后面可以跟 R 类转动，那么就要禁止 R 类的转动后面跟 L 类的转动。

另外每个阶段都有限制的转动，特别是在阶段二，只能使用 { U , D , L 2 , R 2 , F 2 , B 2 } \{U,D,L2,R2,F2,B2\} {U,D,L2,R2,F2,B2}

而在阶段一中，如果 $dist=0ANDtogo1\ne 0$ ，说明是要在阶段一种寻找次优解，以达到总体更优解。这时，即使处于阶段一，但我们禁止使用 { U 1 , U 2 , U 3 , D 1 , D 2 , D 3 , L 2 , R 2 , F 2 , B 2 } \{U1,U2,U3,D1,D2,D3,L2,R2,F2,B2\} {U1,U2,U3,D1,D2,D3,L2,R2,F2,B2} ，这是阶段二使用的转动，但是我们禁止使用。为什么呢？比如一个转动序列 M 可以让初始状态到达目标状态，那么 MU，MR2 等转动序列也能使初始状态到达目标状态，但是这种阶段一的次优解并不能让我们找到更优解。

Ⅲ 状态/变换的逆元

实在不知道去什么名字好，干脆就用数学里面的术语，前面我们求过对称变换的逆变换，这里我们要求一个魔方状态逆元或者说一个普通变换的逆变换。这里的实现方式不同于对称变换那里两两验证，直接来看码：

CubieCube inverse_cubiecube(CubieCube cc){
    CubieCube inv;
	for (Edge edg = UR;edg <= BR; edg++) 		  //棱边位置
        inv.eo[cc.eo[edg].e].e = edg;
	for (Edge edg = UR;edg <= BR; edg++)          //棱边方向
        inv.eo[edg].o = cc.eo[inv.eo[edg].e].o;

	for (Corner crn = URF;crn <= DRB; crn++)      //角块位置
        inv.co[cc.co[crn].c].c = crn;
	for (Corner crn = URF; crn <= DRB;crn++)      //角块方向
        inv.co[crn].o = (3 - cc.co[ccInv.co[crn].c].o) % 3   //角块方向要发生变化
}

里面的点点有一点点多，不要绕晕了，模拟一下怎么变化的应该很好理解。说一下为什么棱边的方向没有变化，而角块的方向要变化。究其原因与角块有 3 个方向和对块边的定义有关，对于角块如果我们只从扭转与未扭转的角度来看，角块方向也是没有发生变化的，可以简单代几个数试一试，最后一个式子作用只是让 1 变 2，2 变 1，0 还是 0。如果还是不太清楚可以对一个魔方分别做相逆操作，然后根据方向定义来观察一下块边的方向如何变化的。

Ⅳ 复合变换/转动时角块方向问题

由于加入了镜像变换，对角块的方向进行了重新定义，镜像的方向取值范围为 $[3,5]$ 。这里顺时针扭转变为 5 ，逆时针扭转变为 4 ，这里我看国内国外有关资料都没有指出来，而根据实际代码操作，魔方转动模拟出来的结果应该是这样。

corner_multiply 函数是来处理变换时角块的位置和方向变化的函数，上述我们只处理了普通状态，加入镜像后会有一些变动，这里完善，还是直接来看码：

void cornerMultiply(const CubieCube *a, const CubieCube *b, CubieCube *ab){ 
	for (Corner crn = URF; crn <= DRB; crn++){
		ab->co[crn].c = a->co[b->co[crn].c].c;     //角块位置变化
		char oriA = a->co[b->co[crn].c].o;         //原本方向
        char oriB = b->co[crn].o);                 //方向怎么变化
        /***根据普通还是镜像分别讨论***/
        if (oriA < 3 && oriB < 3){            //a, b都处于普通状态
            ori = oriA + oriB;      
            if (ori >= 3) ori -= 3;           //结果普通状态
	    }else if (oriA < 3 && oriB >= 3){     //b处于镜像状态
            ori = oriA + oriB;
            if (ori >= 6) ori-=3;	          //结果镜像状态
		}else if (oriA >= 3 && oriB < 3){     //a处于镜像状态
            ori = oriA - oriB;
            if (ori < 3) ori += 3;		      //结果镜像状态
		}else if (oriA >= 3 && oriB >= 3){    //a，b都处于镜像状态
            ori = oriA - oriB;
            if (ori < 0) ori += 3;		      //结果普通状态
		}
    }
}

镜像状态做普通变换，普通状态下做镜像变换，结果为镜像状态，镜像状态下做镜像变化相当于有变回去了所以结果为普通状态，应该很好理解。

主要是上述对于不同情况下 oriA，oriB 的加减操作和最后结果的 ori 加 3 还是 减 3 可能很迷惑，嘿，很迷惑但我不细说，主要是因为关于这个方向的问题实在是不太好表述，要想搞清楚最好的方法就是根据简单的镜像变换，普通变换模拟一下，然后拿一个实物魔方作为参考，看看是不是这样。要详说模拟操作的话内容大多就是离散里面置换操作，重复且枯燥，所以不说了，感兴趣的自己去模拟一下吧。

Ⅴ 多线程多方向搜索提高效率

可以使用多线程并行搜索，并行搜索不是对同一个方向同一种情况多次搜索，而是搜索不同的方向。比如搜索另外的轴向，搜索逆反状态。搜索另外轴向意思就是在搜索之前先对魔方进行 S_URF3 对称变换(不一定是绕 URF3 角转动，其他也行)，同理搜索逆反状态呢就是先对魔方进行逆变换。

记住初始使用的什么变换，搜索结束后得到的解法再做逆变换即可得到原本状态的解法。实验证明，开 6 个线程，三个轴向和逆变换的不同组合效率是最高的。

九、结尾参考

科先巴的二阶段算法大概就是这样，基本上囊括了所有的重难点吧。纵观下来，算法核心在于各种数据结构的构建，要创建各种各样的表，特别是剪枝表的建立，直接影响到搜索性能和解法优劣。

本文参考的外文资料，感兴趣的可以去看看原版资料

http://kociemba.org/cube.htm，科先巴的个人网站，详细的讲述二阶段算法，里面也有程序和源码，看源码的建议看 python 版本的，各个部分都十分清晰明了。

https://www.jaapsch.net/puzzles/，国外一大神的网站，里面有各种魔方的理论，程序源码。

http://www.cube20.org，宣布上帝之数是 20 的一个网站，下面也有一个关于二阶段算法的源码和相应的 pdf 讲解，看完理解透应该就成神了，我是坚持看完之后差点没背过去，感兴趣想变强的可以试试。

https://en.wikipedia.org/wiki/Optimal_solutions_for_Rubik’s_Cube，Wiki百科

https://www.speedsolving.com/wiki/index.php/Kociemba’s_Algorithm，Wiki百科

https://ruwix.com/the-rubiks-cube 也是很不错的一个网站，里面有各种理论也有在线的程序等等很全

The Diameter of the Rubik’s Cube Group Is Twenty，魔方领域影响力最大的一篇论文，得出结论上帝之数为 20，值得一看。

还有一些就不举出来了，上面的都是个人认为比较好的一些网站资源，想看原版的理论详解可以去看看，都是英文的，国内在这方的资料文献的确太少太少了，几乎没有，在下不才，也希望尽一点绵薄之力，同样的有什么错误还请批评指正。有些网站还有资源是需要梯子的，没有的朋友可以在后台回复魔方，我把上述的一些资料源码还有程序放在里面的，方便下载。

本文关于科先巴的二阶段算法就到这里了，科先巴算法一般是得不到最优解的，有没有什么算法能够每次都得到最优解呢，答案是有的，就是篇一的时候也提到过的上帝算法——Krof’s 算法，也没什么神秘的，是基于科先巴的二阶段算法实现的，理解了二阶段算法，应该就很好理解，咱们下篇再详细叙述。

最后再谈一点题外话，对于二阶段算法有感而发吧，其实这二阶段算法不仅才计算机领域魔方领域能给我们带来启发，平时生活学习也有帮助，那句话怎么说来着，先定个小目标，一步步解决。还记得文章封面，关于二阶段算法的那个生动比喻吗？漫无目的航行想要彼岸那是遥遥无期，不如先到特定水域完成一个个目标，最终成功到达彼岸。虽然可能得不到最优解，但人生这条路谁又能够一帆风顺的走直线呢？只要刻苦努力奋斗，争取人生每个阶段都得到相应的最优解，合起来得到一个较优解那也是很不错的嘛对吧。