本文深入探讨了幻方的构造方法，包括奇数阶、双偶数阶和单偶数阶幻方的生成算法。针对不同阶数的幻方，提出了罗伯法、海尔法和斯特雷奇法。此外，通过剪枝优化的深度优先搜索算法，能够在较短的时间内构造特定阶数的幻方。同时，研究了如何通过行列交换、旋转等操作增加幻方的多样性。对于高阶幻方，虽然无法穷举所有可能，但通过特定的幻方构造技巧，能够生成大量幻方实例。
题目
幻方是一种将数字安排在正方形格子中，使每行、列和对角线上的数字和都相等的方法。
- 证明：如果一个n阶幻方存在的话，所讨论的和一定等于n*(n^2+1)/2
- 设计一个穷举算法，生成阶数为n的所有幻方
- 更好的生成幻方的算法
题解
证明：如果一个n阶幻方存在的话，所讨论的和一定等于n*(n^2+1)/2
对于一个n阶的幻方，其中的数字为1 ~ n^2，总和为(1+n ^ 2)*n ^ 2/2，每一行或者每一列或者每一条对角线的元素个数都为n，则元素的和为
(1+n^2)*n^2/2 / n = (1+n^2)*n/2
得证
设计一个穷举算法，生成阶数为n的所有幻方
方案一：
对于一个n阶方阵，穷举所有的填数方法，使用DFS来实现，时间复杂度O((n ^ 2) ^ (n ^ 2)))
display(): 自定义函数，用于打印方阵
void display(int n){
cout << "--" << "第" << count_matrix << "个" << "----" << endl;
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++){
cout << matrix[i][j] << " ";
}
cout << endl;
}
cout << "----------" << endl;
}
match():自定义函数，用来判断某种填数方案是否符合幻方的条件
bool match(int n){
int count;
int sum = (1 + n * n) * n / 2;
//行
for (int i = 0; i < n; i++){
count = 0;
for (int j = 0; j < n; j++){
count = count + matrix[i][j];
}
if (count != sum){
return false;
}
}
//列
for (int i = 0; i < n; i++){
count = 0;
for (int j = 0; j < n; j++){
count = count + matrix[j][i];
}
if (count != sum){
return false;
}
}
//左上到右下对角线
count = 0;
for (int i = 0; i < n; i++){
count = count + matrix[i][i];
}
if (count != sum){
return false;
}
//右上到左小对角线
count = 0;
for (int i = 0; i < n; i++){
count = count + matrix[i][n-1-i];
}
if (count != sum){
return false;
}
return true;
}
number_used[]: int型一维数组，用于记录某个数是否已经被填入方阵中，在进行回溯时，需要恢复其状态
matrix[][]: int型二维数组，记录方阵的填数情况
number_max: 方阵中所能填的最大的数，值为n^2
count_matrix: 统计幻方的数量
void dfs1(int k){
if (k == numer_max){
if (match(n)){
count_matrix = count_matrix + 1;
display(n);
}
}
else{
for (int i = 1; i <= numer_max; i++){
if (number_used[i] == 0){
matrix[k / n][k % n] = i;
number_used[i] = 1;
dfs1(k + 1);
number_used[i] = 0;
}
}
}
}
运行结果：
（不存在二阶幻方，从三阶幻方开始实验）
可以在0.027秒的时间内求出所有的3阶幻方
但无力求出哪怕是一种4阶幻方的方案
方案二：
对方案一的DFS进行剪枝，将判断方阵是否符合幻方特征的过程拆解开来，放入到方阵构造的过程中去，这样可以尽早地让程序退出一条错误的路径（而不是像方案一那样，尽管很早就知道这条路是错误的，依旧要将整个方阵填完）
SUM: 行、列、对角线元素的和
第一处if: 每当一行的数填完后，统计这一行的数的和，若在统计过程中和已经大于SUM或者统计结束后和不等于SUM，表明填法错误，需要回溯。
第二处if: 求当前将要填的位置的前一 个位置所在列的元素的和，若行数不 等于方阵的阶但和大于SUM，则填 法错误，需要回溯;若行数等于方阵 的阶但和不等于SUM，则填法错 误，需要回溯。
当所有数都填完后(k == number_max)，调用match()函数 来检验方阵是否满足幻方的条件。
void dfs2(int k){
int SUM = (1 + n * n) * n / 2;
int temp = 0;
if ((k % n == 0) && (k != 0)){
temp = 0;
for (int i = (k / n - 1) * n; i < k; i++){
temp = temp + matrix[i / n][i % n];
if (temp > SUM){
return;
}
}
if (temp != SUM){
return;
}
}
if (k != 0){
temp = 0;
int kk = k - 1;
for (int i = 0; i * n + kk % n <= kk; i++){
temp = temp + matrix[i][kk % n];
if (i == n - 1){
if (temp != SUM){
return;
}
}
}
if (temp > SUM){
return;
}
}
if (k == numer_max){
if (match(n)){
count_matrix = count_matrix + 1;
display(n);
}
}
else{
for (int i = 1; i <= numer_max; i++){
if (number_used[i] == 0){
matrix[k / n][k % n] = i;
number_used[i] = 1;
dfs2(k + 1);
number_used[i] = 0;
}
}
}
}
运行结果：
剪枝过后的DFS效率有了显著的提 升，搜索所有可能的3阶幻方用时为 0.0008秒(未优化的DFS需要0.027 秒)。
剪枝过后的DFS可以在能够接受的时 间范围内求得所有的4阶幻方，所需 时间为98秒。
更好的生成幻方的算法
查询相关资料，学习了一些著名的幻方构造法
首先需要将幻方的类型进行分类，分为奇数阶、双偶数阶(可以被4整 除)、单偶数阶(不可以被4整除但可以被2整除)这三种幻方。
奇数阶幻方：罗伯法
把1放在第一行正中，并按以下规律排列剩下的(n*n-1)个数:
- 每个数放在前一个数的右上一格
- 若超出顶行，则放在底行，仍要在右边一列
- 若超出最右列，则放在最左列，仍要在上一行
- 若超出顶行且超出最右列，则放在底行且最左列
- 如果下一格已有数，则放在前一个数的下一行同一列的格内。
//奇数阶幻方：罗伯法（连续摆数法）
void CreateOddMagicSquare(int n){
int number = 1;
int x = 0;
int y = n / 2;
while (number <= n * n){
matrix[x][y] = number;
int x_try = x;
int y_try = y;
x_try--;
y_try++;
if (x_try < 0){
x_try = x_try + n;
}
if (y_try == n){
y_try = n - y_try;
}
if (matrix[x_try][y_try] == 0){
x = x_try;
y = y_try;
}
else{
x = x + 1;
if (x == n){
x = 0;
}
}
number++;
}
}
matrix[][]:int型的二维数组，存储方阵的填数情况
先判断是否超出顶行
再判断是否超出最右列
之后判断要放的位置是否已有数
时间复杂度O(n^2)
双偶数阶幻方：海尔法
- 先把数字按顺序填，然后将方阵均匀分成4块
- 将每个小方阵对⻆线上的数字，换成和它互补的数
互补:在n阶幻方中，和为n^2+1的两个数为互补数
//双偶数阶幻方：海尔法
void CreateDoubleEvenMagicSquare(int n){
int two_number_sum = 1 + n * n;
int number = 1;
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++){
matrix[i][j] = number;
number++;
}
}
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++){
if (i % 4 == 0 && abs(i - j) % 4 == 0){
for (int k = 0; k < 4; k ++){
matrix[i + k][j + k] = abs(two_number_sum - matrix[i + k][j + k]);
}
}
else if (i % 4 == 3 && (i + j) % 4 == 3){
for (int k = 0; k < 4; k++){
matrix[i - k][j + k] = abs(two_number_sum - matrix[i - k][j + k]);
}
}
}
}
}
two_number_sum:两个互补的数的和
首先将方阵按顺序填数
之后对分割后的方阵对⻆线上的数进行处理
时间复杂度O(n^2)
单偶数阶幻方：斯特雷奇法
单偶数阶幻方就是当n不可以被4整 除时的偶数阶幻方，即4k+2阶幻 方。半阶数m = n/2。
- 把n阶矩阵分为A，B，C，D四个象限(分别对应象限2，1，3，4)，这样每一个象限肯定是奇数阶，阶数为m。按照奇数阶幻方的构造方法用1~(n/2)^2在A象限填写数字，并用A象限初始化其他象限，其中B=A+2m ^ 2，C=A+3m ^ 2，D=A+m^2。如图一。
- 从A象限的中间行、中间格开始，按自左向右的方向，标出k格。A象限的其它行则标出最左边的k格。将这些格，和C象限相对位置上的数，互换位置。如图2。
- 在B象限任一行的中间格，自右向左，标出k-1列。(有一个例外，6阶幻方由于k-1=0，所以不用做B、D象限的数据交换操作)将B象限标出的这些数，和D象限相对位置上的数进行交换，就形成幻方。
//单偶数阶幻方：斯特雷奇法
void CreateSingleEvenMagicSquare(int n){
int k = (n - 2) / 4;
int m = n / 2;
CreateOddMagicSquare(m);
for (int i = 0; i < m; i++){
for (int j = 0; j < m; j++){
matrix[i][j + m] = matrix[i][j] + 2 * m * m; // 第一象限
matrix[i+m][j] = matrix[i][j] + 3 * m * m; //第三象限
matrix[i+m][j+m] = matrix[i][j] + m * m; //第四象限
}
}
for (int i = 0; i < m; i++){
if (i != m / 2){
for (int j = 0; j < k; j++){
int temp = matrix[i][j];
matrix[i][j] = matrix[i + m][j];
matrix[i + m][j] = temp;
}
}
}
for (int i = 0; i < k; i++){
int temp = matrix[m / 2][m / 2 + i];
matrix[m / 2][m / 2 + i] = matrix[m / 2 + m][m / 2 + i];
matrix[m / 2 + m][m / 2 + i] = temp;
}
for (int i = 0; i < m; i++){
for (int j = 0; j < k - 1; j++){
int temp = matrix[i][m + m / 2 + j];
matrix[i][m + m / 2 + j] = matrix[i + m][m + m / 2 + j];
matrix[i + m][m + m / 2 + j] = temp;
}
}
}
先构造一个n/2阶的奇数幻方
通过这个奇数幻方构造剩下的三个小方阵
根据斯特雷奇法对方阵中的数进行处理
时间复杂度O(n^2)
由此我们已经能够构造任意阶的幻方
但并不能获得某一阶所有的幻方
需要在获得的幻方的基础上，进行一些操作使得幻方的数量增加
双偶数阶幻方个数暴增
幻方中的每两个关于幻方中心中心对称的两个数均为互补数，因而关于幻方中 心中心对称的两列或者两行上的对应位置数字也是互补数。因此关于幻方对称 轴映射交换任意行列都是可以构成幻方的。
行列交换
行列映射交换的方式有以下四种(由于行列的交换都是一样的，这里以列交换进行分析):左右对称交换、左右交叉交换、单侧映射交换、混合交换(即左右对称交换、左右交叉交换和单侧交换混合)
当阶数很大的时候，混合交换的情况就会变得非常复杂，似乎是无法找到所有的情况。经过分析，所有情况都可以通过以下两个步骤获得:
第一步:确定左右交换的列(行)，进行左右对称交换;
第二步:对对称交换后的幻方进行单侧映射交换。
这两个步骤其实就是排列组合的过程。第一步，利用组合对左侧找到所有组合情况，即所有左右交换的情况，包含交换0列，1列…n/2列，对组合结果进行左右对称交换，可得到2^(n/2)种。 第二步，对左侧进行全排列，找到所有排列的情况，即所有单侧映射交换的情况，进行单侧映射交换，可得到(n/2)!种。
组合问题的代码实现
//求组合
void Combination(int n){
count_index_Combination = 1;
for (int i = 0; i < n; i++){
count_index_Combination = count_index_Combination * 2;
}
for (int i = 0; i <= count_index_Combination - 1; i++){
int temp = n-1;
int index = i;
do{
DoubleEven_index_Combination[temp][i] = index % 2;
temp--;
index = index / 2;
}while(index != 0);
}
}
count_index_Combination:记录组合的数量
DoubleEven_index_Combination[][]:int型的二维数组，第一维下标记录对应的方阵的行(列)标，第二位下标记录组合数，当数组的值为1，表示该元素被选中将进行左右交换
由于选中为1，未选中为0，将问题转换为有几个不同的n/2位二进制数的问题，因此遍历0~2^(n/2)-1即可，每次对2取余数，再整除2
全排列的代码实现
//求全排列
void FullPermutation(int begin, int end){
if (begin == end - 1){
for (int i = 0; i < end; i ++){
DoubleEven_index_FullPermutation[i][count_index_FullPermutation] = DoubleEven_index_temp_FullPermutation[i];
}
count_index_FullPermutation++;
}else{
for (int i = begin; i < end; i++){
int temp = DoubleEven_index_temp_FullPermutation[begin];
DoubleEven_index_temp_FullPermutation[begin] = DoubleEven_index_temp_FullPermutation[i];
DoubleEven_index_temp_FullPermutation[i] = temp;
FullPermutation(begin + 1, end);
temp = DoubleEven_index_temp_FullPermutation[begin];
DoubleEven_index_temp_FullPermutation[begin] = DoubleEven_index_temp_FullPermutation[i];
DoubleEven_index_temp_FullPermutation[i] = temp;
}
}
}
DoubleEven_index_temp_FullPermutation[]:存储一种全排列方案，初始化时按顺序填数
DoubleEven_index_FullPermutation[][]:存储对于某个方阵所有的全排列方案
使用DFS的方式找到所有的全排列方案(交换元素)
旋转
四个矩阵，后面一个是前面一个顺时针旋转90度获得的。
通过行列交换可以由a图得到c图，但无法获得b、c图，即无法使得同列关系的数字1，3变为同行的关系。但是可以通过a、c图顺时针旋转分别得到b、d图。因此在由交换行列获得的幻方的基础上，可以再通过旋转获得更多的幻方。
双偶数阶幻方预期个数
全排列：(n/2)!
组合：2^(n/2)
行（列）：平方
旋转：乘2
结果：
0.005秒求出128个4阶幻方，38秒求出294912个8阶幻方，对于12阶幻方则无法在可接受时间内求出所有该算法能够构造的幻方。
奇数阶幻方个数暴增
对于奇数阶幻方，发现当阶数为素数时，可以通过以下方法在原有的构造好的奇数阶幻方的基础上再得到n-3个构造方案(这些构造方案无法通过行、列变换获得)
对于n/2+1 ~ n列，每列相对于左边的一列向上移动k行(k为1~n-3中的值)
对于n/2 ~ 1列，每列相对于右边的一列向下移动k行(k为1~n-3中的值)
奇数阶幻方也可以和双偶数阶幻方一样，利用行、列交换和旋转来增加方案数，方法类似，不再赘述。
奇数阶幻方预期方案数
结果：
和双偶数阶幻方预期方案数相比，区别在于增加了一项(n-2)^(isPrimeNum(n))，用于计算素数阶幻方通过对称移动行、列增加的方案数，并且将这一项放在了平方的括号之内，意味着这个操作可以同时对行、列起作用，但实际实验发现部分行进行此操作生成的方案是错误的。
单偶数阶幻方个数暴增
对于单偶数阶幻方的构造过程，还可以使用矩阵的加法及数乘来进行实现。其中由0、1、2、3组成的矩阵我们称之为叠加矩阵，另一个由奇数幻方生成的因数矩阵我们称之为源矩阵，最右边的矩阵我们称之为生成矩阵。
因此我们可以对四个子矩阵都采用求奇数阶幻方的方法来构造，之后使用叠加矩阵使得值符合大矩阵的填法。
再加上使用旋转的技巧，可以使得数量再次乘4(因为大矩阵没有进行行列交换)
单偶数阶幻方预期方案数
由于实际编码实现过程中，实现四个子方阵皆独立生成代码过于庞大，因此仅实现了子方阵A的幻方生成，剩下三个子方阵由斯特雷奇法构造。
因此实际预期方案数为(退化为n/2阶奇数幻方的方案数):
结果：
由于构造素数阶幻方暴增方法中部分生成幻方无效，因此实际的方案数会小于预期方案数
总结
使用剪枝后的DFS可以在一分钟内求出所有的4阶幻方。
通过查阅资料找到的方案，可以很好地解决构造特定阶数幻方的问题，但无法求出某阶幻方的所有可能方案。
对于构造幻方的方法，三种类别幻方的构造算法的时间复杂度均为O(n^2)
因此在一分钟内，理论上可以生成非常高阶的幻方。
（本作业为算法课的大作业，在一番努力后，最终除了爆搜，也没有找到更好的能求出某一阶所有幻方方案的算法，是一大遗憾，此文的算法和方案仅供参考——可能会有不少错误）
代码
DFS代码：
#include <iostream>
#include <string.h>
#include <math.h>
#include <ctime>
using namespace std;
clock_t start_time, end_time;
int n;
int numer_max;
int matrix[10][10];
int number_used[100];
int count_matrix;
void init(){
for (int i = 0; i < 100; i++){
number_used[i] = 0;
}
}
void display(int n){
cout << "--" << "第" << count_matrix << "个" << "----" << endl;
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++){
cout << matrix[i][j] << " ";
}
cout << endl;
}
cout << "----------" << endl;
}
bool match(int n){
int count;
int sum = (1 + n * n) * n / 2;
//行
for (int i = 0; i < n; i++){
count = 0;
for (int j = 0; j < n; j++){
count = count + matrix[i][j];
}
if (count != sum){
return false;
}
}
//列
for (int i = 0; i < n; i++){
count = 0;
for (int j = 0; j < n; j++){
count = count + matrix[j][i];
}
if (count != sum){
return false;
}
}
//左上到右下对角线
count = 0;
for (int i = 0; i < n; i++){
count = count + matrix[i][i];
}
if (count != sum){
return false;
}
//右上到左小对角线
count = 0;
for (int i = 0; i < n; i++){
count = count + matrix[i][n-1-i];
}
if (count != sum){
return false;
}
return true;
}
void dfs1(int k){
if (k == numer_max){
if (match(n)){
count_matrix = count_matrix + 1;
display(n);
}
}
else{
for (int i = 1; i <= numer_max; i++){
if (number_used[i] == 0){
matrix[k / n][k % n] = i;
number_used[i] = 1;
dfs1(k + 1);
number_used[i] = 0;
}
}
}
}
void dfs2(int k){
int SUM = (1 + n * n) * n / 2;
int temp = 0;
if ((k % n == 0) && (k != 0)){
temp = 0;
for (int i = (k / n - 1) * n; i < k; i++){
temp = temp + matrix[i / n][i % n];
if (temp > SUM){
return;
}
}
if (temp != SUM){
return;
}
}
if (k != 0){
temp = 0;
int kk = k - 1;
for (int i = 0; i * n + kk % n <= kk; i++){
temp = temp + matrix[i][kk % n];
if (i == n - 1){
if (temp != SUM){
return;
}
}
}
if (temp > SUM){
return;
}
}
if (k == numer_max){
if (match(n)){
count_matrix = count_matrix + 1;
display(n);
}
}
else{
for (int i = 1; i <= numer_max; i++){
if (number_used[i] == 0){
matrix[k / n][k % n] = i;
number_used[i] = 1;
dfs2(k + 1);
number_used[i] = 0;
}
}
}
}
int main(int argc, const char * argv[]) {
start_time = clock();
count_matrix = 0;
n = 3;
init();
numer_max = pow(n, 2);
dfs2(0);
cout << "共" << count_matrix << "个" << endl;
end_time = clock();
double total_time = (double)(end_time - start_time) / CLOCKS_PER_SEC;
cout << "耗时：" << total_time << endl;
return 0;
}
幻方构造生成代码：
#include <iostream>
#include <string.h>
#include <math.h>
#include <ctime>
using namespace std;
clock_t start_time, end_time;
#define MAXN 10000
int matrix[100][100];
int DoubleEven_index_FullPermutation[10][100];
int DoubleEven_index_Combination[10][100];
int DoubleEven_index_temp_FullPermutation[10];
int count_index_FullPermutation;
int count_index_Combination;
int count_matrix_number;
void init(int n){
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++){
matrix[i][j] = 0;
}
}
count_matrix_number = 0;
}
void DoubleEven_init(){
for (int i = 0; i < 10; i++){
for (int j = 0; j < 100; j++){
DoubleEven_index_FullPermutation[i][j] = i;
DoubleEven_index_Combination[i][j] = 0;
}
}
for (int i = 0; i < 10; i++){
DoubleEven_index_temp_FullPermutation[i] = i;
}
count_index_FullPermutation = 0;
count_index_Combination = 0;
}
void display(int n, int matrix[][100]){
cout << "----------" << endl;
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++){
cout << matrix[i][j] << " ";
}
cout << endl;
}
cout << "----------" << endl;
}
bool match(int n, int matrix[][100]){
int count;
int sum = (1 + n * n) * n / 2;
//行
for (int i = 0; i < n; i++){
count = 0;
for (int j = 0; j < n; j++){
count = count + matrix[i][j];
}
if (count != sum){
return false;
}
}
//列
for (int i = 0; i < n; i++){
count = 0;
for (int j = 0; j < n; j++){
count = count + matrix[j][i];
}
if (count != sum){
return false;
}
}
//左上到右下对角线
count = 0;
for (int i = 0; i < n; i++){
count = count + matrix[i][i];
}
if (count != sum){
return false;
}
//右上到左小对角线
count = 0;
for (int i = 0; i < n; i++){
count = count + matrix[i][n-1-i];
}
if (count != sum){
return false;
}
return true;
}
//奇数阶幻方：罗伯法（连续摆数法）
void CreateOddMagicSquare(int n){
int number = 1;
int x = 0;
int y = n / 2;
while (number <= n * n){
matrix[x][y] = number;
int x_try = x;
int y_try = y;
x_try--;
y_try++;
if (x_try < 0){
x_try = x_try + n;
}
if (y_try == n){
y_try = n - y_try;
}
if (matrix[x_try][y_try] == 0){
x = x_try;
y = y_try;
}
else{
x = x + 1;
if (x == n){
x = 0;
}
}
number++;
}
}
//双偶数阶幻方：海尔法
void CreateDoubleEvenMagicSquare(int n){
int two_number_sum = 1 + n * n;
int number = 1;
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++){
matrix[i][j] = number;
number++;
}
}
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++){
if (i % 4 == 0 && abs(i - j) % 4 == 0){
for (int k = 0; k < 4; k ++){
matrix[i + k][j + k] = abs(two_number_sum - matrix[i + k][j + k]);
}
}
else if (i % 4 == 3 && (i + j) % 4 == 3){
for (int k = 0; k < 4; k++){
matrix[i - k][j + k] = abs(two_number_sum - matrix[i - k][j + k]);
}
}
}
}
}
//单偶数阶幻方：斯特雷奇法
void CreateSingleEvenMagicSquare(int n){
int k = (n - 2) / 4;
int m = n / 2;
CreateOddMagicSquare(m);
for (int i = 0; i < m; i++){
for (int j = 0; j < m; j++){
matrix[i][j + m] = matrix[i][j] + 2 * m * m; // 第一象限
matrix[i+m][j] = matrix[i][j] + 3 * m * m; //第三象限
matrix[i+m][j+m] = matrix[i][j] + m * m; //第四象限
}
}
for (int i = 0; i < m; i++){
if (i != m / 2){
for (int j = 0; j < k; j++){
int temp = matrix[i][j];
matrix[i][j] = matrix[i + m][j];
matrix[i + m][j] = temp;
}
}
}
for (int i = 0; i < k; i++){
int temp = matrix[m / 2][m / 2 + i];
matrix[m / 2][m / 2 + i] = matrix[m / 2 + m][m / 2 + i];
matrix[m / 2 + m][m / 2 + i] = temp;
}
for (int i = 0; i < m; i++){
for (int j = 0; j < k - 1; j++){
int temp = matrix[i][m + m / 2 + j];
matrix[i][m + m / 2 + j] = matrix[i + m][m + m / 2 + j];
matrix[i + m][m + m / 2 + j] = temp;
}
}
}
//求全排列
void FullPermutation(int begin, int end){
if (begin == end - 1){
for (int i = 0; i < end; i ++){
DoubleEven_index_FullPermutation[i][count_index_FullPermutation] = DoubleEven_index_temp_FullPermutation[i];
}
count_index_FullPermutation++;
}else{
for (int i = begin; i < end; i++){
int temp = DoubleEven_index_temp_FullPermutation[begin];
DoubleEven_index_temp_FullPermutation[begin] = DoubleEven_index_temp_FullPermutation[i];
DoubleEven_index_temp_FullPermutation[i] = temp;
FullPermutation(begin + 1, end);
temp = DoubleEven_index_temp_FullPermutation[begin];
DoubleEven_index_temp_FullPermutation[begin] = DoubleEven_index_temp_FullPermutation[i];
DoubleEven_index_temp_FullPermutation[i] = temp;
}
}
}
//求组合
void Combination(int n){
count_index_Combination = 1;
for (int i = 0; i < n; i++){
count_index_Combination = count_index_Combination * 2;
}
for (int i = 0; i <= count_index_Combination - 1; i++){
int temp = n-1;
int index = i;
do{
DoubleEven_index_Combination[temp][i] = index % 2;
temp--;
index = index / 2;
}while(index != 0);
}
}
//求双偶幻方的所有方案
void displayDoubleEvenMagicSquare(int n){
for (int i = 0; i < count_index_FullPermutation; i++){
int temp_matrix1[100][100];
int temp_matrix2[100][100];
for (int j = 0; j < n / 2; j++){
for (int k = 0; k < n; k++){
temp_matrix1[j][k] = matrix[DoubleEven_index_temp_FullPermutation[j]][k];
temp_matrix1[n - 1 - j][k] = matrix[n - 1 - DoubleEven_index_temp_FullPermutation[j]][k];
}
}
for (int j = 0; j < count_index_FullPermutation; j++){
for (int k = 0; k < n / 2; k++){
for (int l = 0; l < n; l++){
temp_matrix2[l][k] = temp_matrix1[l][DoubleEven_index_temp_FullPermutation[k]];
temp_matrix2[l][n - 1 - k] = temp_matrix1[l][n - 1 - DoubleEven_index_temp_FullPermutation[k]];
}
}
for (int l = 0; l < count_index_Combination; l++){
for (int p = 0; p < n / 2; p++){
if (DoubleEven_index_Combination[p][l] == 1){
int ttemp[1000];
for (int q = 0; q < n; q++){
ttemp[q] = temp_matrix2[p][q];
}
for (int q = 0; q < n; q++){
temp_matrix2[p][q] = temp_matrix2[n - 1 - p][q];
}
for (int q = 0; q < n; q++){
temp_matrix2[n - 1 - p][q] = ttemp[q];
}
}
}
for (int p = 0; p < count_index_Combination; p++){
for (int q = 0; q < n / 2; q++){
if (DoubleEven_index_Combination[q][p] == 1){
int ttemp[1000];
for (int r = 0; r < n; r++){
ttemp[r] = temp_matrix2[r][q];
}
for (int r = 0; r < n; r++){
temp_matrix2[r][q] = temp_matrix2[r][n - 1 - q];
}
for (int r = 0; r < n; r++){
temp_matrix2[r][n - 1 - q] = ttemp[r];
}
}
}
count_matrix_number++;
cout << "-------" << count_matrix_number << "-------" << endl;
for (int ii = 0; ii < n; ii++){
for (int jj = 0; jj < n; jj++){
cout << temp_matrix2[ii][jj] << " ";
}
cout << endl;
}
cout << "---------------" << endl;
count_matrix_number++;
cout << "-------" << count_matrix_number << "-------" << endl;
for (int ii = 0; ii < n; ii++){
for (int jj = 0; jj < n; jj++){
cout << temp_matrix2[n - 1 - jj][ii] << " ";
}
cout << endl;
}
cout << "---------------" << endl;
}
}
}
}
cout << "共" << count_matrix_number << "个" << endl;
}
//判断是否为素数
bool judgeIfPrimeNumber(int x){
if (x == 2){
return true;
}
for (int i = 2; i * i <= x; i++){
if (x % i == 0){
return false;
}
}
return true;
}
void expandOddMagicSquare(int n, int matrix[][100]){
for (int i = 0; i < count_index_FullPermutation; i++){
int temp_matrix1[100][100];
int temp_matrix2[100][100];
for (int j = 0; j < n / 2; j++){
for (int k = 0; k < n; k++){
temp_matrix1[j][k] = matrix[DoubleEven_index_temp_FullPermutation[j]][k];
temp_matrix1[n - 1 - j][k] = matrix[n - 1 - DoubleEven_index_temp_FullPermutation[j]][k];
}
}
for (int j = 0; j < n; j++){
temp_matrix1[n / 2][j] = matrix[n / 2][j];
}
for (int j = 0; j < count_index_FullPermutation; j++){
for (int k = 0; k < n / 2; k++){
for (int l = 0; l < n; l++){
temp_matrix2[l][k] = temp_matrix1[l][DoubleEven_index_temp_FullPermutation[k]];
temp_matrix2[l][n - 1 - k] = temp_matrix1[l][n - 1 - DoubleEven_index_temp_FullPermutation[k]];
}
}
for (int k = 0; k < n; k++){
temp_matrix2[k][n / 2] = temp_matrix1[k][n / 2];
}
for (int l = 0; l < count_index_Combination; l++){
for (int p = 0; p < n / 2; p++){
if (DoubleEven_index_Combination[p][l] == 1){
int ttemp[1000];
for (int q = 0; q < n; q++){
ttemp[q] = temp_matrix2[p][q];
}
for (int q = 0; q < n; q++){
temp_matrix2[p][q] = temp_matrix2[n - 1 - p][q];
}
for (int q = 0; q < n; q++){
temp_matrix2[n - 1 - p][q] = ttemp[q];
}
}
}
for (int p = 0; p < count_index_Combination; p++){
for (int q = 0; q < n / 2; q++){
if (DoubleEven_index_Combination[q][p] == 1){
int ttemp[1000];
for (int r = 0; r < n; r++){
ttemp[r] = temp_matrix2[r][q];
}
for (int r = 0; r < n; r++){
temp_matrix2[r][q] = temp_matrix2[r][n - 1 - q];
}
for (int r = 0; r < n; r++){
temp_matrix2[r][n - 1 - q] = ttemp[r];
}
}
}
count_matrix_number++;
cout << "-------" << count_matrix_number << "-------" << endl;
for (int ii = 0; ii < n; ii++){
for (int jj = 0; jj < n; jj++){
cout << temp_matrix2[ii][jj] << " ";
}
cout << endl;
}
cout << "---------------" << endl;
count_matrix_number++;
cout << "-------" << count_matrix_number << "-------" << endl;
for (int ii = 0; ii < n; ii++){
for (int jj = 0; jj < n; jj++){
cout << temp_matrix2[n - 1 - jj][ii] << " ";
}
cout << endl;
}
cout << "---------------" << endl;
}
}
}
}
}
//求奇偶幻方的所有方案
void displayOddMagicSquare(int n){
if (judgeIfPrimeNumber(n)){
for (int i = 0; i <= n - 3; i++){
int temp_matrix[100][100];
for (int j = 0; j < n; j++){
for (int k = 0; k < n; k++){
temp_matrix[j][k] = matrix[j][k];
}
}
for (int j = 1; j <= n / 2; j++){
int temp[100];
for (int k = 0; k < n; k++){
temp[k] = temp_matrix[k][j + n / 2];
}
for (int k = 0; k < n; k++){
temp_matrix[k][j + n / 2] = temp[(k + i * j) % n];
}
for (int k = 0; k < n; k++){
temp[k] = temp_matrix[k][n / 2 - j];
}
for (int k = 0; k < n; k++){
temp_matrix[(k + i * j) % n][n / 2 - j] = temp[k];
}
}
for (int j = 0; j <= n - 3; j++){
for (int k = 1; k <= n / 2; k++){
int temp[100];
for (int p = 0; p < n; p++){
temp[p] = temp_matrix[k + n / 2][p];
}
for (int p = 0; p < n; p++){
temp_matrix[k + n / 2][p] = temp[(p + j * k) % n];
}
for (int p = 0; p < n; p++){
temp[p] = temp_matrix[n / 2 - k][p];
}
for (int p = 0; p < n; p++){
temp_matrix[n / 2 - k][(p + j * k) % n] = temp[p];
}
}
if (match(n, temp_matrix)){
expandOddMagicSquare(n, temp_matrix);
}
}
}
}else{
expandOddMagicSquare(n, matrix);
}
cout << "共" << count_matrix_number << "个" << endl;
}
void expandSingleEvenMagicSquare(int n, int matrix[][100]){
for (int i = 0; i < count_index_FullPermutation; i++){
int temp_matrix1[100][100];
int temp_matrix2[100][100];
for (int j = 0; j < n / 2; j++){
for (int k = 0; k < n; k++){
temp_matrix1[j][k] = matrix[DoubleEven_index_temp_FullPermutation[j]][k];
temp_matrix1[n - 1 - j][k] = matrix[n - 1 - DoubleEven_index_temp_FullPermutation[j]][k];
}
}
for (int j = 0; j < n; j++){
temp_matrix1[n / 2][j] = matrix[n / 2][j];
}
for (int j = 0; j < count_index_FullPermutation; j++){
for (int k = 0; k < n / 2; k++){
for (int l = 0; l < n; l++){
temp_matrix2[l][k] = temp_matrix1[l][DoubleEven_index_temp_FullPermutation[k]];
temp_matrix2[l][n - 1 - k] = temp_matrix1[l][n - 1 - DoubleEven_index_temp_FullPermutation[k]];
}
}
for (int k = 0; k < n; k++){
temp_matrix2[k][n / 2] = temp_matrix1[k][n / 2];
}
for (int l = 0; l < count_index_Combination; l++){
for (int p = 0; p < n / 2; p++){
if (DoubleEven_index_Combination[p][l] == 1){
int ttemp[1000];
for (int q = 0; q < n; q++){
ttemp[q] = temp_matrix2[p][q];
}
for (int q = 0; q < n; q++){
temp_matrix2[p][q] = temp_matrix2[n - 1 - p][q];
}
for (int q = 0; q < n; q++){
temp_matrix2[n - 1 - p][q] = ttemp[q];
}
}
}
for (int p = 0; p < count_index_Combination; p++){
for (int q = 0; q < n / 2; q++){
if (DoubleEven_index_Combination[q][p] == 1){
int ttemp[1000];
for (int r = 0; r < n; r++){
ttemp[r] = temp_matrix2[r][q];
}
for (int r = 0; r < n; r++){
temp_matrix2[r][q] = temp_matrix2[r][n - 1 - q];
}
for (int r = 0; r < n; r++){
temp_matrix2[r][n - 1 - q] = ttemp[r];
}
}
}
int temp_matrix3[100][100];
for (int q = 0; q < n; q++){
for (int r = 0; r < n; r++){
temp_matrix3[q][r] = temp_matrix2[q][r];
}
}
int kk = (n * 2 - 2) / 4;
for (int q = 0; q < n; q++){
for (int r = 0; r < n; r++){
temp_matrix3[q][r + n] = temp_matrix3[q][r] + 2 * n * n; // 第一象限
temp_matrix3[q+n][r] = temp_matrix3[q][r] + 3 * n * n; //第三象限
temp_matrix3[q+n][r+n] = temp_matrix3[q][r] + n * n; //第四象限
}
}
for (int q = 0; q < n; q++){
if (q != n / 2){
for (int r = 0; r < kk; r++){
int temp = temp_matrix3[q][r];
temp_matrix3[q][r] = temp_matrix3[q + n][r];
temp_matrix3[q + n][r] = temp;
}
}
}
for (int q = 0; q < kk; q++){
int temp = temp_matrix3[n / 2][n / 2 + q];
temp_matrix3[n / 2][n / 2 + q] = temp_matrix3[n / 2 + n][n / 2 + q];
temp_matrix3[n / 2 + n][n / 2 + q] = temp;
}
for (int q = 0; q < n; q++){
for (int r = 0; r < kk - 1; r++){
int temp = temp_matrix3[q][n + n / 2 + r];
temp_matrix3[q][n + n / 2 + r] = temp_matrix3[q + n][n + n / 2 + r];
temp_matrix3[q + n][n + n / 2 + r] = temp;
}
}
count_matrix_number++;
cout << "-------" << count_matrix_number << "-------" << endl;
for (int ii = 0; ii < 2 * n; ii++){
for (int jj = 0; jj < 2 * n; jj++){
cout << temp_matrix3[ii][jj] << " ";
}
cout << endl;
}
cout << "---------------" << endl;
count_matrix_number++;
cout << "-------" << count_matrix_number << "-------" << endl;
for (int ii = 0; ii < 2 * n; ii++){
for (int jj = 0; jj < 2 * n; jj++){
cout << temp_matrix3[2 * n - 1 - jj][ii] << " ";
}
cout << endl;
}
cout << "---------------" << endl;
}
}
}
}
}
void displaySingleEvenMagicSquare(int n){
n = n / 2;
if (judgeIfPrimeNumber(n)){
for (int i = 0; i <= n - 3; i++){
int temp_matrix[100][100];
for (int j = 0; j < n; j++){
for (int k = 0; k < n; k++){
temp_matrix[j][k] = matrix[j][k];
}
}
for (int j = 1; j <= n / 2; j++){
int temp[100];
for (int k = 0; k < n; k++){
temp[k] = temp_matrix[k][j + n / 2];
}
for (int k = 0; k < n; k++){
temp_matrix[k][j + n / 2] = temp[(k + i * j) % n];
}
for (int k = 0; k < n; k++){
temp[k] = temp_matrix[k][n / 2 - j];
}
for (int k = 0; k < n; k++){
temp_matrix[(k + i * j) % n][n / 2 - j] = temp[k];
}
}
for (int j = 0; j <= n - 3; j++){
for (int k = 1; k <= n / 2; k++){
int temp[100];
for (int p = 0; p < n; p++){
temp[p] = temp_matrix[k + n / 2][p];
}
for (int p = 0; p < n; p++){
temp_matrix[k + n / 2][p] = temp[(p + j * k) % n];
}
for (int p = 0; p < n; p++){
temp[p] = temp_matrix[n / 2 - k][p];
}
for (int p = 0; p < n; p++){
temp_matrix[n / 2 - k][(p + j * k) % n] = temp[p];
}
}
if (match(n, temp_matrix)){
expandSingleEvenMagicSquare(n, temp_matrix);
}
}
}
}else{
expandSingleEvenMagicSquare(n, matrix);
}
cout << "共" << count_matrix_number << "个" << endl;
}
int main(int argc, const char * argv[]){
int n;
n = 10;
init(n);
if (n % 2 == 1){
start_time = clock();
// CreateOddMagicSquare(n);
// display(n, matrix);
DoubleEven_init();
CreateOddMagicSquare(n);
FullPermutation(0, n / 2);
Combination(n / 2);
displayOddMagicSquare(n);
end_time = clock();
double total_time = (double)(end_time - start_time) / CLOCKS_PER_SEC;
cout << "耗时：" << total_time << endl;
}else if (n % 4 == 0){
start_time = clock();
// CreateDoubleEvenMagicSquare(n);
// display(n, matrix);
DoubleEven_init();
CreateDoubleEvenMagicSquare(n);
FullPermutation(0, n / 2);
Combination(n / 2);
displayDoubleEvenMagicSquare(n);
end_time = clock();
double total_time = (double)(end_time - start_time) / CLOCKS_PER_SEC;
cout << "耗时：" << total_time << endl;
}else{
start_time = clock();
// CreateSingleEvenMagicSquare(n);
// display(n, matrix);
DoubleEven_init();
CreateOddMagicSquare(n / 2);
FullPermutation(0, (n / 2) / 2);
Combination((n / 2) / 2);
displaySingleEvenMagicSquare(n);
end_time = clock();
double total_time = (double)(end_time - start_time) / CLOCKS_PER_SEC;
cout << "耗时：" << total_time << endl;
}
return 0;
}
扫一扫
406
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
