diff options
| -rw-r--r-- | bignum.lua | 465 | ||||
| -rw-r--r-- | cmd/viewer.lua | 114 | ||||
| -rw-r--r-- | conf.lua | 5 | ||||
| -rw-r--r-- | main.lua | 10 | ||||
| -rw-r--r-- | mandelbrot.glsl | 67 | ||||
| -rw-r--r-- | mandelbrot.lua | 21 |
6 files changed, 682 insertions, 0 deletions
diff --git a/bignum.lua b/bignum.lua new file mode 100644 index 0000000..886f3f2 --- /dev/null +++ b/bignum.lua @@ -0,0 +1,465 @@ +local ffi = require "ffi" +local bit = require "bit" + +local M = {} + +ffi.cdef [[ + typedef struct bigint { + int8_t sign; + uint32_t size; + uint32_t values[?]; + } bigint; +]] + +local int_mt = {} + +local function normalize(int) + while int.size > 1 do + if int.values[int.size - 1] ~= 0 then + break + end + int.size = int.size - 1 + end + if int.size == 1 and int.values[0] == 0 then + int.sign = 0 + end + return int +end + +local function unsigned_add(a, b) + local size = math.max(a.size, b.size) + 1 + local result = ffi.new("bigint", size) + + local carry = 0ULL + for i = 0, size - 1 do + local sum = carry + if i < a.size then + sum = sum + a.values[i] + end + if i < b.size then + sum = sum + b.values[i] + end + result.values[i] = sum + carry = bit.rshift(sum, 32) + end + + result.size = size + result.sign = 1 + return result +end + +local function unsigned_sub(a, b) + local result = ffi.new("bigint", a.size) + local carry = 0LL + for i = 0, a.size - 1 do + local diff = carry + a.values[i] + if i < b.size then + diff = diff - b.values[i] + end + if diff < 0 then + diff = diff + 2^32 + carry = -1LL + else + carry = 0LL + end + result.values[i] = diff + end + + result.size = a.size + result.sign = 1 + return result +end +-- some sort of bug or something occurs with this! +require"jit".off(unsigned_sub) + +local function unsigned_cmp(a, b) + if a.size > b.size then + return 1 + elseif a.size < b.size then + return -1 + end + for i = a.size - 1, 0, -1 do + if a.values[i] > b.values[i] then + return 1 + elseif b.values[i] > a.values[i] then + return -1 + end + end + return 0 +end + +local function signed_cmp(a, b) + if a.sign > b.sign then + return 1 + elseif a.sign < b.sign then + return -1 + end + if a.sign > 0 then + return unsigned_cmp(a, b) + else + return -unsigned_cmp(a, b) + end +end + +local function lshift(int, n) + if n == 0 then return int end + if int == 0 then return int end + assert(n > 0) + + local size = int.size + n / 32 + 1 + local result = ffi.new("bigint", size) + result.size = size + + for i = 0, int.size - 1 do + local val = bit.lshift(ffi.cast("uint64_t", int.values[i]), n % 32) + local index = i + math.floor(n / 32) + result.values[index] = bit.bor(val, result.values[index]) + result.values[index + 1] = bit.rshift(val, 32) + end + result.sign = int.sign + + normalize(result) + return result +end + +function int_mt.__add(a, b) + a = M.is_bigint(a) and a or M.int(a) + b = M.is_bigint(b) and b or M.int(b) + + local result + if a.sign < 0 and b.sign < 0 then + result = unsigned_add(a, b) + result.sign = -1 + elseif a.sign < 0 or b.sign < 0 then + local neg = a.sign < 0 and a or b + local pos = a.sign < 0 and b or a + if unsigned_cmp(pos, neg) > 0 then + result = unsigned_sub(pos, neg) + else + result = unsigned_sub(neg, pos) + result.sign = -1 + end + else + result = unsigned_add(a, b) + end + normalize(result) + return result +end + +function int_mt.__sub(a, b) + return a + -b +end + +function int_mt.__mul(a, b) + a = M.is_bigint(a) and a or M.int(a) + b = M.is_bigint(b) and b or M.int(b) + + local result = ffi.new("bigint", a.size + b.size + 1) + result.size = a.size + b.size + + for i = 0, b.size - 1 do + local carry = 0ULL + for j = 0, a.size - 1 do + local val = + carry + result.values[i + j] + a.values[j] * b.values[i] + carry = bit.rshift(val, 32) + result.values[i + j] = val + end + result.values[i + a.size] = carry + end + result.sign = a.sign * b.sign + + normalize(result) + return result +end + +function int_mt.__unm(int) + local result = ffi.new("bigint", int.size) + result.size = int.size + result.sign = -int.sign + for i = 0, int.size - 1 do + result.values[i] = int.values[i] + end + return result +end + +function int_mt.__tostring(int) + local s = {"0x"} + if int.sign < 0 then table.insert(s, 1, "-") end + + for i = int.size - 1, 0, -1 do + if i < int.size - 1 then + table.insert(s, ("%08x"):format(int.values[i])) + else + table.insert(s, ("%x"):format(int.values[i])) + end + end + return table.concat(s) +end + +function int_mt.__eq(a, b) + a = M.is_bigint(a) and a or M.int(a) + b = M.is_bigint(b) and b or M.int(b) + + return signed_cmp(a, b) == 0 +end + +function int_mt.__lt(a, b) + a = M.is_bigint(a) and a or M.int(a) + b = M.is_bigint(b) and b or M.int(b) + + return signed_cmp(a, b) < 0 +end + +ffi.metatype("bigint", int_mt) + +function M.int(i) + if M.is_bigint(i) then + local int = ffi.new("bigint", i.size) + int.size = i.size + int.sign = i.sign + for j = 0, i.size - 1 do + int.values[j] = i.values[j] + end + return int + end + + i = math.floor(i) + local size = 1 + if i ~= 0 and math.abs(i) > 2^32 then + size = math.max(math.floor(math.log(math.abs(i), 2) / 32), 1) + end + local int = ffi.new("bigint", size) + + if i < 0 then + int.sign = -1 + i = -i + elseif i == 0 then + int.sign = 0 + int.size = 1 + return int + else + int.sign = 1 + end + + int.size = 0 + while i > 0 do + int.values[int.size] = i % (2^32) + i = math.floor(i / 2^32) + int.size = int.size + 1 + end + + return int +end + +function M.is_bigint(v) + return type(v) == "cdata" and ffi.typeof(v) == ffi.typeof"bigint" +end + +local float_mt = {} + +local function truncate(f, precision) + local precision = precision or f.precision + + if f.mantissa == 0 then f.exp = 0 return f end + if f.mantissa.size == 1 then return f end + if precision < 1 then return M.float(0) end + + local to = math.max(f.mantissa.size - precision, 0) + while f.mantissa.values[to] == 0 do + to = to + 1 + end + + if to ~= 0 then + local new = M.float() + new.exp = f.exp + to * 32 + new.precision = precision + + new.mantissa = ffi.new("bigint", f.mantissa.size - to) + new.mantissa.size = f.mantissa.size - to + new.mantissa.sign = f.mantissa.sign + for i = to, f.mantissa.size - 1 do + new.mantissa.values[i - to] = f.mantissa.values[i] + end + + return new + else + return f + end +end + +local function cmp_float(a, b) + if a.mantissa.sign > b.mantissa.sign then + return 1 + elseif a.mantissa.sign < b.mantissa.sign then + return -1 + end + + local a_mag = a.mantissa.size + math.ceil(a.exp / 32) + local b_mag = b.mantissa.size + math.ceil(b.exp / 32) + if a_mag > b_mag then + return 1 * a.mantissa.sign + elseif a_mag < b_mag then + return -1 * a.mantissa.sign + end + + return (a - b).mantissa.sign +end + +local recursion = true +function float_mt.__add(a, b) + a = M.is_bigfloat(a) and a or M.float(a) + b = M.is_bigfloat(b) and b or M.float(b) + + local f = M.float() + f.precision = math.min(a.precision, b.precision) + f.exp = math.min(a.exp, b.exp) + + local a = truncate(a, math.floor(f.precision - (a.exp - f.exp) / 32)) + local b = truncate(b, math.floor(f.precision - (b.exp - f.exp) / 32)) + + local a_shifted = lshift(a.mantissa, a.exp - f.exp) + local b_shifted = lshift(b.mantissa, b.exp - f.exp) + f.mantissa = a_shifted + b_shifted + + return truncate(f) +end + +function float_mt.__sub(a, b) + return a + -b +end + +function float_mt.__mul(a, b) + a = M.is_bigfloat(a) and a or M.float(a) + b = M.is_bigfloat(b) and b or M.float(b) + + local f = M.float() + f.exp = a.exp + b.exp + f.mantissa = a.mantissa * b.mantissa + f.precision = math.min(a.precision, b.precision) + return truncate(f) +end + +function float_mt.__unm(f) + local new = M.float() + new.exp = f.exp + new.mantissa = -f.mantissa + new.precision = f.precision + return new +end + +function float_mt.__eq(a, b) + a = M.is_bigfloat(a) and a or M.float(a) + b = M.is_bigfloat(b) and b or M.float(b) + + return cmp_float(a, b) == 0 +end + +function float_mt.__lt(a, b) + a = M.is_bigfloat(a) and a or M.float(a) + b = M.is_bigfloat(b) and b or M.float(b) + + return cmp_float(a, b) < 0 +end + +local function digits_in_range(f, min, max, anchor) + anchor = anchor or min + + local digits = {0} + for i = min, max do + if (i - anchor) % 4 == 0 and i ~= min then + table.insert(digits, 0) + end + + local b + if i >= 0 and i < f.mantissa.size * 32 then + b = bit.rshift(f.mantissa.values[i / 32], i % 32) or 0 + else + b = 0 + end + b = bit.band(b, 1) + b = bit.lshift(b, (i - anchor) % 4) + + local digit = digits[#digits] + digit = bit.bor(digit, b) + digits[#digits] = digit + end + return digits +end + +function float_mt.__tostring(f) + local before_point = digits_in_range(f, -f.exp, f.mantissa.size * 32 - 1) + while before_point[#before_point] == 0 and #before_point > 1 do + table.remove(before_point, #before_point) + end + + local after_point = digits_in_range(f, 0, -f.exp - 1, -f.exp) + while after_point[1] == 0 and #after_point > 1 do + table.remove(after_point, 1) + end + + local result = {"0x"} + if f.mantissa < 0 then + table.insert(result, 1, "-") + end + local hex = "0123456789abcdef" + for i = #before_point, 1, -1 do + local digit = before_point[i] + table.insert(result, hex:sub(digit + 1, digit + 1)) + end + table.insert(result, ".") + for i = #after_point, 1, -1 do + local digit = after_point[i] + table.insert(result, hex:sub(digit + 1, digit + 1)) + end + return table.concat(result) +end + +function M.float(f, precision) + local new = setmetatable({}, float_mt) + new.precision = precision or math.huge + + if f == 0 or f == nil then + new.mantissa = M.int(0) + new.exp = 0 + return new + elseif M.is_bigfloat(f) then + new.mantissa = f.mantissa + new.exp = f.exp + return new + elseif M.is_bigint(f) then + new.mantissa = f + new.exp = 0 + return new + end + + local m, e = math.frexp(f) + new.mantissa, new.exp = M.int(m * 2^53), e - 53 + return truncate(new) +end + +function M.is_bigfloat(v) + return getmetatable(v) == float_mt +end + +function M.tonumber(n) + local result + if type(n) == "number" then + result = n + elseif M.is_bigint(n) then + result = 0 + for i = 0, n.size - 1 do + result = result + n.values[i] * 2^(32 * i) + end + result = result * n.sign + elseif M.is_bigfloat(n) then + result = 0 + for i = 0, n.mantissa.size - 1 do + result = result + n.mantissa.values[i] * 2^(32 * i + n.exp) + end + result = result * n.mantissa.sign + end + return result +end + +return M diff --git a/cmd/viewer.lua b/cmd/viewer.lua new file mode 100644 index 0000000..05925a5 --- /dev/null +++ b/cmd/viewer.lua @@ -0,0 +1,114 @@ +local bignum = require "bignum" +local mandelbrot = require "mandelbrot" + + +local shader = love.graphics.newShader "mandelbrot.glsl" +local canvas = love.graphics.newCanvas(love.graphics.getDimensions()) +local redraw = true + +local iter = 256 +local view_pos = {bignum.float(0), bignum.float(0)} +local zoom = 300 +local antialias = 1 + +local function draw_brot(zoom, x, y) + love.graphics.setShader(shader) + shader:send("iterations", iter) + shader:send("zoom", zoom) + shader:send("antialias", antialias) + local orbit = mandelbrot.orbit(iter, x, y) + shader:send("orbit", unpack(orbit)) + + local w, h = love.graphics.getDimensions() + local ox, oy = -w / 2 / zoom, -h / 2 / zoom + shader:send("offset", {ox, oy}) + + love.graphics.rectangle("fill", 0, 0, love.graphics.getDimensions()) + love.graphics.setShader() + +-- love.graphics.push() +-- love.graphics.translate(w / 2, h / 2) +-- love.graphics.scale(zoom) +-- love.graphics.translate(-bignum.tonumber(x), -bignum.tonumber(y)) +-- love.graphics.setLineWidth(1 / zoom) +-- love.graphics.setColor(1, 0.1, 1) +-- for i = 2, #orbit do +-- local x1, y1 = unpack(orbit[i - 1]) +-- local x2, y2 = unpack(orbit[i]) +-- love.graphics.line(x1, y1, x2, y2) +-- end +-- love.graphics.setColor(1, 1, 1) +-- love.graphics.pop() +end + +function love.draw() + local w, h = love.graphics.getDimensions() + local cw, ch = canvas:getDimensions() + if cw < w or ch < h then + canvas:release() + canvas = love.graphics.newCanvas(w, h) + end + if cw ~= w or ch ~= h then + redraw = true + end + + if redraw then + love.graphics.setCanvas(canvas) + draw_brot(zoom, unpack(view_pos)) + redraw = false + end + + love.graphics.setShader() + love.graphics.setCanvas() + love.graphics.draw(canvas) + + love.graphics.setColor(0, 0, 0, 0.5) + love.graphics.rectangle("fill", 0, 0, w, 25) + love.graphics.setColor(1, 1, 1) + local info = ("p: (%s, %s) z: %s i: %s a: %dX"):format( + view_pos[1], view_pos[2], zoom, iter, antialias * antialias) + love.graphics.print(info, 5, 5) + love.graphics.setColor(1, 1, 1) + + local x, y = unpack(view_pos) +end + +function love.mousemoved(x, y, dx, dy) + if love.mouse.isDown(1) then + view_pos[1], view_pos[2] = + view_pos[1] - dx / zoom, view_pos[2] - dy / zoom + --print("p", unpack(view_pos)) + redraw = true + end +end + +function love.wheelmoved(_, sy) + local x, y = unpack(view_pos) + local w, h = love.graphics.getDimensions() + local mx, my = love.mouse.getPosition() + mx, my = mx - w / 2, my - h / 2 + + x, y = x + mx / zoom, y + my / zoom + zoom = zoom + sy * (zoom / 15) + x, y = x - mx / zoom, y - my / zoom + view_pos[1], view_pos[2] = x, y + + redraw = true +end + +function love.keypressed(key) + if key == "up" then + if iter < 4096 then + iter = iter * 2 + redraw = true + end + elseif key == "down" then + if iter > 1 then + iter = iter / 2 + redraw = true + end + elseif key == "a" then + antialias = 1 + (antialias % 3) + redraw = true + end +end diff --git a/conf.lua b/conf.lua new file mode 100644 index 0000000..5697142 --- /dev/null +++ b/conf.lua @@ -0,0 +1,5 @@ +function love.conf(t) + t.identity = "brot" + t.window.title = "brot" + t.window.resizable = true +end diff --git a/main.lua b/main.lua new file mode 100644 index 0000000..291f490 --- /dev/null +++ b/main.lua @@ -0,0 +1,10 @@ +function love.load(args) + local cmd = "viewer" + if not love.filesystem.isFused() and args[1] == '-r' then + cmd = args[2] + end + local path = "cmd/"..cmd..".lua" + local fn = assert( + loadstring(assert(love.filesystem.read(path)), "@"..path)) + return fn(unpack(args, 3)) +end diff --git a/mandelbrot.glsl b/mandelbrot.glsl new file mode 100644 index 0000000..42279bd --- /dev/null +++ b/mandelbrot.glsl @@ -0,0 +1,67 @@ +uniform vec2 orbit[4096]; +uniform vec2 offset = vec2(0, 0); +uniform float zoom = 300; +uniform int iterations = 256; +uniform int antialias = 0; + +vec2 mul(vec2 a, vec2 b) { + return vec2(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x); +} + +vec2 sq(vec2 a) { + return vec2(a.x*a.x - a.y*a.y, 2 * a.x * a.y); +} + +float abssq(vec2 a) { + return a.x*a.x + a.y*a.y; +} + +int brot(vec2 c) { + vec2 dz = vec2(0, 0); + int ri = 0; + for (int i = 0; i < iterations; i++) { + dz = mul(mul(vec2(2, 0), dz), orbit[ri]) + sq(dz) + c; + ri++; + vec2 z = orbit[ri] + dz; + if (abssq(z) > 4) return i; + if (abssq(z) < abssq(dz) || abssq(orbit[ri]) > 4) { + dz = z; + ri = 0; + } + } + return -1; +} + +vec3 hsv2rgb(vec3 c) +{ + vec4 K = vec4(1.0, 2.0 / 3.0, 1.0 / 3.0, 3.0); + vec3 p = abs(fract(c.xxx + K.xyz) * 6.0 - K.www); + return c.z * mix(K.xxx, clamp(p - K.xxx, 0.0, 1.0), c.y); +} + +vec3 color_at(vec2 screen_pos) { + vec2 pos = offset + screen_pos / zoom; + int i = brot(pos); + if (i >= 0) { + float it = i; + float ri = it / iterations; + float h = mod(0.4 + log(1 + it) / 12.0, 1.0); + float s = 0.6 + sin(it / 30) * 0.2; + float v = 0.1 + pow(ri, 1.01) * 2.0; + return hsv2rgb(vec3(h, s, v)); + } else { + return vec3(0, 0, 0); + } +} + +vec4 effect(vec4 color, Image t, vec2 img_pos, vec2 screen_pos) { + vec3 result = vec3(0, 0, 0); + for (float x = 0; x < antialias; x++) { + for (float y = 0; y < antialias; y++) { + vec2 pos = screen_pos - 0.5 + vec2(x, y) / antialias; + result += color_at(pos) / (antialias * antialias); + } + } + return vec4(result, 1); +} + diff --git a/mandelbrot.lua b/mandelbrot.lua new file mode 100644 index 0000000..a94f46f --- /dev/null +++ b/mandelbrot.lua @@ -0,0 +1,21 @@ +local bignum = require "bignum" + +local M = {} + +local function sq(x, y) + return x*x - y*y, 2 * x * y +end + +function M.orbit(iter, x, y) + local o = {} + local zx, zy = bignum.float(0, 32), bignum.float(0, 32) + for i = 1, iter do + local rx, ry = bignum.tonumber(zx), bignum.tonumber(zy) + table.insert(o, {rx, ry}) + zx, zy = sq(zx, zy) + zx, zy = zx + x, zy + y + end + return o +end + +return M |