用c编写求信息论中各种熵和信息量的程序

#include<stdio.h>
#include<math.h>
#define u 20
int i,j,n,m;
float Ix[u],I[u][u],H_X,H_Y,H_XY,I_XY,H_XpY,Pypx[u][u],Px[u],H_YpX,Py[u],Pxpy[u][u],Pxy[u][u];
/*定义的Ix[]=I(x[i])自信息；I[][]=I(x[i];y[j])互信息；H_X=H(X)平均自信息；H_XY=H(XY)联合熵;
I_XY=I(X;Y)平均互信息；H_XpY=H(X|Y)、H_YpX=H(Y|X)条件熵;Pypx[i][j]=P(y[j]|x[i])条件概率;
Px[i]=P(x[i])发x[i]的概率；H_XpY=H(Y/X)条件熵;Py[j]=P(y[j])收到y[j]的概率；Pxpy[i][j]=P(x[i]/y[j])条件概率；
Pxy[i][j]=P(x[i]y[j])联合概率*/
/*定义以2为底的对数函数*/
float log2(float x)
{
 float z;
 z=(float)(log(x)/log(2));
 return z;
}
/*求自信息的函数*/
float self_inf(float x)
{
  float z;
  z=-log2(x);
  return z;
}
/*求互信息的函数*/
float mut_inf(int x,int y)
{
 float z;
 z=log2(Pxpy[x][y]/Px[x]);
 return z;
}
/*求平均自信息的函数*/
float entropy(float *x,int n)
{
  float z=0;
  for(i=1;i<=n;i++)
  {
   z+=(*(x+i))*log2(1/(*(x+i)));
  }
  return z;
}
/*求联合熵的函数*/
float joint_entropy(float (*p)[u])
{
 float z=0;
 for(i=1;i<=n;i++)
  for(j=1;j<=m;j++)
   {
    z+=*(p[i]+j)*log2(1/(*(p[i]+j)));
   }
 return z; 
}
main()
{
 float s=0;
 printf("/nplease in put the dimension of 'X' and 'Y'/n");
 scanf("%d %d",&n,&m);
 printf("/nThe dimension of X is n=%d/nThe dimension of Y is m=%d/nPlease input the condition probability:P(y[j]/x[i]),",n,m);
 printf("(after you input one number please click the 'enter')/n");
/*条件概率P(y[j]/x[i])赋值*/
 for(i=1;i<=n;i++)
  {for(j=1;j<=m;j++)
   {printf("P(y[%d]/x[%d])=",j,i);
          scanf("%f",&Pypx[i][j]);
         }
  }
 printf("please input P(x[i]):(after you input one number please click the 'enter')/n");
    for(i=1;i<=n;i++)
 {printf("P(x[%d])=",i);
  scanf("%f",&Px[i]);
 }
/*判断输入X的概率是否正确，不正确则退出程序*/
 for(i=1;i<=n;i++)
 {if((Px[i]<0)||(Px[i]>1))
   {printf("Please input right value of probability/n");
      goto End_exe;
        }
  s+=Px[i];
 }
 if((s-1)*(s-1)>0.001)
    {printf("Please input the right value of probability/n");
 goto End_exe;
 }
 for(i=1;i<=n;i++)
 {s=0;
  for(j=1;j<=m;j++)
  {s+=Pypx[i][j];
  }
  if((s-1)*(s-1)>0.001)
  {printf("Please input the right value of probability/n");
   goto End_exe;
  }
 }
/*计算Py[j]=P(y[j])收到y[j]的概率；Pxpy[i][j]=P(x[i]/y[j])条件概率；Pxy[i][j]=P(x[i]y[j])联合概率*/
 for(i=1;i<=n;i++)
 {for(j=1;j<=m;j++)
  {Pxy[i][j]=Pypx[i][j]*Px[i];
   printf("Px[%d]y[%d]=%f/n",i,j,Pxy[i][j]);
     }
    }
    for(j=1;j<=m;j++)
    {
     s=0;
     for(i=1;i<=n;i++)
     s+=Pxy[i][j];
     Py[j]=s;
     printf("P(y[%d])=%f/n",j,Py[j]);
    }
    for(i=1;i<=n;i++)
    {for(j=1;j<=m;j++)
     {Pxpy[i][j]=Pxy[i][j]/Py[j];
      printf("P(x[%d]/y[%d])=%f/n",i,j,Pxpy[i][j]);
     }
    }
/*结束 计算Py[j]=P(y[j])收到y[j]的概率；Pxpy[i][j]=P(x[i]/y[j])条件概率；Pxy[i][j]=P(x[i]y[j])联合概率*/
/*输出自信息*/
    printf("********************************************************************************/nThe self-information :/n");
    for(i=1;i<=n;i++)
 {
  Ix[i]=self_inf(Px[i]);
   printf("Ix[%d]=%f/n",i,Ix[i]);
 }
/*输出互信息*/
    printf("/nThe mutal-information :/n");
    for(i=1;i<=n;i++)
 {
     for(j=1;j<=m;j++)
  {
      I[i][j]=mut_inf(i,j);
   printf("I(x[%d],y[%d])=%f/n",i,j,I[i][j]);
     }
    }
/*输出平均自信息（熵）*/
    H_X=entropy(Px,n);
    printf("/nThe entropy of X :H(X)=%f/n",H_X);
    H_Y=entropy(Py,m);
    printf("/nThe entropy of Y :H(Y)=%f/n",H_Y);
/*输出联合熵*/
    H_XY=joint_entropy(Pxy);
    printf("/nThe joint entropy of X and Y :H(XY)=%f/n",H_XY);
/*输出平均互信息*/
    I_XY=H_X+H_Y-H_XY;
    printf("/nThe average mutual information of X and Y :I(X;Y)=%f/n",I_XY);
/*输出条件熵*/
    H_XpY=H_X-I_XY;
    H_YpX=H_Y-I_XY;
    printf("/nThe conditional entropy:/nH(X/Y)=%f/t/tH(Y/X)=%f/n",H_XpY,H_YpX);
End_exe:;
}