C语言编程构造拉丁方阵和正交拉丁方阵组

C语言编程构造拉丁方阵和正交拉丁方阵组

将1,2,--,n这n个数填入n*n矩阵，使得每行每列的数两两不同(都是1,2,--,n的全排列)，这样的n阶方阵是拉丁方阵。如果一对n阶拉丁方阵对应的元素构成的有序对两两不同,则称这一对n阶拉丁方阵是正交的。现要编程用计算机构造出所有的n阶拉丁方阵和n阶正交拉丁方阵组，该怎么做？

一个显然的事实是,n阶拉丁方阵中各行i元素（i=1,2---n）是两两不同列的。所以按a1,---,an(为1,2--n的置换)顺序向各行填入ai(i=1,2,--,n)，使其两两不同列,就能得到一个正交拉丁方阵。我们可以按1,2,---,n的顺序填入各数得到一系列正交拉丁方阵,然后从中剔除可以由其中的其他拉丁方阵置换得到的拉丁方阵，经剔除后剩下拉丁方阵姑且叫基拉丁方阵。各基拉丁方阵通过置换{1,2,--,n}->{a1,a2,--,an}就能得到所有n阶拉丁方阵。且每个基拉丁方阵通过置换得到的所有拉丁方阵（包括基拉丁方阵在内）两两不成交,若两个不同的基拉丁方阵相互正交,则对于由这两个基拉丁方阵通过置换{1,2,--,n}->{a1,a2,--,an}衍生出的所有拉丁方阵（包括这两个基拉丁方阵）来说，衍生自其中一个拉丁方正的拉丁方阵和衍生自另一基拉丁方阵的拉丁方阵相互正交，反之亦然。因此只要找出所有的正交基拉丁方阵对就可找出所有的正交拉丁方阵组。

具体实现时，我们需要一个存放找到的拉丁方阵的二维数组Lad(回溯试探操作就对它进行),二维数组position-其第i行第j列的元素表示试探过程中数字j在Lad第i行的位置的列标，我称它为占位矩阵，以及二维数组hold，其第i行第j列的元素为1时表示Lad的第i行第j列已被填充,为0时表示Lad的第i行第j列未被填充。我们用回溯法反复试探，将n个数字全部填入Lad中，使其为拉丁方阵。通过回溯法找到的拉丁方阵及对应的占位矩阵被存放在B1类型的链表中。随后，根据以上分析，我们从找到的的拉丁方阵中剔除可由其中其他拉丁方阵置换得到的拉丁方阵,从而得到存放在B1类型链表中的所有基拉丁方阵,剔除过程中用到了一条性质:若两拉丁方阵对应的占位矩阵的所有列构成的集合相等,则两拉丁方阵中的每一个可以通过置换得到另一个。

得到基拉丁方阵后，就可以轻而易举地通过置换输出所有的n阶拉丁方阵，然后通过判断任意两个基拉丁方阵是否正交来求得并输出（对两个基拉丁方阵作置换）正交拉丁方阵组，问题就解决了。

以下是构造拉丁方阵和正交拉丁方阵组的C语言代码。程序还是有一些问题的,n超过4后等了一段时间后仍然没有结果输出。经过调试发现n=5时回溯操作能顺利完成,但当我越过剔除操作运行至输出操作时发现等了一段时间后程序还是无法运行至输出操作。按步进方式执行剔除操作时没有遇到程序崩溃出错等问题，但无论我按住按键多长时间剔除操作都无法执行完，所以推测是n=5时找出的拉丁方阵过多,程序运行的时间开销太大，短时间内无法执行完毕，没办法。笔者智商捉急，想不到好的方法能解决这个问题，希望有大神能提出改进方案，优化程序，缩短运行时间,目前这个程序至多对n=4能输出正确结果,n>=5时等了一段时间就是不出结果,仔细检查觉得应该不是死循环导致这个问题，暂时也只能这样了，郁闷

