Feverish insensate spurt:

- Added basic AST calculator for derivatives.
  - Added direct derivation via Faa Di Bruno formula.
  - Added ansatz method.
  - Factored out polynomial interpolation functionality.

ansatz.lua

@@ -0,0 +1,21 @@
+--Solve the functional equation f' o f = id by ansatz:
+--assume f(x) = ax ^ b. Then
+--f' o f (x) = ba^b * x ^ ( b ( b-1 ) )
+
+local phi = 1 + math.sqrt( 5 ) / 2
+local ihp = 1 - math.sqrt( 5 ) / 2
+local function SeriesCoefficients( )
+  
+  local a, b = math.pow( phi, ihp ), phi
+  local c, d = a, b
+  for i = 1, 100 do
+    
+    --print( c, d )
+    c = c * d
+    d = d - 1
+
+  end
+  
+end
+
+SeriesCoefficients()

diff.lua

@@ -0,0 +1,110 @@
+--Bespoke symbolic differentiation.
+--Little prototypes for building a syntax tree containing only 
+--multiplication, addition, and composition,
+--evaluating its value, and expanding it according to the chain rule.
+
+local op; 
+local function d( x ) return x:diff() end
+local function D( x ) return setmetatable( { ["n"] = x }, op.atom ) end
+
+local computedDerivatives = {}
+
+local __call = D
+local __add = function( a, b ) return setmetatable( { a, b }, op.add ) end
+local __mul = function( a, b ) return setmetatable( { a, b }, op.mul ) end
+local __pow = function( a, b ) return setmetatable( { a, b }, op.compose ) end
+
+op = {
+
+  atom = { 
+    __add = __add,
+    __mul = __mul,
+    __pow = __pow,
+    __call = __call,
+    __tostring = function( self ) return tostring( self.n ) end,
+
+    __index = {
+
+      eval = function( self )
+        return computedDerivatives[ self.n ]
+      end,
+
+      diff = function( self )
+        return D( self.n + 1 )
+      end,
+    }},
+
+  compose = { 
+    __add = __add,
+    __mul = __mul,
+    __pow = __pow,
+    __call = __call,
+    __tostring = function( self ) return "("..tostring( self[1] ).." o "..tostring( self[2] )..")" end,
+
+    __index = {
+      
+      eval = function( self )
+        --All compositions are of the form f^(n) o f,
+        --whose value at a fixed point p is f^(n)(p)
+        return self[1]:eval()
+      end,
+
+      diff = function( self ) 
+        --All compositions are of the form f^(n) o f,
+        --whose derivatives are always ( f^(n+1) o f  )
+        return d( self[1] ) ^ D(0) * D(1) 
+      end,
+    }},
+
+  add = {
+    __add = __add,
+    __mul = __mul,
+    __pow = __pow,
+    __call = __call,
+    __tostring = function( self ) return tostring( self[1] ).." + "..tostring( self[2] ) end,
+
+    __index = {
+      
+      name = "add",
+
+      eval = function( self )
+        return self[1]:eval() + self[2]:eval()
+      end,
+
+      diff = function( self )
+        return d(self[1]) + d(self[2])
+      end,
+    }},
+
+  mul = {
+    __add = __add,
+    __mul = __mul,
+    __pow = __pow,
+    __call = __call,
+    __tostring = function( self ) return tostring( self[1] ).." x "..tostring( self[2] ) end,
+
+    __index = {
+
+      eval = function( self )
+        return self[1]:eval() * self[2]:eval()
+      end,
+
+      diff = function( self )
+        return d(self[1]) * self[2] + self[1] * d(self[2])
+      end,
+    }},
+
+
+}
+
+op.New = function( nodeType, a, b )
+  return setmetatable( { a, b }, op[nodeType] or error( "Node Type Not Recognized" ) )
+end
+
+local a = D(1) ^ D(0)
+for i = 1, 5 do
+  print( a )
+  a = a:diff()
+end
+
+return op

fdb.lua

