--     n      m      f(n,m )
sum_ :: Int -> Int -> (Int -> Int -> Int) -> Int
sum_ k j f
    | j <= k = f k j + sum_ k (j+1) f
    | otherwise = 0 
 
-- from http://stackoverflow.com/questions/28896626/tail-recursive-binomial-coefficient-function-in-haskell
binom :: Int -> Int -> Int
binom n 0 = 1
binom 0 k = 0
binom n k = binom (n-1) (k-1) * n `div` k
 
fak :: Int -> Int
fak n
    | n > 1 = n* fak ( n-1)
    | otherwise = 1
 
stirling2 :: Int -> Int -> Int
stirling2 n k = (sum_ k 0 (\ k j -> (pot (-1) (k-j)) * (binom k j) * (pot j n) )) `div` (fak k)
 
pot :: Int -> Int -> Int
pot base exp
    | exp >= 1 = base * pot base (exp-1)
    | otherwise = 1

      PROGRAM main
        integer stirling2;
        write (*,*) stirling2(4,2)
      END PROGRAM
 
      integer RECURSIVE FUNCTION fak ( n ) RESULT ( erg )
        integer n
        if (n <= 1) then
          erg = 1
          return
        endif
 
        erg = n * fak ( n-1 ) 
      END FUNCTION
 
      integer RECURSIVE FUNCTION binom ( n, k ) RESULT ( erg )
        integer :: n,k, ergebnis, i
        if ( k .EQ. 0 ) then
          erg = 1
          return
        else if ( 2*k > n ) then
          erg = binom ( n, n-k )
          return
        endif
 
        ergebnis = n - k + 1
        do i=2, k
          ergebnis = ergebnis * (n - k + i)
          ergebnis = ergebnis / i
        enddo
 
        erg = ergebnis
        return
 
      END FUNCTION
 
      integer FUNCTION stirling2 ( n, k )
        integer :: n,k;
        integer :: summe=0;
        integer :: j
        integer :: fak, binom
 
        do j = 0, k
          summe = summe + ((-1)**(k-j)) * binom ( k,j ) * (j**n)
        enddo
 
        stirling2 = summe / fak (k)
 
      END FUNCTION