#include <stdio.h>
#include <malloc.h>
#include <string.h>
#define N 4
struct B
{
    int **Pa;    //结构体B中Pa指针用于指向找到的拉丁方阵
    int **Pb;    //结构体B中Pb指针用于指向该拉丁方阵对应的占位矩阵
    struct B *next;
};
typedef struct B B1;
void L(int **p1, int *p2, int n);  //求所有全排列的递归函数
void output(B1 *psnew, int fac, int **factor);  //输出基拉丁方阵对应的所有拉丁方阵的函数
int fill(int Lad[][N], int position[][N], int hold[][N], int i, bool TF);  //将i填充入方阵的函数,填充成功返回1，否则返回0
int find(int i, int position[][N], int hold[][N], int k);  //寻找方阵的k+1行可以放置i的位置,找到返回列标，否则返回0
int place(int i, int j, int position[][N], int hold[][N], int k);  //函数，判断方阵的k+1行,j列是否可以放置i，可以返回1，否则返回0
int B(int n);  //求阶乘的递归函数

void main()
{
    int Lad[N][N]={0}, position[N][N]={0}, hold[N][N]={0};  //回溯试探针对Lad方阵进行，试探成功后Lad中保存有找到的拉丁方阵.position为占位矩阵,hold为标志方阵中各位置是否已被填充的矩阵
    int que[N];
    int **factor;  //指向存放所有全排列的方阵
    int fac;
    bool TF;  //为true代表上一个数填充成功，应该开始填充下一个，为false表示当前数填充失败,应该返回上一个数进行回溯
    int i, j, k, flag;  //i在第一个while循环中表示填充的数字,flag表示两占位矩阵列集合是否相等
    B1 *head, *tail, *psnew, *pm;

    head=(B1 *) malloc(sizeof(B1));
    tail=head;
    head->next=NULL;

    for (i=0; i<N; i++)
        que[i]=i+1;

    factor=(int **) malloc((fac=B(N))*sizeof(int *)); 
    for (i=0; i<fac; i++)
        factor[i]=(int *) malloc(N*sizeof(int));  //建立用来保存所有全排列的矩阵

    L(factor, que, N);  //将所有全排列填入上述矩阵

    TF=true;
    i=1;      //TF,i初始化

    while (1)  //回溯试探开始,往Lad方阵中按一定规则填充各数
    {
        if (i==1)  //回溯至或开始填充第一个数
        {
            if (fill(Lad, position, hold, i, TF)==0)
            {
                break;  //第一个数填充失败,所有依次填1,2,--n的拉丁方阵均已找到，循环结束
            }
            i++;  //自增1，准备填充第二个数
            TF=true;  //准备填充第二个数
            continue;  //开始新一轮循环
        }
        else
        {
            if (fill(Lad, position, hold, i, TF)==0)  //第i个数填充失败
            {
                i--;
                TF=false;  //回溯至上一个数
                continue;
            }
            else
            {
                if (i!=N)    //第i个数填充成功,但所有数没有填完
                {
                  i++;
                  TF=true;    //准备填充下一个数
                  continue;
                }
                else    //所有数填充成功，找到一个拉丁方阵
                {
                    psnew=(B1 *) malloc(sizeof(B1));
                    psnew->next=NULL;
                    psnew->Pa=(int **) malloc(N*sizeof(int *));  //建立B1类型节点,令节点的Pa,Pb指向该拉丁方阵和对应的占位矩阵
                    psnew->Pb=(int **) malloc(N*sizeof(int *));

                    for (j=0; j<N; j++)
                    {
                        psnew->Pa[j]=(int *) malloc(N*sizeof(int));
                        psnew->Pb[j]=(int *) malloc(N*sizeof(int));

                        for (k=0; k<N; k++)
                        {
                            *(psnew->Pa[j]+k)=Lad[j][k];              //用找到的拉丁方阵和占位矩阵填充Pa,Pb指向的拉丁方阵和占位矩阵
                            *(psnew->Pb[j]+k)=position[j][k];
                        }
                    }

                    tail->next=psnew;
                    tail=psnew;
                    TF=false;        //回溯
                    continue;
                }
            }
        }
    }

    if (head->next!=NULL)  //从找到的依次填入1,2---n的拉丁方阵中剔除可由其中其他拉丁方阵置换得到的拉丁方阵
    {
        for (psnew=head->next; psnew->next!=NULL; psnew=psnew->next)
        {
            for (tail=psnew->next, pm=psnew; tail!=NULL; tail=tail->next, pm=pm->next)
            {
                flag=0;
               
loop:          for (i=0; i<N; i++)
                {
                    for (j=0; j<N; j++)
                    {
                        for (k=0; k<N; k++)
                        {
                            if (*(psnew->Pb[k]+i)!=*(tail->Pa[k]+j))
                            {
                                flag=1;
                                break;            //判断B1类型链表中两拉丁方阵对应的占位矩阵的所有列构成的两个集合是否相等
                            }
                        }
                        if (flag==1)
                        {
                            flag=0;
                            continue;
                        }
                        flag=1;
                        break;
                    }
                    if (flag==1)
                    {
                        flag=0;
                        continue;
                    }
                    flag=1;
                    break;
                }

                if (flag==0) // flag==1 则两占位矩阵列集合不相等，否则相等
                {
                    for (i=0; i<N; i++)          //两占位矩阵的列集合相等，相应的两拉丁方阵可以通过置换相互得到,故删除tail指向的拉丁方阵，占位矩阵和节点
                    {
                        free(tail->Pa[i]);
                        free(tail->Pb[i]);
                    }

                    free(tail->Pa);
                    free(tail->Pb);
                    pm->next=tail->next;  //让pm的指针域指向tail所指节点的后继节点
                    free(tail);    //删除tail所指节点

                    if (pm->next!=NULL)
                    {
                        tail=pm->next;  //tail指向被删节点的后继节点
                        goto loop;  //立即开始该后继节点和psnew指向的节点的比较
                    }
                    else
                    {
                        tail=pm;
                        if (psnew->next==NULL)
                            goto exit;      //所有比较删除工作均已完成，退出循环
                    }
                }
            }
        }
exit: ;

      if (head->next==NULL)  //没有基拉丁方阵,所以没有拉丁方阵
      {
          printf("不存在N阶拉丁方阵\n");
          printf("不存在N阶正交拉丁方阵组\n");
      }
      else
      {
          printf("所有N阶拉丁方阵为:\n");                           
          for (psnew=head->next; psnew!=NULL; psnew=psnew->next)  //输出所有拉丁方阵
          {
              output(psnew, fac, factor);
          }

          psnew=head->next;
          if (psnew->next==NULL)  //只有一个基拉丁方阵，所以不存在正交拉丁方阵组
              printf("不存在N阶正交拉丁方阵组\n");
          else
          {
              k=0;
              for (psnew=head->next; psnew->next!=NULL; psnew=psnew->next)    //判断正交拉丁方阵组的存在性,并输出正交拉丁方阵组
              {
                  for (tail=psnew->next; tail!=NULL; tail=tail->next)
                  {
                      memset(hold[0], 0, N*N*sizeof(int));  //判断正交前把计数矩阵清零
                      flag=0;

                      for (i=0; i<N; i++)
                      {
                          for (j=0; j<N; j++)
                          {
                              hold[*(psnew->Pa[i]+j)-1][*(tail->Pa[i]+j)-1]++;

                              if (hold[*(psnew->Pa[i]+j)-1][*(tail->Pa[i]+j)-1]>1)  //判断基拉丁方阵是否正交
                              {
                                  flag=1;
                                  break;
                              }
                          }
                          if (flag==1)
                              break;
                      }
                      if (flag==0)    //找到一对正交拉丁方阵组输出
                      {
                          k=1;
                          printf("正交拉丁组:\n");
                          printf("第一组\n");
                          output(psnew, fac, factor);
                          printf("第二组\n");
                          output(tail, fac, factor);
                      }
                  }
              }
              if (k==0)
                  printf("不存在N阶正交拉丁方阵组\n");
          }
      }
    }
}

