PostGIS Mandelbrot

postgis-mandelbrot.sql

with recursive
  dims (r1, r2, i1, i2, s, it, p) as (values (-2::float, 1::float, -1.5::float, 1.5::float, 0.01::float, 100, 256.0::float)),
  sprites (s) as (values (st_makepolygon(st_geomfromtext('linestring (0 0, 0 1, 1 1, 1 0, 0 0)')))),
  n1 (r, i) as (select r, i from dims, generate_series((r1 / s)::int, (r2 / s)::int) r, generate_series((i1 / s)::int, (i2 / s)::int) i),
  n2 (r, i) as (select r::float * s::float, i::float * s::float from dims, n1),
  l (cr, ci, zr, zi, g, it, p) as (
    select r, i, 0::float, 0::float, 0, it, p from n2, dims
    union all
    select cr, ci, zr*zr - zi*zi + cr, 2 * zr * zi + ci, g + 1, it, p from l where g < it and zr*zr + zi*zi < p
  ),
  x (cr, ci, zr, zi, g) as (
    select distinct on (cr, ci) cr, ci, zr, zi, g
    from l
    order by cr, ci, g desc
  )
select st_union(st_translate(sprites.s, round(cr / dims.s), round(ci / dims.s)))
from x, sprites, dims
where zr*zr + zi*zi > p