diff --git a/meson.build b/meson.build
index 73ada15..f9d7997 100644
--- a/meson.build
+++ b/meson.build
@@ -17,8 +17,14 @@ sources = [
   'src/ppm.c',
 ]
 
+cc = meson.get_compiler('c')
+deps = [
+  cc.find_library('m', required : false),
+]
+
 executable(
   'buddhabrot',
   sources : sources,
-  c_args : '-D_XOPEN_SOURCE',
+  c_args : '-D_XOPEN_SOURCE=500',
+  dependencies : deps,
 )
diff --git a/src/buddhabrot.c b/src/buddhabrot.c
index a20cc5f..51ac2aa 100644
--- a/src/buddhabrot.c
+++ b/src/buddhabrot.c
@@ -12,12 +12,15 @@
 #include "image.h"
 
 struct data {
-    float min_x;
-    float max_x;
-    float min_y;
-    float max_y;
+    struct window {
+        float min_x;
+        float max_x;
+        float min_y;
+        float max_y;
+    } window;
 
     size_t max_iter;
+    size_t repeats;
 
     size_t h;
     size_t w;
@@ -37,32 +40,104 @@ static size_t data_index(size_t h, size_t w, const struct data *d) {
     return h * d->w + w;
 }
 
-static float scale(float min, float max, size_t x, size_t width) {
-    return min + x * (max - min) / (width - 1);
-}
-
-static size_t inv_scale(float min, float max, float x, size_t width) {
-    return (x - min) / (max - min) * (width - 1);
-}
-
 static bool compute_orbit(struct my_complex z,
         struct orbit *o, size_t max_iter) {
     float r = 0.0;
     float i = 0.0;
 
     for (o->len = 0; o->len < max_iter; ++o->len) {
+        o->values[o->len].r = r;
+        o->values[o->len].i = i;
+
         float xtemp = r * r - i * i + z.r;
         i = 2 * r * i + z.i;
         r = xtemp;
 
         if (r * r + i * i >= 4)
-            return true; // Diverges
-
-        o->values[o->len].r = r;
-        o->values[o->len].i = i;
+            return false; // Diverges
     }
 
-    return false;
+    return true;
+}
+
+static float orbit_probability(const struct window *w, const struct orbit *o) {
+    for (size_t i = 0; i < o->len; ++i) {
+        if (w->min_x <= o->values[i].r && o->values[i].r <= w->max_x
+            && w->min_y <= o->values[i].i && o->values[i].i <= w->max_y)
+            return 1;
+    }
+
+    return FLT_MIN; // Cannot be 0 because of sampler
+}
+
+static float transition_probability(const struct orbit *src,
+        const struct orbit *dst, size_t max_iter) {
+    return (1 - (max_iter - src->len) / max_iter)
+        / (1 - (max_iter - dst->len) / max_iter);
+}
+
+static bool mutate_near(const struct orbit *current,
+        struct orbit *new, size_t max_iter) {
+    struct my_complex next = current->values[1]; // Always the first one after 0
+
+    const double magnifictation = 1; // FIXME: proportional to magnification
+
+    const double r1 = (1.f / magnifictation) * 0.0001;
+    const double r2 = (1.f / magnifictation) * 0.1;
+
+    const double phi =  M_PI * 2 * random() / LONG_MAX;
+    const double r = r2 * exp(-log(r2 / r1) * random() / LONG_MAX);
+
+    next.r += r * cos(phi);
+    next.i += r * sin(phi);
+
+    return compute_orbit(next, new, max_iter);
+}
+
+static bool mutate_any(struct orbit *new, size_t max_iter) {
+    struct my_complex next = {
+        2 - 4. * random() / LONG_MAX,
+        2 - 4. * random() / LONG_MAX,
+    };
+
+    while (next.r * next.r + next.i * next.i >= 4) {
+        next.r = 2 - 4. * random() / LONG_MAX;
+        next.i = 2 - 4. * random() / LONG_MAX;
+    }
+
+    return compute_orbit(next, new, max_iter);
+}
+
+static void mutate(const struct data *d, struct orbit *current) {
+    struct orbit *new =
+        calloc(1, d->max_iter * sizeof(*new->values) + sizeof(*new));
+    if (!new)
+        err(EXIT_FAILURE, "could not allocate orbit data");
+
+    // Only consider values that diverge
+    if (random() < LONG_MAX / 5)
+        while (!mutate_any(new, d->max_iter))
+            continue;
+    else
+        while (!mutate_near(current, new, d->max_iter))
+            continue;
+
+    const float prob_cur = orbit_probability(&d->window, current);
+    const float prob_new = orbit_probability(&d->window, new);
+
+    const float t_new_cur = transition_probability(new, current, d->max_iter);
+    const float t_cur_new = transition_probability(current, new, d->max_iter);
+
+    const float alpha = (prob_new * t_cur_new) / (prob_cur * t_new_cur);
+
+    if (alpha > (float)random() / LONG_MAX)
+        memcpy(current, new, d->max_iter * sizeof(*new->values) + sizeof(*new));
+
+    free(new);
+}
+
+static size_t inv_scale(float min, float max, float x, size_t width) {
+    return (x - min) / (max - min) * width;
 }
 
 static void compute_buddahbrot(struct data *d) {
@@ -71,25 +146,19 @@ static void compute_buddahbrot(struct data *d) {
     if (!o)
         err(EXIT_FAILURE, "could not allocate orbit data");
 
-    for (size_t i = 0; i < d->h; ++i) {
-        for (size_t j = 0; j < d->w; ++j) {
-            struct my_complex z = {
-                scale(d->min_x, d->max_x, j, d->w),
-                scale(d->min_y, d->max_y, i, d->h),
-            };
-            // Only keep if it diverges
-            if (!compute_orbit(z, o, d->max_iter))
-                continue;
+    while (true)
+        if (mutate_any(o, d->max_iter))
+            break; // Initialize with correct value
 
-            for (size_t i = 0; i < o->len; ++i) {
-                size_t w = inv_scale(d->min_x, d->max_x,
-                        o->values[i].r, d->w);
-                size_t h = inv_scale(d->min_y, d->max_y,
-                        o->values[i].i, d->h);
-                if (w >= d->w || h >= d->h)
-                    continue; // We've diverged, keep going just in case
-                d->buf[data_index(h, w, d)] += 1;
-            }
+    for (size_t i = 0; i < d->repeats; ++i) {
+        mutate(d, o); // Computes its orbit
+
+        for (size_t i = 0; i < o->len; ++i) {
+            size_t w = inv_scale(d->window.min_x, d->window.max_x,
+                    o->values[i].r, d->w);
+            size_t h = inv_scale(d->window.min_y, d->window.max_y,
+                    o->values[i].i, d->h);
+            d->buf[data_index(h, w, d)] += 1;
         }
     }
 
@@ -108,7 +177,7 @@ static void fillout_image(const struct data *d, struct image *i) {
     }
 }
 
-void buddhabrot(struct image *i, size_t max_iter) {
+void buddhabrot(struct image *i, size_t max_iter, size_t repeats) {
     struct data *d = calloc(1, sizeof(*d->buf) * i->w * i->h + sizeof(*d));
     if (!d)
         err(EXIT_FAILURE, "could not allocate buddhabrot data");
@@ -116,12 +185,13 @@ void buddhabrot(struct image *i, size_t max_iter) {
     d->h = i->h;
     d->w = i->w;
     d->max_iter = max_iter;
+    d->repeats = repeats;
 
     // FIXME" make it user-definable
-    d->min_x = -2.5;
-    d->max_x = 1.0;
-    d->min_y = -1.0;
-    d->max_y = 1.0;
+    d->window.min_x = -2.5;
+    d->window.max_x = 1.0;
+    d->window.min_y = -1.0;
+    d->window.max_y = 1.0;
 
     compute_buddahbrot(d);
     fillout_image(d, i);
diff --git a/src/buddhabrot.h b/src/buddhabrot.h
index 73f4772..0d50af9 100644
--- a/src/buddhabrot.h
+++ b/src/buddhabrot.h
@@ -5,6 +5,6 @@
 
 struct image; // Forward declaration
 
-void buddhabrot(struct image *i, size_t max_iter);
+void buddhabrot(struct image *i, size_t max_iter, size_t repeats);
 
 #endif /* !BUDDHABROT_H */
diff --git a/src/main.c b/src/main.c
index 251ded6..0d28d12 100644
--- a/src/main.c
+++ b/src/main.c
@@ -4,6 +4,7 @@
 
 #include "buddhabrot.h"
 #include "image.h"
+#include "buddhabrot.h"
 #include "mandelbrot.h"
 #include "options.h"
 #include "ppm.h"
@@ -17,7 +18,7 @@ int main(int argc, char *argv[]) {
 
     switch (opt.render) {
     case BUDDHABROT:
-        buddhabrot(image, opt.max_iter);
+        buddhabrot(image, opt.max_iter, 10000); // FIXME: allow user option
         break;
     case MANDELBROT:
         mandelbrot(image, opt.max_iter);