int fill(int Lad[][N], int position[][N], int hold[][N], int i, bool TF)
{
    int k, j;

    if (TF==true)
        k=0;    //从头开始填充i,行号置1
    else
        k=N-1;  //在当前数字i回溯,行号置尾行号

    while (1)  //对i的回溯试探开始
    {
        if (k==0)  //回溯至或开始填充第一行
        {
            for (j=position[k][i-1]+1; j<=N; j++)  //寻找第一行填充位置
            {
                if (hold[k][j-1]==0)
                  break;
            }

            if (j>N)  //第一行没有填充位置,填充失败，作必要的清理，返回0
            {
                Lad[k][position[k][i-1]-1]=0;
                hold[k][position[k][i-1]-1]=0;
                position[k][i-1]=0;
                return 0;
            }

            if (position[k][i-1]!=0)  //第一行已填充i
            {
                Lad[k][position[k][i-1]-1]=0;
                hold[k][position[k][i-1]-1]=0;  //必要的清理
            }

            position[k][i-1]=j;  Lad[k][position[k][i-1]-1]=i;  hold[k][position[k][i-1]-1]=1;  //填充i
            k++;  continue;  //试探下一行
        }
        else
        {
            if (find(i, position, hold, k)==0)  //在k+1行找不到i的填充位置
            {
                if (position[k][i-1]!=0)    //第k+1行已填充i
                {
                    Lad[k][position[k][i-1]-1]=0;
                    hold[k][position[k][i-1]-1]=0;  //必要的清理
                    position[k][i-1]=0;
                }
                k--;    //回溯至上一行
                continue;
            }
            else
            {
                if (position[k][i-1]!=0)  //第k+1行已填充i
                {
                    Lad[k][position[k][i-1]-1]=0;
                    hold[k][position[k][i-1]-1]=0;  //必要的清理
                }
                position[k][i-1]=find(i, position, hold, k);
                Lad[k][position[k][i-1]-1]=i;            //填充i
                hold[k][position[k][i-1]-1]=1;

                if (k!=N-1)  //还未将i填充完毕
                {
                    k++;      //继续试探下一行
                    continue;
                }
                else       
                {
                    return (1);  //i已填好，填充成功，返回1
                }
            }
        }
    }
}

