六角幻方（高斯消元法求解）

看了网上很多都是用dfs解决的，于是自己就写了一篇用高斯消元法的解决方法

问题描述

把 1 2 3 … 19 共19个整数排列成六角形状，如下：

要求每个直线上的数字之和必须相等。共有15条直线哦！

再给点线索吧！我们预先填好了2个数字，第一行的头两个数字是：15 13，参见下图，黄色一行为所求。

请你填写出中间一行的5个数字。数字间用空格分开。

这是一行用空格分开的整数，请通过浏览器提交答案，不要填写任何多余的内容（比如说明性的文字等）

解法

把15条直线看成是15条方程，每一个圆看成一个未知数，且未知数的系数为1即六角幻方对应的未知数：

	 x0  x1  x2 x3  x4  x5  x6 x7  x8  x9  x10  x11x12  x13  x14  x15x16  x17  x18 

x0+x1+x2=38
x3+x4+x5+x6=38
x7+x8+x9+x10+x11=38
x12+x13+x14+x15=38
x16+x17+x18=38
…太多了，我就不写了

代码

/*六角幻方的各个数的下标 x0  x1  x2 x3  x4  x5  x6 x7  x8  x9  x10  x11x12  x13  x14  x15x16  x17  x18 
*/ #include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
using namespace std;
#define N 19
#define M 15
int am[M][N],bm[M],col[N],vi[N],ans[N],ans2[N];
int in[M][7]={ 
{0,1,2,-1},	 //-1作为结尾标志
{3,4,5,6,-1},
{7,8,9,10,11,-1},
{12,13,14,15,-1},
{16,17,18,-1},
{0,3,7,-1},
{1,4,8,12,-1},
{2,5,9,13,16,-1},
{6,10,14,17,-1},
{11,15,18,-1},
{2,6,11,-1},
{1,5,10,15,-1},
{0,4,9,14,18,-1},
{3,8,13,17,-1},
{7,12,16,-1}
};
void init(){memset(am,0,sizeof(am));for(int i=0;i<M;i++){for(int j=0;in[i][j]!=-1;j++){am[i][in[i][j]]=1;}bm[i]=38;/*由于1到19的累加和为109，通过对图中随便构成5条不相交的直线，即  {0,1,2,-1},{3,4,5,6,-1},{7,8,9,10,11,-1},{12,13,14,15,-1},{16,17,18,-1},作为5条直线，且由题意可知每条直线总和相等，所以109/5=38，即每条直线必定和为38*/ } for(int i=0;i<N;i++)col[i]=i;//因为高斯消元中有交换列的处理，所以要标志i列为那个未知数 
}
void r_swap(int i,int j){if(i==j) return ;for(int k=0;k<N;k++){swap(am[i][k],am[j][k]);}swap(bm[i],bm[j]);
}
void c_swap(int k){int i;for(i=k;i<N;i++){if(am[k][i]!=0)break;}if(k==i || i>=N) return ;for(int j=0;j<M;j++){swap(am[j][k],am[j][i]);}swap(col[k],col[i]);
}
int gcd(int a,int b){return b==0?a:gcd(b,a%b);
}
void xiao_yuan(int k,int j){int u=am[k][k],v=am[j][k];if(u==0||v==0) return ;int g=gcd(u,v);for(int i=0;i<N;i++)//必须从0开始，因为这一行要乘以u/g，且从0开始的元素并不一定为0 {am[j][i]=(u*am[j][i]-v*am[k][i])/g;/*原本为am[j][i]=am[j][i]-v/u*am[k][i]);但是v/u可能存在精度误差，所以解决这个问题最好是在计算过程中得出的都是整数将 am[j][i]=am[j][i]-v/u*am[k][i])两边乘以u，得 u*am[j][i]=u*am[j][i]-v*am[k][i]由于将方程组中的某一行两边同时乘以一个数，并不会影响方程组的求解，所以可以令am[j][i]=u*am[j][i]=u*am[j][i]-v*am[k][i]，即 am[j][i]=u*am[j][i]-v*am[k][i]由于 u*am[j][i]-v*am[k][i]可能结果很大，所以除以u和v的最大公约数最终得到 am[j][i]=(u*am[j][i]-v*am[k][i])/g的形式 */ }bm[j]=(u*bm[j]-v*bm[k])/g;//此处跟 am[j][i]=(u*am[j][i]-v*am[k][i])/g同理 
}
void show(){/*static int cnt=0;printf("第%d次消元\n",cnt);*/printf("最终消元结果：\n");for(int i=0;i<N;i++){printf("%5d",col[i]);}printf("   bm\n");for(int i=0;i<(N+1)*5;i++)printf("%c",'-');printf("\n");for(int i=0;i<M;i++){for(int j=0;j<N;j++)printf("%5d",am[i][j]);printf("%5d\n",bm[i]);}printf("\n\n");
}
void guass(){for(int i=0;i<M;i++){//show();for(int j=i;j<M;j++){if(am[j][i]!=0){r_swap(i,j);//找到不为零的行就交换 break;}}c_swap(i);//如果没有找到不为零的行，就找不为零的列，然后交换 for(int j=0;j<M;j++){if(j!=i)xiao_yuan(i,j);}			}
}
void output(){//static int cnt=1;//printf("no.%d\n",cnt++); if(ans2[0]!=15||ans2[1]!=13)//不满足题目要求的结果不显示 return ;printf("%*c%2d  %2d  %2d%*c\n\n",4,' ',ans2[0],ans2[1],ans2[2],6,' ');printf("%*c%2d  %2d  %2d  %2d%*c\n\n",2,' ',ans2[3],ans2[4],ans2[5],ans2[6],3,' ');printf("%2d  %2d  %2d  %2d  %2d\n\n",ans2[7],ans2[8],ans2[9],ans2[10],ans2[11]);printf("%*c%2d  %2d  %2d  %2d%*c\n\n",2,' ',ans2[12],ans2[13],ans2[14],ans2[15],3,' ');printf("%*c%2d  %2d  %2d%*c\n\n\n",4,' ',ans2[16],ans2[17],ans2[18],6,' ');//getchar();
}
void check(){for(int i=11;i>=0;i--){int t=0;for(int j=12;j<N;j++){t+=am[i][j]*ans[j];}t=(bm[i]-t)/am[i][i];if(t<1||t>N)return ;for(int j=i+1;j<N;j++){/*判断填进去的数字是否有重复，此处其实可以使用set树，然后查找，即时间复杂度为log2(N)，但是这里N=19，运算量不大，所以不使用该方法，有兴趣的小伙伴可以自己做一做 */ if(t==ans[j])return ;			}ans[i]=t;}for(int i=0;i<N;i++){ans2[col[i]]=ans[i];//将结果放回第i列对应的未知数内 }output();
}
void f(int *ans,int k){if(k<12)//由于最后求出的行最简形矩阵只有12行，所以另外7未知数要事先填进去 {check();return ;}for(int i=1;i<=N;i++){bool flag=true;for(int j=k+1;j<N;j++){if(i==ans[j]){flag=false;break;}}if(flag){ans[k]=i; f(ans,k-1);}}
}
int main(){init();guass();show();f(ans,18);return 0;
}

结果