Big commit, basic functionality + plotting.

backend.lua

@@ -0,0 +1,5 @@
+--Dependencies for displaying things.
+local backend = {}
+
+
+return backend

conf.lua

@@ -0,0 +1,3 @@
+function love.conf( t )
+  t.console = true
+end

enumeration.lua

@@ -0,0 +1,7 @@
+--Lazy enumerator of countable infinite set.
+--Pass a function closure for generating the set element-by-element, in order.
+--This function can return a "simple" value followed by an "exact" value
+--Returns a table including:
+-- an array with all elements generated so far (grows as needed upon access)
+local enum = { __index = function( t, n ) for i = #t + 1, n do t[i] = t.next() end return rawget(t, n) end }
+return function(nextElement) return setmetatable( { ["next"] = nextElement }, enum ) end

enumerations.lua

@@ -0,0 +1,88 @@
+local enum = require 'enumeration'
+
+local function CalkinWilf()
+  local n, d = 1, 1
+  return enum( function()
+    repeat n, d = d, ( 2 * d * math.floor( n / d ) - n + d )
+    until n < d
+    return n / d
+  end )
+end
+
+local function Random()
+  return enum( math.random )
+end
+
+--Breadth first traversal of the open left subtree of the Stern-Brocot tree.
+--Formed by inserting the mediants of previous rationals in the tree.
+local function SternBrocot()
+
+  --Rational numbers.
+  local R
+  do
+    local rationals = {}
+    local rmt = {
+      __lt = function(a, b) return a[1]/a[2] < b[1]/b[2] end,
+      __eq = function(a, b) return a[1] == b[1] and a[2] == b[2] end,
+      __pow = function(a, b) return R( a[1] + b[1], a[2] + b[2] ) end, --mediant
+      __call = function(t) return t[1]/t[2] end,
+      __tostring = function(t) return string.format( "%d / %d", t[1], t[2] ) end }
+
+    R = function( n, d )
+      if not rationals[d] then rationals[d] = {} end
+      local den = rationals[d]
+
+      if den[n] then return den[n] end
+
+      local r = setmetatable( { n , d }, rmt )
+      den[n] = r
+      return r
+    end
+  end
+
+  --Flat list of all entries traversed so far.
+  local seq = { R(0, 1), R(1, 2), R(1, 1) }
+  local idx = 0
+  local len = 3
+
+  local function nextRational()
+    --Return every other entry: these are the most recently added mediants.
+    idx = idx + 2
+    if idx < len then return seq[idx] end
+
+    --We have finished iterating over the current level of the tree.
+    --Create a new level by adding all the new mediants.
+    
+    local nSeq = {}
+    for i = 1, len - 1 do
+      nSeq[ 2 * i - 1 ] = seq[i]
+      nSeq[ 2 * i ] = seq[i] ^ seq[i + 1]
+    end
+    nSeq[ 2 * len - 1 ] = seq[ len ]
+    seq = nSeq
+    len = #seq
+    --print( "RATIONALS: NEW LEVEL", len, unpack( nSeq ) )
+    idx = 0
+    return nextRational()
+  end
+
+  return enum( nextRational )
+end
+
+--Enumeration of dyadic rationals in the open unit interval, ascending lexicographic.
+local function Dyad()
+  local n, d = -1, 2
+  return enum( function()
+      n = n + 2
+      if n > d then n, d = 1, d * 2 end
+      return n / d
+    end )
+end
+
+return { 
+  Dyad = Dyad, 
+  SternBrocot = SternBrocot,
+  CalkinWilf = CalkinWilf,
+  Random = Random,
+  QuadIrr = false,
+}

main.lua

@@ -0,0 +1,92 @@
+local enumerations = require 'enumerations'
+local twintree = require 'twintree'
+local step = coroutine.wrap( twintree.buildIncremental )
+local a = enumerations.Dyad()
+local b = enumerations.Random()
+step( a, b )
+
+--Use LOVE for plotting.
+local love = assert( love )
+
+local canvas = love.graphics.newCanvas()
+local numMapped, pointList = 0
+
+local function swap( i, j )
+  b[i], b[j] = b[j], b[i]
+end
+
+local function reset()
+  print( "changing coroutine" )
+  step = coroutine.wrap( twintree.buildIncremental )
+  print( "swapping" )
+  for i = 1, 10000 do
+    swap( i, math.random( 1, 10000 ) )
+  end
+  
+  pointList = nil
+  numMapped = 0
+  print( "step" )
+  step( a, b )
+end
+
+local function new()
+  pointList = step()
+  numMapped = numMapped + 1
+end
+
+local function paint()
+  local type = type
+  --Make sure the whole point list consists of numbers!
+  --We use the __call metamethod in non-number elements to coerce to numbers.
+  --TODO: make this more generic, maybe amenable to FFI types.
+  for i = 1, #pointList do
+    if type( pointList[i] ) ~= 'number' then pointList[i] = pointList[i]() end 
+  end
+  
+  --print( "Points: " )
+  --for i = 1, #pointList - 1, 2 do
+  --  print(pointList[i], pointList[i + 1] )
+  --end
+  
+  love.graphics.setCanvas( canvas )
+  love.graphics.clear()
+  love.graphics.setColor( 1, 1, 1, 0.5 )
+  love.graphics.print( numMapped )
+  love.graphics.push()
+  love.graphics.scale( love.graphics.getWidth(), love.graphics.getHeight() )
+  love.graphics.translate( 0, 1 )
+  love.graphics.scale( 1, -1 )
+  love.graphics.points( pointList )
+  love.graphics.line( pointList )
+  love.graphics.pop()
+  love.graphics.setCanvas()
+end
+
+function love.wheelmoved()
+  for i = 1, 3 do new() end
+  paint()
+ --[[ for i = 1, #sb do
+    print( i, sb[i], dy[i] )
+  end]]
+end
+
+function love.mousepressed()
+  --[[for i = 1, #sb do
+    print( i, sb[i] )
+  end]]
+end
+
+function love.draw()
+  love.graphics.setCanvas()
+  love.graphics.setColor( 1, 1, 1, 1 )
+  love.graphics.draw( canvas )
+end
+
+function love.keypressed()
+  print( "pressed key!" )
+  return reset() and love.wheelmoved()
+end
+
+function love.load()
+  love.graphics.setLineWidth( 0.001 )
+end

rapport.lua

@@ -0,0 +1,132 @@
+local function Gap(t, i)
+  if not t[i] then return i else return Gap(t, i + 1) end
+end
+
+local function Gaps(t, i)
+  return Gap, t, i
+end
+
+--Data structure supporting fast ordered insertion and search.
+local function Order()
+  local o = {}
+  
+  local list = { val = nil, prev = nil, next = nil }
+  
+  o.Get = function( x )
+    
+  end
+  
+  o.GetAdjacent = function( node )
+    return node.prev.val, node.next.val
+  end
+  
+  o.Insert = function( x )
+    if not list.val then
+      list.val = x
+      return list
+    end
+    
+    local node = list
+    while node and node.val < x do node = node.next end
+    
+    
+  end
+  
+  return o
+end
+
+--Construct an isomorphism between two countable dense total orders
+--these orders are specified by functions which return the next element in some enumeration
+--these elements must be comparable with comparison operators <, >, &c.
+local function Isomorphism( Enumeration, InverseEnumeration, inv )
+
+  --Exported table.
+  local t = {}
+  local enum, order, mapping, inverse
+
+  --Another Isomorphism structure,
+  --which contains the image of this isomorphism and the inverse
+  --isomorphism mapping back to our dense countable linear order.
+  inverse = inv or Isomorphism( InverseEnumeration, nil, t )
+
+  --Array for caching and lazily evaluating the supplied Enumeration,
+  --an iterator that enumerates the elements of a dense countable linear order.
+  enum = { Enumeration() }
+  setmetatable(enum, {__index = function(_, i)
+        for j = #enum + 1, i do
+          enum[j] = Enumeration()
+        end
+        return rawget( enum, i )
+      end})
+  t.enum = enum --Expose as field 
+
+  --Structure supporting fast ordered insertion and search.
+  --Values are i such that mapping[order[i]] is not nil,
+  --sorted according to getmetatable( enum[i] ).__lt
+  order = { 1 }
+
+  --Find the indices i and k such that
+  --enum[i] is the greatest lower bound on enum[j] in the mapping's support
+  --enum[k] is the least upper bound on enum[j] in the mapping's support
+  local function GetBounds( leaf )
+
+  end
+
+  --Array defining isomorphism mapping this dense countable order to the one in inverse.
+  --If mapping[i] = j, then f(enum[i]) = inverse.enum[j]
+  mapping = { 1 }
+  setmetatable( mapping, {__newindex = function(map, k, v)
+    rawset(mapping, k, v)
+    
+    
+      end})
+
+  --Iterates over the natural numbers outside the isomorphism's support,
+  --so returns the least n such that n >= i and mapping[n] == nil.
+  --(Warning, infinite loop! Make sure to break out of this.)
+  local function NextGap(i)
+    if not mapping[i] then return i else return NextGap(i + 1) end --tail-recursive
+  end
+  
+  --Include another point in the mapping's support.
+  --Infinite loop: ))<<>>((
+  function t.ExpandMap()
+    local gapIdx = NextGap(1)
+    local pred, succ = GetBounds( OrderedInsertion( gapIdx ) )
+    mapping[gapIdx] = inverse.AddNextBetween( i, mapping[pred], mapping[succ] )
+    return inverse.ExpandMap()
+  end
+
+  --Find the first j such that mapping[j] == nil,
+  --enum[lower] < enum[j] < enum[upper]
+  --Then map j to i, and add j to the order.
+  function t.AddNextBetween( i, lower, upper )
+    local p, q, r, j = enum[lower], nil, enum[upper], nil
+
+    for k in Gaps( 1 ) do
+      q = enum[k]
+      if  ((not p) or (p < q))
+      and ((not r) or (q < r)) then
+        j = k
+        break
+      end
+    end
+
+    mapping[j] = i
+    OrderedInsertion(j)
+    return j
+  end
+
+  --Get the isomorphism as a dictionary, for plotting.
+  function t.GetMapping( )
+    local map = {}
+    for idx, imageIdx in pairs( mapping ) do
+      map[ enum[ idx ] ] = inverse.enum[ imageIdx ]
+    end
+    return map
+  end
+
+  return t
+end
+
+return setmetatable( { Isomorphism = Isomorphism }, {__call = Isomorphism } )

sequence.lua

@@ -0,0 +1,25 @@
+--Partial sequence s:S->N, S <= N, N = { 1, 2, 3, ... }
+--Consists of:
+  --A sequence stored as a (sparse) array
+  --["cmin"], a key whose value is the smallest element outside the support of the sequence
+local sequence = {}
+sequence.__index = sequence
+function sequence.new( name )
+  return setmetatable( { cmin = 1, name = name }, sequence )
+end
+
+function sequence:insert( x, y )
+  if self.cmin == x then
+    repeat self.cmin = self.cmin + 1 until not self[self.cmin]
+  end
+  self[x] = y
+  --debug:
+  --print( "First unclaimed element:", self.name, self.cmin )
+end
+
+--Maybe leave this for later, we'll inline it for now.
+function sequence:complement()
+  
+end
+
+return sequence

twintree.lua

@@ -0,0 +1,103 @@
+local sequence = require 'sequence'
+
+local maxLevels = 20 --Maximum depth of tree.
+local baseIdx = 1 --Base index of implicit tree. Is it a 0-based or 1-based array?
+local tree = {}
+
+local function insert( aTree, aSeq, aVals, bTree, bSeq, bVals )
+
+  --print( "INSERTING" )
+  local treeIdx = 1 --Pointer into implicit binary tree.
+  local aIdx, bIdx = aSeq.cmin, bSeq.cmin --First elements of a and b outside the support of the sequence
+  local aVal, bVal = aVals[aIdx], bVals[bIdx]
+  local aTest, bTest
+
+  --Search for appropriate bIdx.
+  --Sentinel value.
+  local newCandidate = false
+  while aTree[treeIdx] do
+    if newCandidate then
+      newCandidate = false
+      treeIdx = 1
+      repeat bIdx = bIdx + 1 until ( not bSeq[bIdx] )
+      bVal = bVals[bIdx]
+    end
+
+    --invariants guarantee these aren't empty:
+    --loop invariant, that aTree[treeIdx] is nonempty
+    --struct invariant, that bTree and aTree have identical support
+    aTest, bTest = aTree[treeIdx], bTree[treeIdx] 
+    --Go left if bounded above, go right if bounded below,
+    --if there is a discrepancy then start over with the next bIdx
+    --print( "VALUES:", "A[[", aVal, aTest, "B[[", bVal, bTest )
+    if ( aVal < aTest ) then
+      if ( bVal < bTest ) then 
+        --print( "A LEFT B LEFT", aVal, aTest, bVal, bTest ) --debug
+        treeIdx = 2 * treeIdx
+      else
+        --print( "A LEFT B RIGHT", aVal, aTest, bVal, bTest ) --debug
+        newCandidate = true
+      end
+    elseif ( bVal > bTest ) then 
+      --print( "A RIGHT B RIGHT", aVal, aTest, bVal, bTest ) --debug
+      treeIdx = 2 * treeIdx + 1 
+    else 
+      --print( "A RIGHT B LEFT", aVal, aTest, bVal, bTest ) --debug
+      newCandidate = true
+    end
+  end
+
+  --Insert elements into tree.
+  aTree[treeIdx], bTree[treeIdx] = aVal, bVal
+
+  --Insert indices into sequence.
+  aSeq:insert( aIdx, bIdx )
+  bSeq:insert( bIdx, aIdx )
+
+end
+
+--In-order traversal of the tree. 
+--Returns points of twin-tree, sorted and interlaced, for graphing.
+-- { a_min, b_min, a_next, b_next, ...}
+do
+  local a, b, n, result
+  local function flat( idx )
+    if a[ 2 * idx ] then flat( 2 * idx ) end
+    result[ n ], result[ n + 1 ], n = a[idx], b[idx], n + 2
+    if a[ 2 * idx + 1 ] then flat( 2 * idx + 1 ) end
+  end
+
+  function tree.flat( aTree, bTree, t )
+    result = t or {}
+    a, b = aTree, bTree
+    n = 1
+    flat( 1 )
+    for i = n + 1, #result do result[i] = nil end
+    return result
+  end
+end
+
+function tree.build( n, av, bv )
+  local at, bt = {}, {}
+  local as, bs = sequence.new 'a', sequence.new 'b'
+  for i = 1, n do
+    insert( at, as, av, bt, bs, bv ) --Forward map.
+    insert( bt, bs, bv, at, as, av ) --Backward map.
+  end
+  return at, bt
+end
+
+do
+  local yield = assert( coroutine.yield )
+  function tree.buildIncremental( av, bv )
+    local at, bt = {}, {}
+    local as, bs = sequence.new 'a', sequence.new 'b'
+    local flattened = {}
+    repeat
+      insert( at, as, av, bt, bs, bv ) --Forward map.
+      insert( bt, bs, bv, at, as, av ) --Backward map.
+    until yield( tree.flat( at, bt, flattened ) ) and false
+  end
+end
+
+return tree