@@ -0,0 +1,103 @@
+--Compute derivatives at the fixed point using Faa Di Bruno's formula.
+
+local ORDER = 21
+
+local binCoefs = {}
+
+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 = 13, 1, -1 do 
+  for j = 13, 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 )
+  local d = { p, 1.0 / p } -- q[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
+
+  --Function that evaluates Taylor series.
+  return {
+    Interpolant = function( x )
+      x = x - p
+      local pow = 1
+      local fact = 1
+      local y = p
+      for i = 1, ORDER do
+        pow = pow * x
+        fact = fact * i
+        y = y + pow * d[i] / fact
+      end
+      return y
+    end,
+
+    Derivative = function( x )
+      x = x - p
+      local pow = 1
+      local fact = 1
+      local y = 0
+      for i = 1, ORDER do
+        pow = pow * x
+        fact = fact * i
+        y = y + pow * d[i] / fact
+      end
+    end
+  }
+end
+
+return Series

interpolant.lua

@@ -0,0 +1,35 @@
+--Returns a function which takes a sequence of polynomial coefficients and
+--returns two functions which evaluate the polynomial and its derivative, respectively
+return function( coefs )
+  local fixedPoint = coefs[0] or error( "Must have constant coefficient!" )
+  return function(x) 
+
+    x = x - fixedPoint
+    local y = fixedPoint
+    local pow = 1
+    local fact = 1
+    for i = 1, #coefs do
+      pow = pow * x
+      fact = fact / i
+      y = y + pow * fact * coefs[i]
+    end
+    return y
+  end, 
+
+  function(x) 
+
+    x = x - fixedPoint
+    local y = 0
+    local pow = 1
+    local fact = 1
+    for i = 1, #coefs do
+      y = y + pow * fact * coefs[i]
+      pow = pow * x
+      fact = fact / i
+    end
+    return y
+  end
+
+
+
+end

main.lua

@@ -7,41 +7,78 @@
 --get truncated taylor series expansion of f at p
 --plot to get some idea about values, convergence.
 
+local RESOLUTION = 1000
+local EXTENT = 16
 
+local FaaDiBruno = require( 'fdb' )
+local Ansatz = require( 'ansatz' )
+local Poly = require( 'interpolant' )
+
+local a = 0.2
 local function PlotFunction( f )
 
-  local RESOLUTION = 1000
+
 
   local points = {}
   for n = 1, RESOLUTION do
-    points[ 2 * n - 1 ] = n / RESOLUTION
+    points[ 2 * n - 1 ] = EXTENT * n / RESOLUTION
-    points[ 2 * n ] = f( n / RESOLUTION )
+    points[ 2 * n ] = f( EXTENT * n / RESOLUTION )
+
+    --inverse
+    --points[ 2 * RESOLUTION + 2 * n ] = points[ 2 * n - 1]
+    --points[ 2 * RESOLUTION + 2 * n - 1] = points[ 2 * n ]
   end
 
+  local x, y = love.graphics.getDimensions()
   local tf = love.math.newTransform(
-      0, love.graphics.getHeight(),
+    x / EXTENT, y,-- / EXTENT,
-      0, love.graphics.getWidth(),
+    0, x / EXTENT,
-      -love.graphics.getHeight())
+    -y / EXTENT)
-    
+  local dtf = love.math.newTransform( 0, 0, 0, 1, 1 )
-  love.graphics.setColor( 1, 1, 1, 0.5 )
+
-  love.graphics.setLineWidth( 0.003 )
+  --love.graphics.setColor( 1, 1, 1, 0.5 )
+  love.graphics.setLineWidth( 0.1 )
   love.graphics.setLineJoin( "miter" )
   love.graphics.setLineStyle( "smooth" )
   local draw = love.draw or function() end
 
   love.draw = function()
-    draw()
+    draw() 
-    love.graphics.replaceTransform( tf )
+    love.graphics.push("transform")
-    love.graphics.line( points )
+    love.graphics.translate( 0, y )
-  end
+    love.graphics.scale( x, -y )
+    love.graphics.scale( 1 / EXTENT, 1 / EXTENT )
 
-  love.update = function( dt )
+    --love.graphics.scale( 1.0 / EXTENT, 1.0 / EXTENT )
-    print( love.mouse.getPosition( ) )
+    love.graphics.line( points )
+    love.graphics.pop()
+
+    love.graphics.print( a )
   end
 
 end
 
-PlotFunction( function( x ) return x * x end )
+love.update = function( dt )
-PlotFunction( function( x ) return x * x * x end )
+  if love.keyboard.isScancodeDown("w") then
-PlotFunction( math.sin )
+    a = a + 0.05
-PlotFunction( function( x ) return math.exp( x ) - 1 end )
+    love.graphics.setColor( a, 0, 0, a )
+    love.draw = nil
+    --PlotFunction( FaaDiBruno( a ) )
+    --PlotFunction( function(x)return x end)
+  end
+  
+  if love.keyboard.isScancodeDown("s") then
+    a = a - 0.05
+    love.graphics.setColor( 1, 0, 0, 1 )
+    love.draw = nil
+    --[[local f = FaaDiBruno(a)
+    PlotFunction( f.Interpolant )
+    PlotFunction( f.Derivative )]]
+    PlotFunction( function(x)return x end)
+    local f, df = Poly({[0] = 1.5, 1 + a, 2 })
+    PlotFunction( f )
+    PlotFunction( df )
+  end
+  if love.keyboard.isScancodeDown("q") then EXTENT = EXTENT + 0.5 end
+  if love.keyboard.isScancodeDown("e") then EXTENT = EXTENT * 0.99 end
+end