构造拉丁方阵和正交拉丁方阵组

将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时就是漫长的等待

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

int main()
{
    int Lad[N][N], position[N][N];  //回溯试探针对Lad方阵进行，试探成功后Lad中保存有找到的拉丁方阵.position为占位矩阵,hold为标志方阵中各位置是否已被填充的矩阵
    for (int i = 0; i < N; ++i)
    {
        for (int j = 0; j < N; ++j)
        {
            Lad[i][j] = N;
            position[i][j] = N;
        }
    }
    bool assist[N][N] = { false };
    int que[N];
    int** factor;  //指向存放所有全排列的方阵
    int fac;
    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;

    int factor_result[N];
    factor_result[0] = 1;
    for (int i = 1; i < N; ++i)
    {
        factor_result[i] = factor_result[i - 1] * (i + 1);
    }

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

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

    i = 0;
    while (true)  //回溯试探开始,往Lad方阵中按一定规则填充各数
    {
        if (i == 0)  //回溯至或开始填充第一个数
        {
            if (fill(Lad, position, i, assist) == false)
            {
                break;   //第一个数填充失败,所有依次填1,2,--n的拉丁方阵均已找到，循环结束
            }
            i++;  //自增1，准备填充第二个数
        }
        else
        {
            if (fill(Lad, position, i, assist) == false)  //第i个数填充失败
            {
                i--;           //回溯至上一个数 
            }
            else
            {
                if (i != N - 1)     //第i个数填充成功,但所有数没有填完
                {
                    i++;     //准备填充下一个数
                }
                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;
                }
            }
        }
    }

    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)
                    {
                        int hold[N][N] = { 0 };
                        flag = 0;

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

                                if (hold[*(psnew->Pa[i] + j)][*(tail->Pa[i] + j)] > 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");
            }
        }
    }
    return 0;
}

bool fill(int Lad[][N], int position[][N], int i, bool assist[][N])
{
    int k;
    if (assist[i][0] == false)
        k = 0;     //从头开始填充i,行号置1
    else
    {
        k = N - 1;   //在当前数字i回溯,行号置尾行号
    }

    while (true)
    {
        int temp = find(i, position, k, assist, Lad);
        if (k == 0 && temp == N || position[i][k] != N)
        {
            Lad[k][position[i][k]] = N;
            assist[i][position[i][k]] = false;
            if (temp == N)
            {
                position[i][k] = N;
                if (k == 0)
                    return false;
            }
        }

        if (temp == N)
            --k;
        else
        {
            position[i][k] = temp;
            Lad[k][temp] = i;             //填充i
            assist[i][temp] = true;
            if (k != N - 1)
                ++k;
            else
                return true;
        }
    }
}

int find(int i, int position[][N], int k, bool assist[][N], int Lad[][N])
{
    int j;
    if (position[i][k] == N)
        j = 0;
    else
        j = position[i][k] + 1;

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

bool place(int i, int j, int k, bool assist[][N], int Lad[][N])
{
    if (assist[i][j])  //检查k+1行以上各行的i是否都和k+1行j列不在同一列
        return false;   //有i在同一列,k+1行j列位置不合法

    if (Lad[k][j] != N)  //k+1行j列已被填充,位置不合法
        return false;    //位置不合法
    return true;  //位置合法
}

void L(int** p1, int row, int col, int* p2, int n, int factor_result[], int currentN)  //将p2指向的数组中的n个数的所有全排列存放在p1所指向的二维数组中
{
    int i;
    if (currentN == 0)
    {
        for (i = 0; i < n; i++)
        {
            if (*(p2 + i) != 0)
                break;
        }
        *(p1[row] + col) = *(p2 + i);
    }
    else
    {
        int c, f, j;

        c = factor_result[currentN - 1];
        j = row; f = row;

        for (i = 0; i < n; i++)
        {
            if (*(p2 + i) == 0)
                continue;

            f = f + c;

            *(p2 + i) = 0;
            L(p1, j, col + 1, p2, n, factor_result, currentN - 1);
            *(p2 + i) = i + 1;

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


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) + 1);    //输出基拉丁方阵
        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)));    //输出置换得到的所有其他拉丁方阵
            printf("\n");
        }
        printf("\n");
    }
}