gruff-0.1: compute.cc
#include <cmath>
#include <stdint.h>
#include <stdlib.h>
#include <qd/fpu.h>
#include <qd/dd_real.h>
#include <qd/dd_inline.h>
#include <qd/qd_real.h>
#include <qd/qd_inline.h>
#define likely(x) __builtin_expect((x),1)
#define unlikely(x) __builtin_expect((x),0)
static inline float to_float(float x) { return x; }
static inline float to_float(double x) { return x; }
static inline float to_float(dd_real x) { return to_double(x); }
static inline float to_float(qd_real x) { return to_double(x); }
static inline float sqr(float x) { return x * x; }
static inline double sqr(double x) { return x * x; }
template <typename T>
struct iter {
T zx, zy, dx, dy; int32_t n; int16_t i, j;
};
template <typename T>
static int compute(int *stop, float *out_n, float *out_d, float *out_a, T in_cx, T in_cy, int in_level, int in_iters) {
unsigned int fpu;
fpu_fix_start(&fpu);
struct iter<T> *iters[2];
iters[0] = (struct iter<T> *) malloc(sizeof(struct iter<T>) * 256 * 256);
iters[1] = (struct iter<T> *) malloc(sizeof(struct iter<T>) * 256 * 256);
int active[2];
active[0] = 0;
active[1] = 0;
int from = 0;
int to = 1;
int max_iters = -1;
{ /* initialize border */
struct iter<T> *it = iters[from];
int n = active[from];
{ /* corners */
it->zx = it->zy = it->dx = it->dy = it->n = 0; it->i = 0; it->j = 0; out_n[256 * ((int)it->j) + it->i] = -1; it++; n++;
it->zx = it->zy = it->dx = it->dy = it->n = 0; it->i = 255; it->j = 0; out_n[256 * ((int)it->j) + it->i] = -1; it++; n++;
it->zx = it->zy = it->dx = it->dy = it->n = 0; it->i = 0; it->j = 255; out_n[256 * ((int)it->j) + it->i] = -1; it++; n++;
it->zx = it->zy = it->dx = it->dy = it->n = 0; it->i = 255; it->j = 255; out_n[256 * ((int)it->j) + it->i] = -1; it++; n++;
}
for (int i = 1; i < 255; ++i) { /* edges */
it->zx = it->zy = it->dx = it->dy = it->n = 0; it->i = 0; it->j = i; out_n[256 * ((int)it->j) + it->i] = -1; it++; n++;
it->zx = it->zy = it->dx = it->dy = it->n = 0; it->i = 255; it->j = i; out_n[256 * ((int)it->j) + it->i] = -1; it++; n++;
it->zx = it->zy = it->dx = it->dy = it->n = 0; it->i = i; it->j = 0; out_n[256 * ((int)it->j) + it->i] = -1; it++; n++;
it->zx = it->zy = it->dx = it->dy = it->n = 0; it->i = i; it->j = 255; out_n[256 * ((int)it->j) + it->i] = -1; it++; n++;
}
active[from] = n;
}
const T scale = T(1.0) / (T(32.0) * pow(T(2.0), T(in_level)));
int progress = 1;
int progress2 = 1;
int progressed = 0;
int min_iters = 0;
int step_iters = 64;
int retval = 0;
while (active[from] && (progressed ? progress2 : 1) && step_iters < in_iters) {
progress2 = 0;
progress = 1;
while (progress) {
progress = 0;
int o = 0;
for (int i = 0; i < active[from]; ++i) {
if (*stop) { goto cleanup; }
T zx = iters[from][i].zx;
T zy = iters[from][i].zy;
T dx = iters[from][i].dx;
T dy = iters[from][i].dy;
T cx = in_cx + scale * iters[from][i].i;
T cy = in_cy + scale * iters[from][i].j;
int32_t n = iters[from][i].n;
{ /* iterate */
const T er2 = 65536;
T zx2 = sqr(zx);
T zy2 = sqr(zy);
T z2 = zx2 + zy2;
T z2xy = 2 * zx * zy;
while (likely(n < step_iters && z2 < er2)) {
T zdzx = zx * dx - zy * dy;
T zdzy = zx * dy + zy * dx;
dx = 2 * zdzx + 1;
dy = 2 * zdzy;
zx = zx2 - zy2 + cx;
zy = z2xy + cy;
zx2 = sqr(zx);
zy2 = sqr(zy);
z2xy = 2 * zx * zy;
z2 = zx2 + zy2;
++n;
}
if (z2 >= er2) {
int k = ((int) iters[from][i].j) * 256 + iters[from][i].i;
out_n[k] = to_float(1 + n - log(log(z2) / log(er2))/log(2.0));
out_d[k] = to_float((log(z2) * sqrt(z2 / (sqr(dx) + sqr(dy)))) / scale);
out_a[k] = to_float(atan2(zy, zx));
for (int x = iters[from][i].i - 1; x <= iters[from][i].i + 1; ++x) {
if (x < 0 || 255 < x) continue;
for (int y = iters[from][i].j - 1; y <= iters[from][i].j + 1; ++y) {
if (y < 0 || 255 < y) continue;
k = y * 256 + x;
if (out_n[k] == 0) {
iters[to][o].zx = T(0);
iters[to][o].zy = T(0);
iters[to][o].dx = T(0);
iters[to][o].dy = T(0);
iters[to][o].n = 0;
iters[to][o].i = x;
iters[to][o].j = y;
out_n[k] = -1;
++o;
}
}
}
if (min_iters < n) min_iters = n;
if (max_iters < n) max_iters = n;
++progress;
} else {
iters[to][o].zx = zx;
iters[to][o].zy = zy;
iters[to][o].dx = dx;
iters[to][o].dy = dy;
iters[to][o].n = n;
iters[to][o].i = iters[from][i].i;
iters[to][o].j = iters[from][i].j;
++o;
}
}
}
active[to] = o;
int temp = from; from = to; to = temp;
progress2 = progress2 || progress;
}
step_iters *= 2;
progressed = progressed || progress2;
}
{ /* deinitialize border */
int k;
for (int i = 0; i < 256; ++i) { /* edges */
k = 256 * i + 0; if (out_n[k] < 0) out_n[k] = 0;
k = 256 * i + 255; if (out_n[k] < 0) out_n[k] = 0;
k = 256 * 0 + i; if (out_n[k] < 0) out_n[k] = 0;
k = 256 * 255 + i; if (out_n[k] < 0) out_n[k] = 0;
}
}
retval = max_iters;
cleanup:
free(iters[0]);
free(iters[1]);
fpu_fix_end(&fpu);
return retval;
}
extern "C" {
int compute_f32(int *stop, float *out_n, float *out_d, float *out_a, float in_cx, float in_cy, int in_level, int in_iters) {
return compute(stop, out_n, out_d, out_a, in_cx, in_cy, in_level, in_iters);
}
int compute_f64(int *stop, float *out_n, float *out_d, float *out_a, double in_cx, double in_cy, int in_level, int in_iters) {
return compute(stop, out_n, out_d, out_a, in_cx, in_cy, in_level, in_iters);
}
int compute_f128(int *stop, float *out_n, float *out_d, float *out_a, double *in_cx, double *in_cy, int in_level, int in_iters) {
return compute(stop, out_n, out_d, out_a, dd_real(in_cx), dd_real(in_cy), in_level, in_iters);
}
int compute_f256(int *stop, float *out_n, float *out_d, float *out_a, double *in_cx, double *in_cy, int in_level, int in_iters) {
return compute(stop, out_n, out_d, out_a, qd_real(in_cx), qd_real(in_cy), in_level, in_iters);
}
}