P <- 0
if ((k>0) & (k<n)) 
{ 
  logBinomCoeff <- sum(log(seq(n-k+1, n))) - sum(log(seq(1, k))) 
  logP <- logBinomCoeff + k*log(p) + (n-k)*log(1-p)
  P <- exp(logP)
}else if(k == 0)
{
  P <- (1-p)^n
}else if(k == n)
{
  P <- p^n
}

P <- pbinom(k,n,p)-pbinom(k-1,n,p)