一般离散无记忆信道容量的逐步迭代算法—C++实现

最新推荐文章于 2022-09-19 10:24:04 发布

大弱智鱼

最新推荐文章于 2022-09-19 10:24:04 发布

阅读量1.8k

点赞数 2

分类专栏： C++ 信计上机作业

本文链接：https://blog.csdn.net/qq_42756072/article/details/89856126

版权

C++ 同时被 2 个专栏收录

19 篇文章 1 订阅

订阅专栏

信计上机作业

15 篇文章 1 订阅

订阅专栏

参考教材：《信息论（第四版）——基础理论与应用》

教材给的算法描述让我怀疑作者有没有亲自编程尝试过。

注意到 x^(lny) = y.

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll N = 1e3;
const double eps = 1e-6;

double p[N];
double p_next[N];
double a[N];
int r;
int s;
double P[N][N];
double C;
double C_new;

int main()
{
     //freopen("input.txt", "r", stdin);
     cin>>r>>s;
     for(int i=1; i<=r; i++)
     {
          p[i] = 1.0 / r;
     }
     for(int i=1; i<=r; i++)
     {
          for(int j=1; j<=s; j++)
          {
               cin>>P[i][j];
          }
     }
     int n=1;
//     for(int i=1; i<=r; i++)
//     {
//          cout<<p[i]<<' ';
//     }
//     cout<<endl;
//     for(int i=1; i<=r; i++)
//     {
//          for(int j=1; j<=s; j++)
//          {
//               //cout<<P[i][j]<<' ';
//               printf("%.8f\t", P[i][j]);
//          }
//          cout<<endl;
//     }

     while(n++)
     {
          for(int i=1; i<=r; i++)
          {
               double temp = 1;
               for(int j=1; j<=s; j++)
               {
                    double temp2 = 0;
                    for(int k=1; k<=r; k++)
                    {
                         temp2 = temp2 + p[k] *P[k][j];
                    }
                    temp2 = P[i][j] / temp2;
                    temp2 = pow(temp2, P[i][j]);
                    temp = temp * temp2;
               }
               a[i] = temp;
          }
          C = 0;
          C_new = 0;
          for(int i=1; i<=r; i++)
          {
               C = C + p[i] * a[i];
               C_new = max(C_new, a[i]);
          }
          C = log( C );
          C_new = log(C_new);
          if( fabs(C - C_new) >eps )
          {
               double sum = 0;
               for(int i=1; i<=r; i++)
               {
                    sum = sum + a[i] * p[i];
               }
               for(int i=1; i<=r; i++)
               {
                    p_next[i] = p[i] * a[i] / sum;
                    p[i] = p_next[i];
               }
          }
          else
          {
               break;
          }
     }
     cout<<n-1<<endl;
     for(int i=1; i<=r; i++)
     {
          printf("%12.8f",p[i]);
     }
     cout<<endl;
     printf("%.6f\n", C_new);

     return 0;
}