--[[
 Correctness and benchmark tests of various hamming weight implementations.
 See also http://lua-users.org/wiki/HammingWeight .
 David Manura, 2012-03.
--]]

-- utility section --

-- https://gist.github.com/2064991
local function memoize(f)
  local mt = {}
  local t = setmetatable({}, mt)
  function mt:__index(k)
    local v = f(k); t[k] = v
    return v
  end
  return t
end

-- https://gist.github.com/1414923
local function requireany(...)
  local errs = {}
  for i = 1, select('#', ...) do local name = select(i, ...)
    if type(name) ~= 'string' then return name, nil end
    local ok, mod = pcall(require, name)
    if ok then return mod, name end
    errs[#errs+1] = mod
  end
  error(table.concat(errs, '\n'), 2)
end

-- testing utility
local OP = {}
OP['=='] = function(a, b) return a == b end
local function check(compare, a, b, name)
  local comparef = OP[compare]
  if not comparef(a, b) then
    error(name..': '..compare..tostring(a)..' '..tostring(b))
  end
end

-- benchmarking utility.
-- Note: uses CPU time (clock function).
-- `bench` repeatedly calls `f` for at least `min_seconds` seconds (defaults to 2).
--   `f_count` is number of operations performed in `f` (defaults to 1).
--   Returns number of seconds per operation.
local clock = os.clock
local function bench(min_seconds, f, f_count)
  min_seconds = min_seconds or 2
  f_count = f_count or 1
  local t1 = clock()
  local ntimes = 1
  local count = 0
  while 1 do
    for i=1,ntimes do
      f()
    end
    count = count + ntimes
    local t2 = clock()
    if t2 - t1 >= min_seconds then
      local period = (t2 - t1) / (count * f_count)
      return period
    end
    ntimes = ntimes * 2  -- exponential increase
  end
end
local function dump_bench(name, seconds)
  local ts = {}
  for i=1,#seconds do
    ts[i] = ('%0.1E'):format(seconds[i])
  end
  print(('%10s: %s s'):format(name, table.concat(ts, ' ')))
end


-- main code section --

local bit = requireany('bit', 'bit32')
local rshift = bit.rshift
local band = bit.band
local extract8 = bit.extract or function(n, field, _) -- https://github.com/davidm/lua-bit-numberlua
  return band(rshift(n, field), 0xFF)
end

-- Hamming weight of 32-bit integer.
-- Simple implementation.
local function hw_simple(x)
  local sum = 0
  while x ~= 0 do
    sum = sum + band(x, 1)
    x = rshift(x, 1)
  end
  return sum
end

-- Hamming weight of 32-bit integer.
-- Implementation based on 32-bit version of popcount_2
-- http://en.wikipedia.org/wiki/Hamming_weight .
-- The modulo product (x * h01) in popcount_4 could overflow floating point.
local MASK1 = 0x55555555  -- repeating 01 pattern
local MASK2 = 0x33333333  -- repeating 0011 pattern
local MASK4 = 0x0F0F0F0F  -- repeating 00001111 pattern
local function hw_pc2(x)
  x = x - band(rshift(x, 1), MASK1)              -- 2 bit sums
  x = band(x, MASK2) + band(rshift(x, 2), MASK2) -- 4 bit sums
  x = band(x + rshift(x, 4), MASK4)              -- 8 bit sums
  return band(x + rshift(x, 8) + rshift(x, 16) + rshift(x, 24), 0xFF) -- sum of bytes
end

-- Hamming weight of 32-bit integer.
-- Implementation based on 8-bit table lookups (memoized).
local HW = memoize(hw_simple)
local function hw_t4(x)
  local n0 = extract8(x, 0, 8)
  local n1 = extract8(x, 8, 8)
  local n2 = extract8(x, 16, 8)
  local n3 = extract8(x, 24, 8)
  return HW[n0] + HW[n1] + HW[n2] + HW[n3]
end

-- Hamming weight of 32-bit integer.
-- Implementation based on
-- Peter Wegner 1960 - http://dx.doi.org/10.1145/367236.367286
local function hw_wegner(x)
  local sum = 0
  while x ~= 0 do
    x = band(x, x-1)
    sum = sum + 1
  end
  return sum
end

local function test_hw_value(f, x)
  check('==', hw_simple(x), f(x), x)
end

local function test_hw(f)
  for x=0,1E+3 do
    test_hw_value(f, x)
    test_hw_value(f, 0xFFFFFFFF - x)
  end
end

-- note: test both high and low numbers to avoid bias in Wegner version.
local function make_test_func(f)
  local M = 0xFFFFFFFF
  local count = 0; for _=0,1E+3,7 do count = count + 10 end
  return function()
    local result = 0 -- guard against unused computations optimizing away.
    for i=0,1E+3,7 do
      result = result + f(i) + f(i) + f(i) + f(i) + f(i)
      local i = M-i
      result = result + f(i) + f(i) + f(i) + f(i) + f(i)
    end
    return result
  end, count
end

local results = {}
local function bench_hw(f, name)
  results[name] = results[name] or {}
  table.insert(results[name], bench(2.0, make_test_func(f)))
  dump_bench(name, results[name])
end

-- test correctness
test_hw(hw_pc2)
test_hw(hw_t4)
test_hw(hw_wegner)


-- test performance
for i=1,3 do
  bench_hw(hw_simple, 'hw_simple')
  bench_hw(hw_pc2,    'hw_pc2')
  bench_hw(hw_t4,     'hw_t4')
  bench_hw(hw_wegner, 'hw_wegner')
  print '---'
end

print 'DONE'


--[[
(c) 2012 David Manura.

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
--]]