reset;

param ver symbolic = "a73b";


param M integer >=1;
param DER;


set P within interval [0,1] = setof{p in 0..0.9 by 0.001} round(p,10);


param minN{p in P} = max(M,if p=0 then 0 else floor(log10(DER/M)/log10(p)));

param maxN_probe{p in P, pr in 1..4} = minN[p]*2^pr;

param Binomial_probe{p in P, pr in 1..4} =
  if M>=2 then
    prod{i in 1..M-1}
    (
      (maxN_probe[p,pr]-M+1+i)/i *
        p^((maxN_probe[p,pr]-M+1)/(M-1))*(1-p)
    )
  else
    p^maxN_probe[p,pr];

param maxN{p in P} = maxN_probe[p, min{pr in 1..4: Binomial_probe[p,pr]<=DER/M} pr];


set NN{p in P} = {minN[p]..maxN[p]};


param Binomial1{p in P, Ntr in NN[p]}=
  if M>=2 then
    prod{i in 1..M-1}((Ntr-M+1+i)/i*p^((Ntr-M+1)/(M-1))*(1-p))
  else
    p^Ntr;

param pbyq{p in P} = p/(1-p);

param Binomial{p in P, Ntr in NN[p], n in Ntr-M+1..Ntr} =
  if n=Ntr-M+1
  then
    Binomial1[p,Ntr]
  else
    Binomial[p,Ntr,n-1]*(Ntr-n+1)/n*pbyq[p];

param Failure{p in P, Ntr in NN[p]} = sum{i in Ntr-M+1..Ntr}Binomial[p,Ntr,i];

param goodN{p in P} = min{Ntr in NN[p]: Failure[p,Ntr]<=DER} Ntr;


param codeword{p in P} =
 if goodN[p]-1 in NN[p] and Failure[p,goodN[p]-1]>DER
   then
     goodN[p]-1 +
       (log(Failure[p,goodN[p]-1]) - log(DER)) /
       (log(Failure[p,goodN[p]-1]) - log(Failure[p,goodN[p]]))
   else goodN[p];



param file symbolic = ver & "-out.csv";


let DER := 1e-5;

param log symbolic = "benchmark.log";
param benchmark_rec symbolic;
param benchmark_start;
param benchmark_stop;

let benchmark_start := _ampl_time;
printf "Version %s:\n",ver;

printf ",FEC\n" > (file);
for{m in 1..10}
{
  let M := m;
  print M, DER;
  for{p in P}
  {
    printf "%.10f,%.20f\n",p,
      codeword[p]/M
      > (file);
  }
  printf "\n" > (file);
}
close (file);

let benchmark_stop := _ampl_time;
let benchmark_rec := sprintf("%s, %s, done in %f seconds",
  sub(ctime(),"^[^ ]* *",""),ver,benchmark_stop-benchmark_start);
print benchmark_rec;
print benchmark_rec >> (log);
close (log);