Benchmarks from benchmarks game

.gitignore

@@ -2,6 +2,9 @@ src/lua
 src/luac
 src/luaot
 
+experiments/*_fast.c
+experiments/*_fast.so
+
 *.a
 *.o
 *.so

experiments/binarytrees.lua

@@ -0,0 +1,79 @@
+-- Binary trees benchmark from benchmarks game
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/description/binarytrees.html
+--
+-- This is a small GC-focused benchmark, with not much room for clever tricks.
+-- The code used here is based on the one contributed by Mike Pall:
+--  * The leaf nodes are now store `false` instead of being empty tables.
+--    This is more consistent with what is done in other languages.
+--  * Use print instead of io.write
+--  * Use bitshift operator in Lua 5.4 and Pallene
+--
+-- Expected output (N = 21)
+--   stretch tree of depth 22	      check: 8388607
+--   2097152  trees of depth 4	      check: 65011712
+--   524288	  trees of depth 6	      check: 66584576
+--   131072	  trees of depth 8	      check: 66977792
+--   32768	  trees of depth 10	      check: 67076096
+--   8192	  trees of depth 12	      check: 67100672
+--   2048	  trees of depth 14	      check: 67106816
+--   512      trees of depth 16	      check: 67108352
+--   128      trees of depth 18	      check: 67108736
+--   32       trees of depth 20       check: 67108832
+--   long lived tree of depth 21      check: 4194303
+
+local function BottomUpTree(depth)
+    if depth > 0 then
+        depth = depth - 1
+        local left  = BottomUpTree(depth)
+        local right = BottomUpTree(depth)
+        return { left, right }
+    else
+        return { false, false }
+    end
+end
+
+local function ItemCheck(tree)
+    if tree[1] then
+        return 1 + ItemCheck(tree[1]) + ItemCheck(tree[2])
+    else
+        return 1
+    end
+end
+
+local function Stress(mindepth, maxdepth, depth)
+    local iterations = 1 << (maxdepth - depth + mindepth)
+    local check = 0
+    for _ = 1, iterations do
+        local t = BottomUpTree(depth)
+        check = check + ItemCheck(t)
+    end
+    return { iterations, check }
+end
+
+return function(N)
+    N = N or 10
+
+    local mindepth = 4
+    local maxdepth = math.max(mindepth + 2, N)
+
+    do
+        local stretchdepth = maxdepth + 1
+        local stretchtree = BottomUpTree(stretchdepth)
+        print(string.format("stretch tree of depth %d\t check: %d",
+            stretchdepth, ItemCheck(stretchtree)))
+    end
+
+    local longlivedtree = BottomUpTree(maxdepth)
+
+    for depth = mindepth, maxdepth, 2 do
+        local r = Stress(mindepth, maxdepth, depth)
+        local iterations = r[1]
+        local check      = r[2]
+        print(string.format("%d\t trees of depth %d\t check: %d",
+            iterations, depth, check))
+    end
+
+    print(string.format("long lived tree of depth %d\t check: %d",
+        maxdepth, ItemCheck(longlivedtree)))
+end
+

experiments/coro.lua

@@ -0,0 +1,13 @@
+local c = coroutine.wrap(function()
+    coroutine.yield(10)
+    coroutine.yield(20)
+    coroutine.yield(30)
+    coroutine.yield(40)
+end)
+
+return function()
+    print(c())
+    print(c())
+    print(c())
+    print(c())
+end

experiments/fannkuch.lua

@@ -0,0 +1,112 @@
+-- Fannkuch-redux benchmark from benchmarks game
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/description/fannkuchredux.html
+--
+-- Based on the C version found at:
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/fannkuchredux-gcc-1.html
+--
+-- Code by Hugo Gualandi, translated to Lua from the C version. The thing with
+-- the count vector is complicated, but we need to copy that from the original
+-- program to ensure that the permutation order is the right one.
+--
+-- * note: I made the function return a list of values because Pallene cannot
+--   return multiple values yet, or print the integer checksum.
+--
+-- Expected output (N = 7):
+--    228
+--    Pfannkuchen(7) = 16
+--
+-- Expected output (N = 12):
+--    3968050
+--    Pfannkuchen(12) = 65
+
+
+local function fannkuch(N)
+
+    local initial_perm = {}
+    for i = 1, N do
+        initial_perm[i] = i
+    end
+
+    local perm = {} -- Work copy, allocated once
+
+    local count = {}
+    count[1] = 0
+    local r = N
+
+    local perm_count = 0
+    local max_flips = 0
+    local checksum = 0
+
+    while true do
+
+        -- Flip the pancakes, working on a copy of the permutation
+
+        do
+            for i = 1, N do
+                perm[i] = initial_perm[i]
+            end
+
+            local flips_count = 0
+            local h = perm[1]
+            while h > 1 do
+                 local i = 1
+                 local j = h
+                 repeat
+                    local a = perm[i]
+                    local b = perm[j]
+                    perm[i] = b
+                    perm[j] = a
+                    i = i + 1
+                    j = j - 1
+                until i >= j
+
+                flips_count = flips_count + 1
+                h = perm[1]
+            end
+
+            if flips_count > max_flips then
+                max_flips = flips_count
+            end
+
+            if perm_count % 2 == 0 then
+                checksum = checksum + flips_count
+            else
+                checksum = checksum - flips_count
+            end
+        end
+
+        -- Go to next permutation
+
+        while r > 1 do
+            count[r] = r
+            r = r - 1
+        end
+
+        while true do
+            if r == N then
+                return { checksum, max_flips }
+            end
+
+            local tmp = initial_perm[1]
+            for i = 1, r do
+                initial_perm[i] = initial_perm[i+1]
+            end
+            initial_perm[r+1] = tmp
+
+            local r1 = r+1
+            count[r1] = count[r1] - 1
+            if count[r1] > 0 then break end
+            r = r1
+        end
+        perm_count = perm_count + 1
+    end
+end
+
+return function(N)
+    N = N or 7
+    local ret = fannkuch(N)
+    local checksum = ret[1]
+    local flips    = ret[2]
+    print(checksum)
+    print(string.format("Pfannkuchen(%d) = %d", N, flips))
+end

experiments/fasta.lua

@@ -0,0 +1,141 @@
+-- Fasta benchmark from benchmarks game
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/description/fasta.html
+--
+-- Translated from the Java and Python versions found at
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/fasta-java-2.html
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/fasta-python3-1.html
+--
+-- This differs from the Lua versions of the benchmark in some aspects
+--  * Don't use load/eval metaprogramming
+--  * Repeat fasta now works correctly when #alu < 60
+--  * Use linear search (actually faster than binary search in the tests)
+--  * Use // integer division
+
+local IM   = 139968
+local IA   = 3877
+local IC   = 29573
+
+local seed = 42
+local function random(max)
+    seed = (seed * IA + IC) % IM
+    return max * seed / IM;
+end
+
+local WIDTH = 60
+
+local function print_fasta_header(id, desc)
+    io.write(">" .. id .. " " .. desc .. "\n")
+end
+
+local function repeat_fasta(id, desc, alu, n)
+    print_fasta_header(id, desc)
+
+    local alusize = #alu
+
+    local aluwrap = alu .. alu
+    while #aluwrap < alusize + WIDTH do
+        aluwrap = aluwrap .. alu
+    end
+
+    local lines     = n // WIDTH
+    local last_line = n % WIDTH
+    local start = 0 -- (This index is 0-based bacause of the % operator)
+    for _ = 1, lines do
+        local stop = start + WIDTH
+        io.write(string.sub(aluwrap, start+1, stop))
+        io.write("\n")
+        start = stop % alusize
+    end
+    if last_line > 0 then
+        io.write(string.sub(aluwrap, start+1, start + last_line))
+        io.write("\n")
+    end
+end
+
+local function linear_search(ps, p)
+    for i = 1, #ps do
+        if ps[i]>= p then
+            return i
+        end
+    end
+    return 1
+end
+
+local function random_fasta(id, desc, frequencies, n)
+    print_fasta_header(id, desc)
+
+    -- Prepare the cummulative probability table
+    local nitems  = #frequencies
+    local letters = {}
+    local probs = {}
+    do
+        local total = 0.0
+        for i = 1, nitems do
+            local o = frequencies[i]
+            local c = o[1]
+            local p = o[2]
+            total = total + p
+            letters[i] = c
+            probs[i]   = total
+        end
+        probs[nitems] = 1.0
+    end
+
+    -- Generate the output
+    local col = 0
+    for _ = 1, n do
+        local ix = linear_search(probs, random(1.0))
+        local c = letters[ix]
+
+        io.write(c)
+        col = col + 1
+        if col >= WIDTH then
+            io.write("\n")
+            col = 0
+        end
+    end
+    if col > 0 then
+        io.write("\n")
+    end
+end
+
+local HUMAN_ALU =
+    "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG" ..
+    "GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA" ..
+    "CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT" ..
+    "ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA" ..
+    "GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG" ..
+    "AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC" ..
+    "AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA"
+
+local IUB = {
+    { 'a', 0.27 },
+    { 'c', 0.12 },
+    { 'g', 0.12 },
+    { 't', 0.27 },
+    { 'B', 0.02 },
+    { 'D', 0.02 },
+    { 'H', 0.02 },
+    { 'K', 0.02 },
+    { 'M', 0.02 },
+    { 'N', 0.02 },
+    { 'R', 0.02 },
+    { 'S', 0.02 },
+    { 'V', 0.02 },
+    { 'W', 0.02 },
+    { 'Y', 0.02 },
+}
+
+local HOMO_SAPIENS = {
+    { 'a', 0.3029549426680 },
+    { 'c', 0.1979883004921 },
+    { 'g', 0.1975473066391 },
+    { 't', 0.3015094502008 },
+}
+
+return function(N)
+    N = N or 100
+    repeat_fasta("ONE", "Homo sapiens alu", HUMAN_ALU, N*2)
+    random_fasta('TWO', 'IUB ambiguity codes', IUB, N*3)
+    random_fasta('THREE', 'Homo sapiens frequency', HOMO_SAPIENS, N*5)
+end

experiments/fib.lua

@@ -0,0 +1,9 @@
+local function fib(n)
+    if (n <= 1) then
+        return 1
+    else
+        return fib(n-1) + fib(n-2)
+    end
+end
+
+return fib

experiments/knucleotide.lua

@@ -0,0 +1,88 @@
+-- k-nucleotide norm benchmark from benchmarks game
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/description/knucleotide.html#knucleotide
+--
+-- Original Lua code by Mike Pall:
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/knucleotide-lua-2.html
+--
+-- Original Ruby code:
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/knucleotide-mri-2.html
+--
+-- Modifications by Hugo Gualandi:
+--  * reindented the code
+--  * rewrote in more idiomatic lua
+--
+-- Expected output (N = 5500):
+--    1.274224153
+
+
+--
+-- Reads the last sequence generated by the FASTA benchmark
+--
+local function read_input()
+    for line in io.lines() do
+        if string.sub(line, 1, 6) == ">THREE" then
+            break
+        end
+    end
+    local lines = {}
+    for line in io.lines() do
+        table.insert(lines, line)
+    end
+    return string.upper(table.concat(lines))
+end
+
+--
+-- Compute the frequency table for substrings of size `length`.
+-- The particular iteration order, using a moving reading frame, is mandatory.
+--
+local function compute_freqs(seq, length)
+    local freqs = setmetatable({}, { __index = function() return 0 end })
+    local n = #seq - length + 1
+    for f = 1, length do
+        for i = f, n, length do
+            local s = string.sub(seq, i, i+length-1)
+            freqs[s] = freqs[s] + 1
+        end
+    end
+    return n, freqs
+end
+
+--
+-- Print the most common subsequences, in decreasing order
+--
+local function sort_by_freq(seq, length)
+    local n, freq = compute_freqs(seq, length)
+    
+    local seqs = {}
+    for k, _ in pairs(freq) do
+        seqs[#seqs+1] = k
+    end
+    
+    table.sort(seqs, function(a, b)
+        local fa, fb = freq[a], freq[b]
+        return (fa == fb) and (a < b) or (fa > fb)
+    end)
+    
+    for _, c in ipairs(seqs) do
+        io.write(string.format("%s %0.3f\n", c, freq[c]*100/n))
+    end
+    io.write("\n")
+end
+
+local function single_freq(seq, s)
+    local n, freq = compute_freqs(seq, #s)
+    io.write(freq[s], "\t", s, "\n")
+end
+
+return function()
+    local seq = read_input()
+
+    sort_by_freq(seq, 1)
+    sort_by_freq(seq, 2)
+    
+    single_freq(seq, "GGT")
+    single_freq(seq, "GGTA")
+    single_freq(seq, "GGTATT")
+    single_freq(seq, "GGTATTTTAATT")
+    single_freq(seq, "GGTATTTTAATTTATAGT")
+end

experiments/knucleotidepall.lua

@@ -0,0 +1,69 @@
+-- The Computer Language Benchmarks Game
+-- https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
+-- contributed by Mike Pall
+
+local function kfrequency(seq, freq, k, frame)
+  local sub = string.sub
+  local k1 = k - 1
+  for i=frame,string.len(seq)-k1,k do
+    local c = sub(seq, i, i+k1)
+    freq[c] = freq[c] + 1
+  end
+end
+
+local function freqdefault()
+  return 0
+end
+
+local function count(seq, frag)
+  local k = string.len(frag)
+  local freq = setmetatable({}, { __index = freqdefault })
+  for frame=1,k do kfrequency(seq, freq, k, frame) end
+  io.write(freq[frag] or 0, "\t", frag, "\n")
+end
+
+local function frequency(seq, k)
+  local freq = setmetatable({}, { __index = freqdefault })
+  for frame=1,k do kfrequency(seq, freq, k, frame) end
+  local sfreq, sn = {}, 1
+  for c,v in pairs(freq) do sfreq[sn] = c; sn = sn + 1 end
+  table.sort(sfreq, function(a, b)
+    local fa, fb = freq[a], freq[b]
+    return fa == fb and a > b or fa > fb
+  end)
+  sum = string.len(seq)-k+1
+  for _,c in ipairs(sfreq) do
+    io.write(string.format("%s %0.3f\n", c, (freq[c]*100)/sum))
+  end
+  io.write("\n")
+end
+
+local function readseq()
+  local sub = string.sub
+  for line in io.lines() do
+    if sub(line, 1, 1) == ">" and sub(line, 2, 6) == "THREE" then break end
+  end
+  local lines, ln = {}, 0
+  for line in io.lines() do
+    local c = sub(line, 1, 1)
+    if c == ">" then
+      break
+    elseif c ~= ";" then
+      ln = ln + 1
+      lines[ln] = line
+    end
+  end
+  return string.upper(table.concat(lines, "", 1, ln))
+end
+
+
+return function()
+local seq = readseq()
+frequency(seq, 1)
+frequency(seq, 2)
+count(seq, "GGT")
+count(seq, "GGTA")
+count(seq, "GGTATT")
+count(seq, "GGTATTTTAATT")
+count(seq, "GGTATTTTAATTTATAGT")
+end    

experiments/mandelbrot.lua

@@ -0,0 +1,61 @@
+-- Mandelbrot benchmark from benchmarks game
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/description/mandelbrot.html
+--
+-- Translated from the Java version found at
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/mandelbrot-java-1.html
+--
+--  * I don't use any explicit buffering. In my experiments it speeds up the
+--    program by less than 10%, while considerably increasing the complexity.
+--  * The LuaJIT version needs to be separate, due to the lack of bitwise ops.
+--  * The output is in the Netpbm file format. Use an image viewer to view the picture.
+
+local function mandelbrot(N)
+    local bits  = 0
+    local nbits = 0
+
+    local delta = 2.0 / N
+    for y = 0, N-1 do
+        local Ci = y*delta - 1.0
+        for x = 0, N-1 do
+            local Cr = x*delta - 1.5
+
+            local bit = 0x1
+            local Zr  = 0.0
+            local Zi  = 0.0
+            local Zr2 = 0.0
+            local Zi2 = 0.0
+            for _ = 1, 50 do
+                Zi = 2.0 * Zr * Zi + Ci
+                Zr = Zr2 - Zi2 + Cr
+                Zi2 = Zi * Zi
+                Zr2 = Zr * Zr
+                if Zi2 + Zr2 > 4.0 then
+                    bit = 0x0
+                    break
+                end
+            end
+
+            bits = (bits << 1) | bit
+            nbits = nbits + 1
+
+            if nbits == 8 then
+                io.write(string.char(bits))
+                bits  = 0
+                nbits = 0
+            end
+        end
+
+        if nbits > 0 then
+            bits  = bits << (8 - nbits)
+            io.write(string.char(bits))
+            bits  = 0
+            nbits = 0
+        end
+    end
+end
+
+return function(N)
+    N = N or 100
+    io.write(string.format("P4\n%d %d\n", N, N))
+    mandelbrot(N)
+end

experiments/manorboy.lua

@@ -0,0 +1,24 @@
+-- The man or boy test was proposed by computer scientist Donald Knuth as a
+-- means of evaluating implementations of the ALGOL 60 programming language.
+-- The aim of the test was to distinguish compilers that correctly implemented
+-- "recursion and non-local references" from those that did not.
+-- https://rosettacode.org/wiki/Man_or_boy_test#Lua
+
+local function a(k,x1,x2,x3,x4,x5)
+  local function b()
+    k = k - 1
+    return a(k,b,x1,x2,x3,x4)
+  end
+   if k <= 0 then return x4() + x5() else return b() end
+end
+
+local function K(n)
+  return function()
+    return n
+  end
+end
+
+return function(N)
+    N = N or 10
+    print( a(N, K(1), K(-1), K(-1), K(1), K(0)) )
+end

experiments/nbody.lua

@@ -0,0 +1,168 @@
+-- N-body benchmark from benchmarks game
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/description/nbody.html
+--
+-- Original code by Mike Pall and Geoff Leyland, with modifications by Hugo
+-- Gualandi for correctness, clarity, and compatibility with Pallene.
+--   * Fix implementation of advance() to use two outer loops.
+--   * Use #xs instead of passing "n" as a parameter to all the functions
+--   * Use math.math.sqrt instead of caching math.sqrt in a local variable
+--   * Use 0.0 instead of 0 where appropriate.
+--   * Don't use multiple assignments.
+--
+-- Expected output (N = 1000):
+--   -0.169075164
+--   -0.169087605
+--
+-- Expected output (N = 50000000):
+--   -0.169075164
+--   -0.169059907
+
+local function new_body(x, y, z, vx, vy, vz, mass)
+    return {
+        x = x,
+        y = y,
+        z = z,
+        vx = vx,
+        vy = vy,
+        vz = vz,
+        mass = mass,
+    }
+end
+
+local function offset_momentum(bodies)
+    local n = #bodies
+    local px = 0.0
+    local py = 0.0
+    local pz = 0.0
+    for i= 1, n do
+      local bi = bodies[i]
+      local bim = bi.mass
+      px = px + (bi.vx * bim)
+      py = py + (bi.vy * bim)
+      pz = pz + (bi.vz * bim)
+    end
+
+    local sun = bodies[1]
+    local solar_mass = sun.mass
+    sun.vx = sun.vx - px / solar_mass
+    sun.vy = sun.vy - py / solar_mass
+    sun.vz = sun.vz - pz / solar_mass
+end
+
+local function advance(bodies, dt)
+    local n = #bodies
+    for i = 1, n do
+        local bi = bodies[i]
+        for j=i+1,n do
+          local bj = bodies[j]
+          local dx = bi.x - bj.x
+          local dy = bi.y - bj.y
+          local dz = bi.z - bj.z
+
+          local dist = math.sqrt(dx*dx + dy*dy + dz*dz)
+          local mag = dt / (dist * dist * dist)
+
+          local bjm = bj.mass * mag
+          bi.vx = bi.vx - (dx * bjm)
+          bi.vy = bi.vy - (dy * bjm)
+          bi.vz = bi.vz - (dz * bjm)
+
+          local bim = bi.mass * mag
+          bj.vx = bj.vx + (dx * bim)
+          bj.vy = bj.vy + (dy * bim)
+          bj.vz = bj.vz + (dz * bim)
+        end
+    end
+    for i = 1, n do
+        local bi = bodies[i]
+        bi.x = bi.x + dt * bi.vx
+        bi.y = bi.y + dt * bi.vy
+        bi.z = bi.z + dt * bi.vz
+    end
+end
+
+local function advance_multiple_steps(nsteps, bodies, dt)
+    for _ = 1, nsteps do
+        advance(bodies, dt)
+    end
+end
+
+local function energy(bodies)
+    local n = #bodies
+    local e = 0.0
+    for i = 1, n do
+        local bi = bodies[i]
+        local vx = bi.vx
+        local vy = bi.vy
+        local vz = bi.vz
+        e = e + (0.5 * bi.mass * (vx*vx + vy*vy + vz*vz))
+        for j = i+1, n do
+          local bj = bodies[j]
+          local dx = bi.x-bj.x
+          local dy = bi.y-bj.y
+          local dz = bi.z-bj.z
+          local distance = math.sqrt(dx*dx + dy*dy + dz*dz)
+          e = e - ((bi.mass * bj.mass) / distance)
+        end
+    end
+    return e
+end
+
+local PI = 3.141592653589793
+local SOLAR_MASS = 4 * PI * PI
+local DAYS_PER_YEAR = 365.24
+local bodies = {
+  -- Sun
+  new_body(
+    0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.0,
+    SOLAR_MASS),
+  -- Jupiter
+  new_body(
+     4.84143144246472090e+00,
+    -1.16032004402742839e+00,
+    -1.03622044471123109e-01,
+     1.66007664274403694e-03 * DAYS_PER_YEAR,
+     7.69901118419740425e-03 * DAYS_PER_YEAR,
+    -6.90460016972063023e-05 * DAYS_PER_YEAR,
+     9.54791938424326609e-04 * SOLAR_MASS ),
+  -- Saturn
+  new_body(
+     8.34336671824457987e+00,
+     4.12479856412430479e+00,
+    -4.03523417114321381e-01,
+    -2.76742510726862411e-03 * DAYS_PER_YEAR,
+     4.99852801234917238e-03 * DAYS_PER_YEAR,
+     2.30417297573763929e-05 * DAYS_PER_YEAR,
+     2.85885980666130812e-04 * SOLAR_MASS ),
+  -- Uranus
+  new_body(
+     1.28943695621391310e+01,
+    -1.51111514016986312e+01,
+    -2.23307578892655734e-01,
+     2.96460137564761618e-03 * DAYS_PER_YEAR,
+     2.37847173959480950e-03 * DAYS_PER_YEAR,
+    -2.96589568540237556e-05 * DAYS_PER_YEAR,
+     4.36624404335156298e-05 * SOLAR_MASS ),
+  -- Neptune
+  new_body(
+     1.53796971148509165e+01,
+    -2.59193146099879641e+01,
+     1.79258772950371181e-01,
+     2.68067772490389322e-03 * DAYS_PER_YEAR,
+     1.62824170038242295e-03 * DAYS_PER_YEAR,
+    -9.51592254519715870e-05 * DAYS_PER_YEAR,
+     5.15138902046611451e-05 * SOLAR_MASS ),
+}
+
+return function (N)
+    N = N or 50000000
+    offset_momentum(bodies)
+    print(string.format("%0.9f", energy(bodies)))
+    advance_multiple_steps(N, bodies, 0.01)
+    print(string.format("%0.9f", energy(bodies)))
+end

experiments/revcomp.lua

@@ -0,0 +1,77 @@
+-- Reverse Complement benchmark from benchmarks game
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/description/revcomp.html
+--
+-- Original LUa code:
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/revcomp-lua-4.html
+--
+-- Adapted by Hugo Gualandi
+--  * Removed the eval-based metaprogramming
+--  * Don't do complicated logic with buffers
+--  * Use string.gsub
+
+-- 
+-- DNA complement table
+--
+local R = {}
+do
+    local codes = "ACGTUMRWSYKVHDBN"
+    local cmpls = "TGCAAKYWSRMBDHVN"
+    for i = 1, #codes do
+        local s = string.sub(codes, i, i)
+        local t = string.sub(cmpls, i, i)
+        R[s] = t
+        R[string.lower(s)] = t
+    end
+end
+
+
+--
+-- Print a long DNA sequence in multiple lines,
+-- with 60 columns per line
+--
+local function wrap(str)
+    local n = #str
+    local r = n
+    local i = 1
+    while r >= 60 do
+        io.write(string.sub(str, i, i+59), "\n")
+        i = i + 60
+        r = r - 60
+    end
+    if r > 0 then
+        io.write(string.sub(str, i, n), "\n")
+    end
+end
+
+--
+-- Print the reverse-complement of the input sequences
+--
+local function revcomp(infile)
+    local n = 0
+    local names  = {}
+    local liness = {}
+    
+    local lines
+    for line in infile:lines() do
+        if string.sub(line, 1, 1) == ">" then
+            n = n + 1
+            lines = {}
+            names[n] = line
+            liness[n] = lines
+        else
+            lines[#lines+1] = line
+        end
+    end
+    
+    for i = 1, n do
+        local seq = table.concat(liness[i])
+        local rev = string.reverse(seq)
+        local cpl = string.gsub(rev, '.', R)
+        io.write(names[i], "\n")
+        wrap(cpl)
+    end
+end
+
+return function()
+    revcomp(io.stdin)
+end

experiments/spectralnorm.lua

@@ -0,0 +1,89 @@
+-- Spectral norm benchmark from benchmarks game
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/description/spectralnorm.html
+--
+-- Original C# code:
+-- https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/spectralnorm-csharpcore-1.html
+--
+-- Code by Hugo Gualandi. It is translated directly from the C# specification
+-- with the following main differences:
+--   - The A() function uses 1-based indexing instead of 0-based
+--   - The A() function multiplies by 0.5 instead of doing an integer division.
+--     According to my measurements, this is not a significant difference, and
+--     the multiplication by 0.5 has the advantage that it works on LuaJIT too.
+--   - Some of the "out" tables not initialized with zeroes.
+--
+-- Expected output (N = 5500):
+--    1.274224153
+
+
+-- Return A[i][j], for the infinite matrix A
+--
+--  A = 1/1  1/2  1/4 ...
+--      1/3  1/5  ... ...
+--      1/6  ...  ... ...
+--      ...  ...  ... ...
+local function A(i, j)
+    local ij = i + j
+    return 1.0 / ((ij-1) * (ij-2) * 0.5 + i)
+end
+
+-- Multiply vector v by matrix A
+local function MultiplyAv(N, v, out)
+    for i = 1, N do
+        local s = 0.0
+        for j = 1, N do
+            s = s + A(i,j) * v[j]
+        end
+        out[i] = s
+    end
+end
+
+-- Multiply vector v by matrix A transposed
+local function MultiplyAtv(N, v, out)
+    for i=1, N do
+        local s = 0.0
+        for j = 1, N do
+            s = s + A(j,i) * v[j]
+        end
+        out[i] = s
+    end
+end
+
+-- Multiply vector v by matrix A and then by matrix A transposed
+local function MultiplyAtAv(N, v, out)
+    local u = {}
+    MultiplyAv(N, v, u)
+    MultiplyAtv(N, u, out)
+end
+
+local function Approximate(N)
+    -- Create unit vector
+    local u = {}
+    for i = 1, N do
+        u[i] = 1.0
+    end
+
+    -- 20 steps of the power method
+    local v = {}
+    for _ = 1, 10 do
+        MultiplyAtAv(N, u, v)
+        MultiplyAtAv(N, v, u)
+    end
+
+    local vBv = 0.0
+    local vv  = 0.0
+    for i = 1, N do
+        local ui = u[i]
+        local vi = v[i]
+        vBv = vBv + ui*vi
+        vv  = vv  + vi*vi
+    end
+
+    return math.sqrt(vBv/vv)
+end
+
+return function(N)
+    N = N or 5500
+    local res = Approximate(N)
+    print(string.format("%0.9f", res))
+end