--Compute derivatives at the fixed point using Faa Di Bruno's formula.

local ORDER = 100

local binCoefs = {}
--memoized binonmial coefficient
local function Choose( n, k )
  if k < 1 then return 1 end
  if n < 1 then return 0 end
  if k > n then return 0 end
  local idx = n * ORDER + k
  if binCoefs[idx] then return binCoefs[idx] end

  local z = Choose( n - 1, k - 1) + Choose( n - 1, k )
  binCoefs[ idx ] = z
  return z
end

for i = ORDER, 1, -1 do 
  for j = ORDER, 1, -1 do Choose(i, j) end
end

Choose = function( n, k )
  local c =  binCoefs[ n * ORDER + k ] or ( ( k < 1 ) and 1 or 0 )
  return c
end

local function Series( p, order )
  order = order or ORDER
  order = math.min( order, ORDER )
  local d = { [0] = p, p, 1.0 / p } -- d[i] := f^(i) (p)
  local bellPolynomial = {}

  --Generate incomplete Bell polynomials by recursive formula:
  --d[n+1] = - \sum_{k=1}^{n-1} d[k+1] * B(n,k) / p^n
  --B[n][k] = \sum_{i=1}^{n-k+1} Choose( n-1, i-1 ) * d[i] * B( n-i, k-1 )
  --sum_i=1^(n-k+1)

  local function Bell( n, k )
    --Base cases.
    if n == 0 then 
      return ( k == 0 ) and 1 or 0
    elseif k == 0 then 
      return 0 
    end

    --Cached values.
    local idx = n * order + k
    if bellPolynomial[idx] then return bellPolynomial[idx] end

    --Sum.
    local b = 0
    for i = 1, n - k + 1 do
      local a = Choose( n - 1, i - 1 )
      a = a * d[i]
      a = a * Bell( n - i, k - 1 )
      b = b + a
    end
    bellPolynomial[idx] = b
    return b
  end

  local NextDerivative = function()
    local n = #d
    local new = 0
    for k = 1, n - 1 do
      new = new + d[k + 1] * Bell(n, k)
    end
    new = -new / math.pow( p, n )
    d[n + 1] = new
    return new
  end

  for i = 1, order - 2 do NextDerivative() end

  --First few derivatives evaluated at fixed point p.
  
  return d
end

return Series