最近老师给了一个题目，说是研究一个正常矩阵任意概率置点概率下，双向导通（x,y）的概率（要求有自然边界条件，也就是可以从0->length-1），用代码敲了一下demo，结果发现有个好有趣的结果：不同大小的矩阵，导通概率从某一个概率上升，点越多，上升程度越快（斜率越大），但是都会在（0.592...，35%个点导通）左右相交，符合ziff论文的结果。

　　其实这个结果是统计学的一个结果，不过我们还没学概率论，这个神奇的结论我还不能证明（不过这个证明也之能证明二维的，对于三维以上只能通过模拟得出）。

　　这个结果对于不同的图都有不同的结果，我模拟的这个是最简单的图（自然边界条件矩形图）。

　　下面是demo，没有什么特别的算法，也就是DFS稍微优化一下。

1 #include 
     
        2 #include 
      
         3 #include 
       
          4 #include 
        
           5 #define MARTIX_SIZE 256  6 #define NP 6  7   8 using namespace std;  9  10 typedef struct _tempset 11 { 12     bool flag;//是否被访问过 13     bool dir_x;//规定x一个方向 14     bool dir_y;//规定x一个方向 15     bool is_dot;//是否是一个节点 16 }Set; 17 typedef struct _group 18 { 19     int particle_sum; 20     bool Link_px; 21     bool Link_py; 22 }Group; 23 typedef struct _temp_recur 24 { 25     int x, y, dir; 26 }T_Group; 27 typedef int Position; 28  29 void Inivilize(Set ***const, T_Group **const,Position **const,Position **const); 30 Set **ReFresh_Martix(Set **const, int *const, Position *const, Position *const, int); 31 void Draw_The_Gragh(Set **const); 32 void Destory(Set **const, T_Group *, Position *, Position *); 33 void If_Cls(void); 34 void Show(const int, const int, const int, const double,FILE *); 35 void DFS_Martix(Set **, const int, const int, T_Group *const, Group *, Position *const, Position *const); 36  37 static double poss[NP] = { 0.590, 0.5925, 0.5926, 0.5927, 0.5928, 0.5930}; 38  39 static pair
         
          direction[4]; 40  41 void Inivilize(Set ***const martix_tmp, 42     T_Group **const G_Stack, Position **const Mark_x, Position **const Mark_y) 43 { 44     /*函数名:       Inivilize 45     *调用者函数:    void main(void) 46     *函数形参:      martix_tmp（矩阵指针） 47                     G_Stack（伪递归栈堆指针） 48                     Mark_x，Mark_y（GPA优化的断层标记数组指针） 49     *函数功能:      初始化 50     *返回值：       无 51     */ 52     *martix_tmp = (Set **)malloc(sizeof(Set *)*MARTIX_SIZE); 53     Set *Zone = (Set *)malloc(sizeof(Set)*MARTIX_SIZE*MARTIX_SIZE); 54     for (int i = 0; i < MARTIX_SIZE; i++) 55         (*martix_tmp)[i] = &Zone[i*MARTIX_SIZE]; 56  57     *G_Stack = (T_Group *)malloc(sizeof(T_Group)*MARTIX_SIZE*MARTIX_SIZE);//整形最大值 58     *Mark_x = (Position *)malloc(sizeof(Position)*MARTIX_SIZE); 59     *Mark_y = (Position *)malloc(sizeof(Position)*MARTIX_SIZE); 60  61     //............................................... 62     srand((unsigned)time(NULL)); 63     direction[0].first = -1; direction[0].second = 0; 64     direction[1].first = 1; direction[1].second = 0; 65     direction[2].first = 0; direction[2].second = -1; 66     direction[3].first = 0; direction[3].second = 1; 67 } 68  69 Set **ReFresh_Martix(Set **const Martix, int *const sum_particle, 70     Position *const Mark_x, Position *const Mark_y, int k) 71 { 72     /*函数名:       ReFresh_Martix 73     *调用者函数:    void main(void) 74     *函数形参:      Martix（矩阵指针） 75                     sum_particle（新产生的粒子数） 76                     Mark_x，Mark_y（GPA优化的断层标记数组指针） 77                     k（以第k个概率产生点） 78     *函数功能:      给矩阵产生点，以及赋值 79     *返回值：       Martix（矩阵指针） 80     */ 81     int tmp; 82     memset(Mark_x, 0, sizeof(Position)*MARTIX_SIZE); 83     memset(Mark_y, 0, sizeof(Position)*MARTIX_SIZE); 84     for (int i = 0; i < MARTIX_SIZE; i++) 85     { 86         for (int j = 0; j < MARTIX_SIZE; j++) 87         { 88             tmp = rand(); 89             Martix[i][j].flag = tmp < RAND_MAX * poss[k] ? 1 : 0; 90             Martix[i][j].dir_y = 0; Martix[i][j].dir_x = 0;//强行规定一个方向 91             Martix[i][j].is_dot = 0; 92             if (Martix[i][j].flag) 93             { 94                 (*sum_particle)++; 95                 Mark_x[j]++; Mark_y[i]++; 96                 Martix[i][j].is_dot = 1;  97             } 98         } 99     }100     return Martix;101 }102 103 void Destory(Set **const Martix, T_Group *G_Stack, Position *Mark_x, Position *Mark_y)104 {105     /*函数名:       Destory106     *调用者函数:    void main(void)107     *函数形参:      Martix（矩阵指针）108                     G_Stack（伪递归栈堆指针）109                     Mark_x，Mark_y（GPA优化的断层标记数组指针）110     *函数功能:      销毁111     *返回值：       无112     */113     free(Martix[0]);//去掉集体分配的首地址就可以了114     free(Martix);115     free(G_Stack);//销毁递归栈堆116     free(Mark_x); free(Mark_y);//销毁点集117 }118 119 void If_Cls(void)120 {121     /*函数名:       If_Cls122     *调用者函数:    void main(void)123     *函数形参:      无124     *函数功能:      询问是否清屏125     *返回值：       无126     */127     char choice;128     printf("Do you want to clean the screen?(Y or N)");129     while (1)130     {131         choice = getchar();132         if (choice == 'Y')133         {134             system("cls");135             return;136         }137         else if (choice == 'N')138             return;139         else140         {141             puts("Error Input!Please enter again!\n");142             fflush(stdin);143         }144     }145 }146 147 void Show(const int k, const int link_group_sum, const int loop, const double max_poss, FILE *fp)148 {149     /*函数名:       Show150     *调用者函数:    void main(void)151     *函数形参:      k（以第k个概率产生点）152                     link_group_sum（联通集团个数）153                     loop（循环次数）154                     max（联通最大点个数）155                     p_sum（所有联通集团的粒子总数）156                     fp（指向data.txt的文件指针）157     *函数功能:      显示当前概率集团信息，并且把数据写入data.txt中158     *返回值：       无159     */160     puts("===============================================================================\n");161     printf("\a%cEcho %d:\nPossibility of %.4f\nThe linked graphs` sum is %d\nAccount for %.2f%% in %d graghs\n",162         16, k + 1, poss[k], link_group_sum, 100 * (double)link_group_sum / (double)loop, loop);163     printf("The possibility of(max linked points/sum points)accounts for %.2f%%\n", 100 * max_poss);164     printf("\n");165 166     if (fp != NULL)167     {168         fprintf(fp, "===============================================================================\n");169         fprintf(fp, "%cEcho %d:\nPossibility of %.4f\nThe linked graphs` sum is %d\nAccount for %.2f%% in %d graghs\n",170             7, k + 1, poss[k], link_group_sum, 100 * (double)link_group_sum / (double)loop, loop);171         fprintf(fp, "The possibility of(max linked points/sum points)accounts for %.2f%%\n", 100 * max_poss);172         fprintf(fp, "\n");173     }174 }
         
        
       
      
     

#include "plug.h"static bool x_full, y_full;int main(void)//森林火灾项目的计算连通性的概率{    int loop, sum_particle, link_group_sum, max_particle_sum;    double max_poss, poss_tmp;            Set **martix = NULL;                                        //矩阵      Group tmp;                                                        T_Group *G_Stack;                                            //范围太大，不用递归，手动直接开大栈提高效率    Position *Mark_x = NULL, *Mark_y = NULL;                    //Gpa优化标记数组    FILE *fp = NULL;                                            //文件指针    //......................................................................................    Inivilize(&martix, &G_Stack, &Mark_x, &Mark_y);                //初始化    while (1)    {        printf("Please enter the times of the loops(The martix`s size is %d*%d)\n(The program will stop until you input zero)\n",             MARTIX_SIZE, MARTIX_SIZE);        while (1)        {            scanf("%d", &loop);            if (loop >= 0)break;            puts("DO NOT try to input a negative loop");            system("cls");            puts("Please enter the time of the loop(break until you input 0)");        }        if (loop == 0) break;        fflush(stdin);        puts("The Date will be written in ""data.txt"" in the folder of this program\n");        if ((fp = fopen("data.txt", "r")) == NULL)            puts("""data.txt"" doesn`t exsit,the program will create a new one\n");        fp = fopen("data.txt", "a+");        time_t c1 = clock();                                    //这里先得到一个循环开始的时间        for (int k = 0; k < NP; k++)        {            link_group_sum = 0; max_poss = 0;            /*DFS: 递归每一个概率，和每一张图，然后每一个点深度优先搜索看是否联通             *GPA优化：如果Mark数组出现断层，则此图一定不连通(x_full和y_full)            */            for (int i = 0; i < loop; i++)            {                sum_particle = 0; x_full = 1; y_full = 1;                 martix = ReFresh_Martix(martix, &sum_particle, Mark_x, Mark_y, k);                max_particle_sum = 0;                for (int y = 0; y < MARTIX_SIZE && x_full && y_full; y++)                {                    for (int x = 0; x < MARTIX_SIZE && x_full && y_full; x++)                    {                        if (martix[y][x].flag)                        {                            DFS_Martix(martix, x, y, G_Stack, &tmp, Mark_x, Mark_y);                            if (tmp.Link_px && tmp.Link_py)                            {                                link_group_sum++;                                if (tmp.particle_sum > max_particle_sum)                                    max_particle_sum = tmp.particle_sum;                            }                        }                    }                }                poss_tmp = (double)max_particle_sum / (double)sum_particle;                max_poss = poss_tmp > max_poss ? poss_tmp : max_poss;            }            Show(k, link_group_sum, loop, max_poss, fp);        //联通的概率（loop个图），联通的时候最大联通集团粒子数        }        time_t c2 = clock();                                    //得到循环结束的时间（ms）        printf("The program has run %lldms\n\n", c2 - c1);                fclose(fp);                                                //记得关闭文件流        system("pause");         If_Cls();     }     Destory(martix, G_Stack, Mark_x, Mark_y);                    //销毁栈组        fflush(stdin);    system("cls"); system("pause");    return 0;}void DFS_Martix(Set **martix, const int st_x, const int st_y,     T_Group *const G_Stack, Group *Part_Group, Position *const Mark_x, Position *const Mark_y){    /*函数名:        DFS_Martix     *调用者函数:    void main(void)     *函数形参:      martix（矩阵）                     st_x，st_y（开始坐标横纵坐标）                     G_Stack（伪递归栈堆）                     Part_Group（集团参数指针）                     Mark_x，Mark_y（GPA优化的断层标记数组）     *函数功能:      DFS搜索联通集团     *返回值：       无    */    bool If_Push,                                //出入栈标志        if_link_y = 0, if_link_x = 0;            //图是否联通标志    bool round_dir_x = 0, round_dir_y = 0,        //环图指向标记：规定沿着某个轴沿一个方向突然发生反转时，某个方向就会被标记1        tmp_dir_x, tmp_dir_y;                        int dir_tmp,                                        particle_sum = 0,                        //粒子数        top = 0;                                //栈顶位    Position dir_x, dir_y;    T_Group tmp;    //......................................................................................        //初始化栈..........    martix[st_y][st_x].flag = 0; martix[st_y][st_x].dir_x = 0; martix[st_y][st_x].dir_y = 0;    Mark_x[st_x]--; Mark_y[st_y]--;    tmp.x = st_x, tmp.y = st_y; tmp.dir = dir_tmp = 0;    G_Stack[top++] = tmp; particle_sum++;    //..................    Part_Group->Link_px = 0;    Part_Group->Link_py = 0;    while (top != 0)    {        If_Push = 0;        tmp = G_Stack[top - 1];                    //读取栈顶元素,以及栈顶元素环图两个指向标记        round_dir_x = martix[tmp.y][tmp.x].dir_x;        round_dir_y = martix[tmp.y][tmp.x].dir_y;        //每弹出一个点，就把该节点的x和y做标记,直接从上一个搜索方向之后开始搜索        for (int i = dir_tmp; i < 4; i++)        {            dir_x = tmp.x + direction[i].first;            tmp_dir_x = (dir_x < 0 || dir_x == MARTIX_SIZE) ? !round_dir_x : round_dir_x;            dir_x = (dir_x + MARTIX_SIZE) % MARTIX_SIZE; //循环坐标            dir_y = tmp.y + direction[i].second;            tmp_dir_y = (dir_y < 0 || dir_y == MARTIX_SIZE) ? !round_dir_y : round_dir_y;            dir_y = (dir_y + MARTIX_SIZE) % MARTIX_SIZE; //循环坐标            if (martix[dir_y][dir_x].is_dot && !martix[dir_y][dir_x].flag)            {                //如果环图指向标记出现相反，则图已经联通                if (!if_link_x)                    if_link_x = martix[dir_y][dir_x].dir_x == tmp_dir_x ? 0 : 1;                if (!if_link_y)                    if_link_y = martix[dir_y][dir_x].dir_y == tmp_dir_y ? 0 : 1;            }            if (martix[dir_y][dir_x].flag)            {                If_Push = 1;                Mark_x[dir_x]--; Mark_y[dir_y]--;                round_dir_x = tmp_dir_x; round_dir_y = tmp_dir_y;                martix[dir_y][dir_x].flag = 0;                martix[dir_y][dir_x].dir_x = round_dir_x; martix[dir_y][dir_x].dir_y = round_dir_y;                if (!Mark_x[dir_x]) x_full = 0;                if (!Mark_y[dir_y]) y_full = 0;                tmp.x = dir_x; tmp.y = dir_y; tmp.dir = i;                G_Stack[top++] = tmp; particle_sum++;                break;            }        }        dir_tmp = 0;                            //这里一定要初始化为0，否则接下来的点就无法搜索其他方向        if (!If_Push)                            //如果没有可以入栈的，那说明开始递归了，pop栈顶        {            top--;            dir_tmp = tmp.dir;                    //定义为当前点上一次出去的点的方向        }    }    if (particle_sum >= MARTIX_SIZE)    {        Part_Group->Link_px = if_link_x;        Part_Group->Link_py = if_link_y;    }    Part_Group->particle_sum = particle_sum;}