# factorial
function fac(n)
   return (1..n).prod()
end

# raising factorial
function rf(x,n)
   n.prod(|k| x+k)
end

# falling factorial
function ff(x,n)
   n.prod(|k| x-k)
end

# binomial coefficient
function bc(n,k)
   if k<0 or k>n
      return 0
   elif k==0 or k==n
      return 1
   else
      k=min(k,n-k)
      y=1
      for i in 0..k-1
         y = y*(n-i)//(i+1)
      end
      return y
   end
end

# multinomial coefficient
function mc(k)
   return len(k).prod(|i| bc(i.sum(|j| k[j]),k[i]))
end

# finite difference
function Delta(n,f,x,h)
   return (fn|x| f(x+h)-f(x) end^n)(x)
end

begin public pf
   m = {}
   function partition(n,k)
      if n<k
         return 0
      elif k==1 or n==k
         return 1
      else
         t = [n,k]
         if t in m
            return m[t]
         else
            m[t] = partition(n-1,k-1)+partition(n-k,k)
            return m[t]
         end
      end
   end

   # partition function
   function pf(n,k=null)
      if k is null
         return (1..n).sum(|k| partition(n,k))
      else
         return partition(n,k)
      end
   end
end

begin public stirling1
   m = {}
   function s1(n,k)
      if n==0
         return 1 if k==0 else 0
      elif k==0
         return 0
      else
         t = [n,k]
         if t in m
            return m[t]
         else
            m[t] = (n-1)*s1(n-1,k)+s1(n-1,k-1)
            return m[t]
         end
      end
   end

   # Stirling numbers of the first kind
   function stirling1(n,k)
      if n<0 and k<0
         return stirling2(-k,-n)
      else
         return 0 if n<k else s1(n,k)
      end
   end
end

begin public stirling2
   m = {}
   function s2(n,k)
      if n==0
         return 1 if k==0 else 0
      elif k==0
         return 0
      else
         t = [n,k]
         if t in m
            return m[t]
         else
            m[t] = k*s2(n-1,k)+s2(n-1,k-1)
            return m[t]
         end
      end
   end

   # Stirling numbers of the second kind
   function stirling2(n,k)
      if n<0 and k<0
         return stirling1(-k,-n)
      else
         return 0 if n<k else s2(n,k)
      end
   end
end


begin public euler1
   # Eulerian numbers
   m = {}
   function euler1(n,k)
      if n==0
         return 1 if k==0 else 0
      elif k<0
         return 0
      else
         t=[n,k]
         if t in m
            return m[t]
         else
            m[t] = (k+1)*euler1(n-1,k)+(n-k)*euler1(n-1,k-1)
            return m[t]
         end
      end
   end
end

begin public euler2
   # second-order Eulerian numbers
   m = {}
   function euler2(n,k)
      if n==0
         return 1 if k==0 else 0
      elif k<0
         return 0
      else
         t=[n,k]
         if t in m
            return m[t]
         else
            m[t] = (k+1)*euler2(n-1,k)+(2*n-1-k)*euler2(n-1,k-1)
            return m[t]
         end
      end
   end
end

# Bell numbers
function bell(n)
   (0..n).sum(|k| stirling2(n,k))
end

function permutations(a)
   if len(a)<=1
      return [copy(a)]
   end
   x = a[..0]
   b = []
   for p in permutations(a[1..])
      for i in len(a)
         b.push(p[..i-1]+x+p[i..])
      end
   end
   return b
end

begin
   public combinations, partitions

   function comb(m,s)
      if m==0
         return [s[1..0]]
      elif len(s)==0
         return []
      else
         return(comb(m-1,s[1..]).map(|a| s[..0]+a)
            +comb(m,s[1..]))
      end
   end

   function combinations(m,s)
      if s: Map
         return comb(m,list(s)).map(set)
      elif s: String
         return comb(m,s)
      else
         return comb(m,list(s))
      end
   end

   function part(n,k,max)
      if n==0
         return [[]]
      elif n>k*max
         return []
      elif n==k*max
         return [[max]*k]
      elif k==1
         return [[n]]
      else
         a=[]
         for x in 1..min(max,n-1)
            for t in part(n-x,k-1,x)
               a.push(t+[x])
            end
         end
         return a
      end
   end

   function partitions(n,k)
      part(n,k,n)
   end
end