/**
 ** Copyright 1999
 ** Rick Romero
 ** rickr+@cmu.edu
 ** Permission to use and modify is granted,
 ** but not to redistribute original or modified code.
 **/


#include <iostream.h>
#include <fstream.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>


#include "NTL/RR.h"
#include "NTL/fileio.h"

void Usage(char *name) {

  cerr << "Usage: " << name << " -r rate -N numWords -l bitsPerWord [-p precision] [-d decimalPrecision]\n\n";
  cerr << "-r rate       \tRate of 1's in binary words (prob per bit)\n";
  cerr << "-N numWords   \tNumber of total words observed\n";
  cerr << "-l bitsPerWord\tNumber of bits used to compose a word\n";
  cerr << "-p precision  \tHow much precision do you want to use\n";
  cerr << "              \twhile computing the result. Default is\n";
  cerr << "              \t  to use double's which have 53 bits.\n";
  cerr << "-o outFile    \tName of file to store the output in.\n";
  cerr << "              \t  Default is to send to stdout.\n";
  cerr << "-d decimalPrecision\tNumber of decimal digits to use\n";
  cerr << "              \t  when printing the output results.\n";
  exit(-1);

}



RR GetExpectedEntropy(int N, int n, double r) {
  RR e, b, beta, alpha, a;
  ZZ ncr, ncr2;
  int i, k;


  e = 0;

  for (k = 0; k <= n; k++) {
    // Compute prob of a particular word w/ k bits on occurring
    a = pow(r, k) * pow(1 - r, n - k);

    // Compute n choose k
    ncr = 1;
    for (i = k + 1; i <= n; i++)
      ncr = ncr * i;
    ncr2 = 1;
    for (i = 1; i <= n - k; i++)
      ncr2 = ncr2 * i;
    ncr = ncr / ncr2;

    // Compute prob of any word having k bits on
    beta = a * to_RR(ncr);

    // Prob of not having k bits on
    alpha = 1 - a;

    // Prob of getting N words w/out k bits on
    b = 1 - pow(alpha, to_RR(N));

    if (IsZero(b)) {
      // Do nothing, ie assign zero to the entropy
    } else {
      // Entropy due to words w/ k bits on computed as
      //   -beta log(a / b)
      b = -beta * (log(a) - log(b)) / log(2);
      e = e + b;
    }
  }

  return e;

}


int main(int argc, char *argv[]) {
  int N, n, p;
  double r;
  RR e;
  int optChar;
  int outPrecision = 20;


  // Default to standard double precision
  p = 53;

  if (argc < 6) {
    Usage(argv[0]);
  }

  while ((optChar = getopt(argc, argv, "d:o:r:l:N:p:")) != EOF) {
    switch ((char) optChar) {
    case 'd':
      outPrecision = atoi(optarg);
      break;
    case 'r':
      sscanf(optarg, "%lf", &r);
      break;
    case 'l':
      n = atoi(optarg);
      break;
    case 'N':
      N = atoi(optarg);
      break;
    case 'p':
      p = atoi(optarg);
      break;
    default:
      Usage(argv[0]);
      break;
    }
  }


  if ((p < 1) || (n < 1) || (N < 1) || (r <= 0.0)) {
    cerr << "Bad parameter: r = " << r << "   N = " << N << "   n = " << n;
    cerr << "   p = " << r << "\n";
    Usage(argv[0]);
  }

  RR::SetPrecision(p);
  RR::SetOutputPrecision(outPrecision);

  e = GetExpectedEntropy(N, n, r);

  cout << "Expected Entropy:  " << e << "\n";
  cout << "True Entropy:      " << (-to_RR(n) * (r * log(r) + (1 - r) * log(1 - r)) / log(2));
  cout << "\n";

  return 0;

}