int find(int i, int position[][N], int hold[][N], int k)
 {
    int j=position[k][i-1];

    for (j++; j<=N; j++)        //搜索k+1行下一个可以填充i的位置
        if (place(i, j, position, hold, k)==1)
            return (j);    //找到返回该位置的列标
    }
    return (0);  //没有下一个可以填充i的位置,返回0
 }

int place(int i, int j, int position[][N], int hold[][N], int k)
{
    int m, flag=0;

    for (m=0; m<k; m++)    //检查k+1行以上各行的i是否都和k+1行j列不在同一列
    {
        if (position[m][i-1]==j)
            flag=1;  //有i在同一列,k+1行j列位置不合法
    }

    if (hold[k][j-1]==1)  //k+1行j列已被填充,位置不合法
        flag=1;

    if (flag==1)
        return (0);  //位置不合法
    else
        return (1);  //位置合法
}

void L(int **p1, int *p2, int n)  //将p2指向的数组中的n个数的所有全排列存放在p1所指向的二维数组中
{
    if (n==1)
        *p1[0]=*p2;
    else
    {

      int **p3;
      int *p4;
      int c, d, f, k, m, i, j;

      d=n-1;
      c=B(d);

      p4=(int *) malloc(d*sizeof(int));

      p3=(int **) malloc(c*sizeof(int *));
      for (i=0; i<c; i++)
          p3[i]=(int *) malloc(d*sizeof(int));

      j=0; f=c;

      for (i=0; i<n; i++)
      {

          m=0;
          for (k=0; k<n; k++)
          {
              if (k==i)
                continue;
              *(p4+m)=*(p2+k);
              m++;
          }

          L(p3, p4, d);

          for (k=0; k<c; k++)
            for (m=0; m<d; m++)
                *(p1[k+j]+m+1)=*(p3[k]+m);

          while (j<f)
          {
            *p1[j]=*(p2+i);
            j++;
          }

          if (i!=d)
            f=f+c;
      }
    }
}

int B(int n)  //计算n的阶乘
{
    if (n==0)
        return (1);
    else
        return (n*B(n-1));
}

void output(B1 *psnew, int fac, int **factor)  //输出psnew指向的节点指向的基拉丁方阵和置换得到的所有其他拉丁方阵
{
    int i, j, k;

    for (i=0; i<N; i++)
    {
      for (j=0; j<N; j++)
          printf("%d ", *(psnew->Pa[i]+j));    //输出基拉丁方阵
          printf("\n");
    }
    printf("\n");

    for (k=1; k<fac; k++)
    {
        for (i=0; i<N; i++)
        {
            for (j=0; j<N; j++)
                printf("%d ", *(factor[k]+(*(psnew->Pa[i]+j)-1)));    //输出置换得到的所有其他拉丁方阵
            printf("\n");
        }
        printf("\n");
    }
}