81 lines
1.8 KiB
Lua
81 lines
1.8 KiB
Lua
--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 |