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, }