# binomial.py

import math
import decimal
from decimal import Decimal as D


def binomial(n, p, x):
  return math.comb(n, x) * p**x * (1-p)**(n-x)


def binomialD(n, p, x):
  n_choose_x = D(1)
  for i in range(x):
    n_choose_x *= D(n-i)/D(i+1)
  return n_choose_x * (D(p)**x) * ((D(1) - D(p))**(n-x))


def binomialIter(n, p, x):
  q = 1-p
  bin = q**n
  for i in range(1, x):
    bin *= p/q * (n - x + 1)/x
  return bin


def binomialIterD(n, p, x):
  dx = D(x)
  dp = D(p)
  dq = D(1)-dp
  bin = dq**n
  for i in range(1, x):
    bin *= dp/dq * (n - dx + 1)/dx
  return bin


def sumBinomial():
  binSum = 0
  binLogSum = 0
  binNormSum = 0
  for x in range(n):
    binSum += binomial(n, p, x)
    binLogSum += math.exp(binomialLog(n, p, x))
    binNormSum += binomialNorm(n, p, x)



def binomialLog(n, p, x):
  combNumer = sum([ math.log(i) for i in range(x+1, n+1)])
  combDenom = sum([ math.log(i) for i in range(1, n-x+1)])
  pxLog = x * math.log(p)
  negPxLog = (n - x) * math.log(1-p)
  return (combNumer - combDenom) + pxLog + negPxLog


def binomialLogD(n, p, x):
  combNumer = sum([ D(i).ln() for i in range(x+1, n+1)])
  combDenom = sum([ D(i).ln(i) for i in range(1, n-x+1)])
  pxLog = x * D(p).ln()
  negPxLog = (n - x) * (D(1)-D(p)).ln()
  return (combNumer - combDenom) + pxLog + negPxLog


def binomialNorm(n, p, x):
  # normal approximation; accurate for large n
  # https://en.wikipedia.org/wiki/Normal_distribution
  mean = n*p
  std = math.sqrt(n*p*(1-p))
  return phi((x+0.5-mean)/std) - phi((x-0.5-mean)/std)

def phi(x):
  # Cumulative distribution function for the 
  # standard normal distribution
  return (1 + math.erf(x / math.sqrt(2))) / 2



if __name__ == "__main__":
  n = 500
  x = 2
  p = 0.6
  # n, x = map(int, input("n, x=? ").split())
  # p = float(input("p=? "))
  
  print("Binomial Dist:", binomial(n, p, x))
  print("Binomial using logs:", math.exp(binomialLog(n, p, x)))
  print("Binomial using Normal Approx:", binomialNorm(n, p, x))
  
  binSum = 0
  binLogSum = 0
  binNormSum = 0
  for x in range(n):
    binSum += binomial(n, p, x)
    binLogSum += math.exp(binomialLog(n, p, x))
    binNormSum += binomialNorm(n, p, x)
  
  print("\nBinomial Dist Sum:", binSum)
  print("Binomial Sum using logs:", binLogSum)
  print("Binomial Sum using Normal Approx:", binNormSum)