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

看了网上很多都是用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  x11
 
  x12  x13  x14  x15
  
    x16  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  x11

 x12  x13  x14  x15
 
   x16  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;
}

结果