// ------------------------------------------------------------------------
// File: prob_error_RScode.c
//
// Evaluate the performance of an RS code
// ------------------------------------------------------------------------
// This program is complementary material for the book:
//
// R.H. Morelos-Zaragoza, The Art of Error Correcting Coding, Wiley, 2002.
//
// ISBN 0471 49581 6
//
// This and other programs are available at http://the-art-of-ecc.com
//
// You may use this program for academic and personal purposes only. 
// If this program is used to perform simulations whose results are 
// published in a journal or book, please refer to the book above.
//
// The use of this program in a commercial product requires explicit
// written permission from the author. The author is not responsible or 
// liable for damage or loss that may be caused by the use of this program. 
//
// Copyright (c) 2002. Robert H. Morelos-Zaragoza. All rights reserved.
// ------------------------------------------------------------------------

#include <stdio.h>
#include <math.h>

main(int argc, char *argv[])
{
  int i,j;
  int N,T,m;         /* Length, error correction and bits/symbol */
  double q;          /* Size of Galois field */
  double pb;         /* Input bit error probability */
  double ps;         /* Input symbol error probability */

  double pbo;        /* Output bit error probability */
  double pcd;        /* Probability of correct decoding */
  double picd;       /* Probability of incorrect decoding */

  double comb2(double a, double b);
  double fact(double a);

  if (argc != 5) 
    {
      printf ("\nUsage: error_computation <N> <T> <m> <log10(Input BER)>\n");
      exit(1);
    }
  
  N = strtol(argv[1], NULL, 10);
  T = strtol(argv[2], NULL, 10);
  m = strtol(argv[3], NULL, 10);
  q = pow(2.0,m);
  pb = pow(10.0, (strtol(argv[4], NULL, 10)));

  printf("N = %ld, T = %ld, m = %ld, pb = %e\n", N, T, m, pb);

  /* Symbol error probability */
  ps = 1.0 - pow((1.0-pb), m);
  printf("ps = %e\n", ps);

  /* Output bit error probability */
  pbo = 0.0;
  for (i = T+1; i<=N; i++)
    {
      pbo += ((double)(i+T)/(double)N)*comb2(N,i)*pow(ps,i)*pow(1.0-ps,N-i);
#ifdef DEBUG
      printf("log(factor) = %7.3lf ", log10(i+T) - log10(N) + log10(comb2(N,i))
	     +i*log10(ps) +(N-i)*log10(1.0-ps));
      printf("pbo = %e\n", pbo);
#endif
    }

  pbo *= (pow(2.0,m-1.0)/(pow(2.0,m)-1));
  printf("pbo = %e\n", pbo);

  /* Probability of incorrect decoding */
  pcd = 0.0;
  for (i=0; i<=T; i++)
    {
      /*
      pcd += (comb2(N,i)*pow((q-1.0),i)*pow(ps,i)*pow((1.0-ps),N-i));
      */
      pcd += comb2(N,i)*pow(ps,i)*pow((1.0-ps),N-i);
#ifdef DEBUG
      printf("log(factor) = %7.3lf ", log10(comb2(N,i))+
	     i*log10(ps) +(N-i)*log10(1.0-ps));
      printf("pcd = %e\n", pcd);
#endif
    }
  picd = 1.0 - pcd;
  printf("picd = %e\n", picd);
}

double fact(double a)
{
    double i,tot;

    tot = 1.0;
    for (i=a;i>0.001;i=i-1.0)
    {
        tot = tot * i;
    }
    return(tot);
}

double comb2(double a,double b)
{
    double z,tot;

        if (a<b)
           return(1.0);
    if (b > (a-b))
    {
       tot = 1.0/fact(a-b);
       for (z=a;z>b;z=z-1.0)
       {
          tot = tot * z;
       }
    }
    else
    {
       tot = 1.0/fact(b);
       for (z=a;z>a-b;z=z-1.0)
       {
          tot = tot * z;
       }
    }
    return(tot);
}