rapport/enumerations.lua

local enum = require 'enumeration'
local ffi = require 'ffi'

local function NadFFI( p )
  local rmt = { 
      __lt = function(a, b) return a.n / a.d < b.n / b.d end,
      __eq = function(a, b) return a.n == b.n and a.d == b.d end,
      __call = function(t) return t.n / t.d end,
      __tostring = function(t) return string.format( "<%d/%d>", t.n, t.d ) end
  }
  ffi.cdef "typedef struct { uint32_t n, d; } rational;"
  ffi.metatype( "rational", rmt )
  local list = ffi.new( "rational[?]", 123456 )
  local d, i = 1, 1
  repeat
    d = d * p
    for n = d - 1, 1, -p do
      list[i].n = n
      list[i].d = d
      i = i + 1
    end
  until i > 123456 / p
  return list
end

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

local function Nad( p )
  local n, d = -1, p
  return enum( function()
      n = n + p
      if n > d then n, d = 1, d * p end
      return n / d
    end )
end

return { 
  Dyad = Dyad,
  NadFFI = NadFFI,
  Nad = Nad,
  SternBrocot = SternBrocot,
  CalkinWilf = CalkinWilf,
  Random = Random,
  QuadIrr = false,
}