-- brit2dice.hs
-- Haskell program for calculating proabilities of Britannia battles

-- ak = number of values (out of 6) that the attacker kills on (normally 2)
-- dk = number of values (out of 6) that the defender kills on (normally 2)
-- aa = original number of attackers
-- dd = original number of defenders
-- a = final number of attackers
-- d = final number of defenders

p0(ak,dk,aa,dd,a,d) = p(a,d) where

  p(a,d) | (a,d)==(aa,dd) = 1
  p(a,d) = sum [(pp!!i!!j)*t(i,j,a,d) | i<-[a..aa], j <-[d..dd], (i,j)/=(a,d)]

  pp = [[p(a,d)|d<-[0..dd]]|a<-[0..aa]]

  t(i,j,a,d) | i<j-d || j<i-a = 0
  t(i,j,a,d) = qa(i,j,d) * qd(j,i,a) / (6^(i+j) - (6-ak)^i * (6-dk)^j)

  qa(i,j,0) = sum [ak^m * (6-ak)^(i-m) * k(i,m) | m<-[j..i]]
  qa(i,j,d) = ak^(j-d) * (6-ak)^(i-j+d) * k(i,j-d)

  qd(j,i,0) = sum [dk^m * (6-dk)^(j-m) * k(j,m) | m<-[i..j]]
  qd(j,i,a) = dk^(i-a) * (6-dk)^(j-i+a) * k(j,i-a)

  k(n,m) = fromInt(product [n-m+1..n] `div` product [1..m])