int width, height;
CLKernel mandelbrot;
CLKernel dmandelbrot;
+ CLKernel ddmandelbrot;
Calculate(int width, int height) throws CLException {
devs = new CLDevice[]{CLPlatform.getBestDevice(CL_DEVICE_TYPE_ALL)};
dmandelbrot = prog.createKernel("dmandelbrot");
dmandelbrot.setArg(0, image);
dmandelbrot.setArg(1, width, height);
+
+ ddmandelbrot = prog.createKernel("ddmandelbrot");
+ ddmandelbrot.setArg(0, image);
+ ddmandelbrot.setArg(1, width, height);
}
void release() {
- CLObject.release(dmandelbrot, mandelbrot, prog, image, q, cl);
+ CLObject.release(ddmandelbrot, dmandelbrot, mandelbrot, prog, image, q, cl);
CLObject.release(devs);
}
} catch (CLException ex) {
ex.printStackTrace();
}
+ }
+ public void rundd(ByteBuffer buffer, DD ox, DD oy, DD scale, double arg, int m) {
+ try {
+ ddmandelbrot.setArg(2, ox.x, ox.y, oy.x, oy.y);
+ ddmandelbrot.setArg(3, scale.x, scale.y);
+ ddmandelbrot.setArg(4, sin(arg), cos(arg));
+ ddmandelbrot.setArg(5, m);
+
+ q.enqueue2DKernel(ddmandelbrot, 0, 0, width, height, 8, 8, null, null);
+ q.enqueueReadBuffer(image, true, 0, buffer.limit(), buffer, null, null);
+ q.finish();
+ } catch (CLException ex) {
+ ex.printStackTrace();
+ }
}
}
}
+ enum Mode {
+ FLOAT,
+ DOUBLE,
+ DD,
+
+ }
+
// some from here, i made up the names http://paulbourke.net/fractals/mandelbrot/
static Target[] targets = {
new Target("jewel", -0.7746806106269039, -0.1374168856037867, 30),
};
Target target = targets[0];
- boolean isdouble = true;
+ Mode mode = Mode.DD;
double alpha = 0.005;
double beta = 0;
double gamma = 80;
getParameters().getUnnamed().forEach((c) -> {
switch (c) {
case "-f":
- isdouble = false;
+ mode = mode.FLOAT;
break;
case "-d":
- isdouble = true;
+ mode = mode.DOUBLE;
+ break;
+ case "-dd":
+ mode = mode.DD;
break;
case "--sync":
case "-s":
image = new WritableImage(width, height);
if (beta == 0.0)
- beta = isdouble ? target.beta : Math.min(12, target.beta);
+ beta = mode != Mode.FLOAT ? target.beta : Math.min(12, target.beta);
ImageView iv = new ImageView(image);
Text fps = new Text();
double scale = alpha / Math.exp(frac * beta);
int depth = (int)(-Math.log(scale) * gamma);
- recalc(target, isdouble, scale, arg, depth);
+ DD dscale = DD.div(new DD(alpha), DD.exp(new DD(frac * beta)));
+ int ddepth = (int)DD.neg(DD.mul(DD.log(dscale), gamma)).x;
+
+ System.out.printf("scale: %f dscale: %f\n", scale, dscale.x);
+ System.out.printf("depth: %d ddepth: %d\n", depth, ddepth);
+
+ if (mode == Mode.DD) {
+ recalc(target, mode, dscale, arg, ddepth);
+ } else {
+ recalc(target, mode, scale, arg, depth);
+ }
}
}
double arg = 0;
int lastiter;
- void recalc(Target target, boolean isdouble, double scale, double arg, int m) {
+ void recalc(Target target, Mode mode, double scale, double arg, int m) {
double x = target.cx - width * scale * 0.5;
double y = target.cy - height * scale * 0.5;
busy = true;
queue.submit(() -> {
- if (isdouble)
+ switch (mode) {
+ case DOUBLE:
+ case DD:
calc.rund(buffer, x, y, scale, scale, arg, m);
- else
+ break;
+ case FLOAT:
calc.run(buffer, x, y, scale, scale, arg, m);
+ }
Platform.runLater(() -> {
nframes++;
});
});
} else {
-
- if (isdouble)
+ switch (mode) {
+ case DOUBLE:
+ case DD:
calc.rund(buffer, x, y, scale, scale, arg, m);
- else
+ break;
+ case FLOAT:
calc.run(buffer, x, y, scale, scale, arg, m);
+ }
+ nframes++;
+ image.getPixelWriter()
+ .setPixels(0, 0, width, height, PixelFormat.getByteBgraInstance(), buffer, width * 4);
+ }
+ }
+
+ void recalc(Target target, Mode mode, DD scale, double arg, int m) {
+ //double x = target.cx - width * scale * 0.5;
+ //double y = target.cy - height * scale * 0.5;
+
+ DD x = DD.sub(new DD(target.cx), DD.mul_pwr2(DD.mul(scale, width), 0.5));
+ DD y = DD.sub(new DD(target.cy), DD.mul_pwr2(DD.mul(scale, height), 0.5));
+
+ lastiter = m;
+
+ if (!synchronous) {
+ if (busy)
+ return;
+ busy = true;
+
+ queue.submit(() -> {
+ calc.rundd(buffer, x, y, scale, arg, m);
+
+ Platform.runLater(() -> {
+ nframes++;
+ image.getPixelWriter()
+ .setPixels(0, 0, width, height, PixelFormat.getByteBgraInstance(), buffer, width * 4);
+ busy = false;
+ });
+ });
+ } else {
+ switch (mode) {
+ case DD:
+ calc.rundd(buffer, x, y, scale, arg, m);
+ break;
+ }
nframes++;
image.getPixelWriter()
.setPixels(0, 0, width, height, PixelFormat.getByteBgraInstance(), buffer, width * 4);
image[dx + dy * size.x] = convert_uchar4(colour * 255.0f);
}
}
+
+// double double?
+// https://github.com/lumianph/gpuprec/tree/master/gqd
+typedef double2 dd_real;
+
+#define _GQD_SPLITTER (134217729.0) // = 2^27 + 1
+#define _GQD_SPLIT_THRESH (6.69692879491417e+299) // = 2^996
+
+dd_real make_d(double x) {
+ return (dd_real)( x, 0.0 );
+}
+
+dd_real make_dd(double x, double y) {
+ return (dd_real)( x, 0.0 );
+}
+
+
+//computs fl( a + b ) and err( a + b ), assumes |a| > |b|
+
+double quick_two_sum(double a, double b, double *err) {
+
+ if (b == 0.0) {
+ *err = 0.0;
+ return (a + b);
+ }
+
+ double s = a + b;
+ *err = b - (s - a);
+
+ return s;
+}
+
+double two_sum(double a, double b, double *err) {
+
+ if ((a == 0.0) || (b == 0.0)) {
+ *err = 0.0;
+ return (a + b);
+ }
+
+ double s = a + b;
+ double bb = s - a;
+ *err = (a - (s - bb)) + (b - bb);
+
+ return s;
+}
+
+//computes fl( a - b ) and err( a - b ), assumes |a| >= |b|
+double quick_two_diff(double a, double b, double *err) {
+ if (a == b) {
+ *err = 0.0;
+ return 0.0;
+ }
+
+ double s;
+
+ /*
+ if(fabs((a-b)/a) < GPU_D_EPS) {
+ s = 0.0;
+ err = 0.0;
+ return s;
+ }
+ */
+
+ s = a - b;
+ *err = (a - s) - b;
+ return s;
+}
+
+//computes fl( a - b ) and err( a - b )
+double two_diff(double a, double b, double *err) {
+ if (a == b) {
+ *err = 0.0;
+ return 0.0;
+ }
+
+ double s = a - b;
+
+ /*
+ if(fabs((a-b)/a) < GPU_D_EPS) {
+ s = 0.0;
+ err = 0.0;
+ return s;
+ }
+ */
+
+ double bb = s - a;
+ *err = (a - (s - bb)) - (b + bb);
+ return s;
+}
+
+/* double-double = double + double */
+dd_real ddadd(double a, double b) {
+ double s, e;
+ s = two_sum(a, b, &e);
+ return make_dd(s, e);
+}
+
+dd_real negative(const dd_real a) {
+ return make_dd(-a.x, -a.y);
+}
+
+/* double-double + double */
+dd_real dd_add_d(const dd_real a, double b) {
+ double s1, s2;
+ s1 = two_sum(a.x, b, &s2);
+ s2 += a.y;
+ s1 = quick_two_sum(s1, s2, &s2);
+ return make_dd(s1, s2);
+}
+
+dd_real dd_add_dd(const dd_real a, dd_real b) {
+ double s, e;
+
+ s = two_sum(a.x, b.x, &e);
+ e += (a.y + b.y);
+ s = quick_two_sum(s, e, &e);
+ return make_dd(s, e);
+}
+/*********** Subtractions *********/
+
+dd_real dd_sub_dd(const dd_real a, const dd_real b) {
+
+ double s, e;
+ s = two_diff(a.x, b.x, &e);
+ //return make_dd(s, e);
+ e += a.y;
+ e -= b.y;
+ s = quick_two_sum(s, e, &e);
+ return make_dd(s, e);
+}
+
+
+// Computes high word and lo word of a
+void dd_split(double a, double *hi, double *lo) {
+ double temp;
+ if (a > _GQD_SPLIT_THRESH || a < -_GQD_SPLIT_THRESH) {
+ a *= 3.7252902984619140625e-09; // 2^-28
+ temp = _GQD_SPLITTER * a;
+ *hi = temp - (temp - a);
+ *lo = a - *hi;
+ *hi *= 268435456.0; // 2^28
+ *lo *= 268435456.0; // 2^28
+ } else {
+ temp = _GQD_SPLITTER * a;
+ *hi = temp - (temp - a);
+ *lo = a - *hi;
+ }
+}
+
+/* Computes fl(a*b) and err(a*b). */
+double two_prod(double a, double b, double *err) {
+ double a_hi, a_lo, b_hi, b_lo;
+ double p = a * b;
+
+ dd_split(a, &a_hi, &a_lo);
+ dd_split(b, &b_hi, &b_lo);
+
+ //err = (a_hi*b_hi) - p + (a_hi*b_lo) + (a_lo*b_hi) + (a_lo*b_lo);
+ *err = (a_hi * b_hi) - p + (a_hi * b_lo) + (a_lo * b_hi) + (a_lo * b_lo);
+
+ return p;
+}
+
+/* Computes fl(a*a) and err(a*a). Faster than the above method. */
+double dd_two_sqr(double a, double *err) {
+ double hi, lo;
+ double q = a * a;
+
+ dd_split(a, &hi, &lo);
+ *err = ((hi * hi - q) + 2.0 * hi * lo) + lo * lo;
+ return q;
+}
+
+/*********** Squaring **********/
+dd_real dd_sqr(const dd_real a) {
+ double p1, p2;
+ double s1, s2;
+ p1 = dd_two_sqr(a.x, &p2);
+ p2 += (2.0 * a.x * a.y);
+ p2 += (a.y * a.y);
+ s1 = quick_two_sum(p1, p2, &s2);
+ return make_dd(s1, s2);
+}
+
+dd_real d_sqr(double a) {
+ double p1, p2;
+ p1 = dd_two_sqr(a, &p2);
+ return make_dd(p1, p2);
+}
+
+dd_real ldexp_dd(const dd_real a, int exp) {
+ return make_dd(ldexp(a.x, exp), ldexp(a.y, exp));
+}
+
+/* double-double * double, where double is a power of 2. */
+dd_real dd_mul_pwr2(const dd_real a, double b) {
+ return make_dd(a.x * b, a.y * b);
+}
+
+/* double-double * double-double */
+dd_real dd_mul_dd(const dd_real a, const dd_real b) {
+ double p1, p2;
+
+ p1 = two_prod(a.x, b.x, &p2);
+ p2 += (a.x * b.y + a.y * b.x);
+
+ p1 = quick_two_sum(p1, p2, &p2);
+
+ return make_dd(p1, p2);
+}
+
+/* double-double * double */
+dd_real dd_mul_d(const dd_real a, double b) {
+ double p1, p2;
+
+ p1 = two_prod(a.x, b, &p2);
+ p2 += a.y * b;
+ p1 = quick_two_sum(p1, p2, &p2);
+ return make_dd(p1, p2);
+}
+
+/* double-double < double */
+int dd_lt_d(const dd_real a, double b) {
+ return (a.x < b || (a.x == b && a.y < 0.0));
+}
+
+kernel void
+__attribute__((reqd_work_group_size(LWS_X, LWS_Y, 1)))
+ddmandelbrot(global uchar4 *image, int2 size, double4 origin, double2 scale, double2 sincos, int m) {
+ int dx = get_global_id(0);
+ int dy = get_global_id(1);
+
+ if (dx < size.x && dy < size.y) {
+ double2 s = convert_double2(size);
+ double2 c;
+
+ c = (double2)(dx, dy) - s * 0.5;
+ c = (double2)(c.x * sincos.x + c.y * sincos.y,
+ c.x * sincos.y - c.y * sincos.x);
+ c += s * 0.5;
+
+ dd_real cx = make_d(c.x);
+ dd_real cy = make_d(c.y);
+
+ //c *= scale;
+ //c += origin;
+
+ cx = dd_mul_dd(cx, scale);
+ cy = dd_mul_dd(cy, scale);
+
+ cx = dd_add_dd(cx, origin.s01);
+ cy = dd_add_dd(cy, origin.s23);
+
+ dd_real zx = make_d(0.0);
+ dd_real zy = make_d(0.0);
+
+ int n = 0;
+
+ do {
+ dd_real tx = dd_add_dd(dd_sub_dd(dd_sqr(zx), dd_sqr(zy)), cx);
+ dd_real ty = dd_add_dd(dd_mul_pwr2(dd_mul_dd(zx, zy), 2), cy);
+
+ zx = tx;
+ zy = ty;
+
+ //z = (double2)(z.x * z.x - z.y * z.y, 2.0f * z.x * z.y) + c;
+ n++;
+ } while (dd_lt_d(dd_add_dd(dd_sqr(zx), dd_sqr(zy)), 4) && n < m);
+
+ // we use a simple cosine palette to determine color:
+ // http://iquilezles.org/www/articles/palettes/palettes.htm
+ float t = n * 10.0f / m;
+ float3 d = (float3)(0.5f, 0.5f, 0.5f);
+ float3 e = (float3)(0.5f, 0.5f, 0.5f);
+ float3 f = (float3)(1.0f, 1.0f, 1.0f);
+ float3 g = (float3)(0.00f, 0.33f, 0.67f);
+ float4 colour = (float4)( d + e * cos(6.28318f * (f * t + g)), 1.0f);
+
+ if (convert_int(n) == m)
+ colour = (float4)(0, 0, 0, 1);
+
+ // store the rendered mandelbrot set into a storage buffer:
+ image[dx + dy * size.x] = convert_uchar4(colour * 255.0f);
+ }
+}