#include <stdio.h>
#include <math.h>
#include <stdlib.h>
/*#include <unistd.h>*/
#include <values.h>
#define delta 1e-6
main()
{
int i,j,k,r,s;
float *p_i=NULL;
float **p_ji=NULL;
float **phi_ij=NULL;
float C,C_pre,validate;
float *sum=NULL;
float p_j;
fscanf(stdin,"%d",&r);
fscanf(stdin,"%d",&s);
p_i=(float *) calloc(r,sizeof(float));
p_ji=(float **)calloc(r,sizeof(float));
for(i=0;i<r;i++)
p_ji[i]=(float *)calloc(s,sizeof(float));
phi_ij=(float **)calloc(r,sizeof(float));
for(i=0;i<r;i++)
phi_ij[i]=(float *)calloc(s,sizeof(float));
for(i=0;i<r;i++)
for(j=0;j<s;j++)
fscanf(stdin,"%f",&p_ji[i][j]);
for(i=0;i<r;i++)
{
validate=0.0;
for(j=0;j<s;j++)
validate+=p_ji[i][j];
if(fabs(validate-1.0)>delta)
{
fprintf(stdout,"\ninvalid input data\n");
}
}
for(i=0;i<r;i++)
{
p_i[i]=1.0;
}
C=-MAXFLOAT;
k=0;
sum=(float *)calloc(r,sizeof(float));
do
{
k++;
for(j=0;j<s;j++)
{
p_j=0.0;
for(i=0;i<r;i++)
p_j+=p_i[i]*p_ji[i][j];
if(fabs(p_j)>=delta)
for(i=0;i<r;i++)
phi_ij[i][j]=p_i[i]*p_ji[i][j]/p_j;
else
for(i=0;i<r;i++)
phi_ij[i][j]=0.0;
}
p_j=0.0;
for(i=0;i<r;i++)
{
sum[i]=0.0;
for(j=0;j<s;j++)
{
if(fabs(phi_ij[i][j])>=delta)
sum[i]+=p_ji[i][j]*log10(phi_ij[i][j])/log10(2.0);
}
sum[i]=pow(2.0,sum[i]);
p_j+=sum[i];
}
for(i=0;i<r;i++)
{
p_i[i]=sum[i]/p_j;
}
C_pre=C;
C=log10(p_j)/log10(2.0);
}
while(fabs(C-C_pre)/C>delta);
free(sum);
sum=NULL;
fprintf(stdout,"the iteration numbre is %d.\n\n",k);
fprintf(stdout,"The capacity of the channel is %.6f bit/symbol.\n\n",C);
fprintf(stdout,"the best input probablity distribution ia :\n");
for(i=0;i<r;i++)
fprintf(stdout,"%.6f",p_i[i]);
fprintf(stdout,"\n");
}
