add doubledouble version

Makefile

@@ -1,4 +1,5 @@
 all:
-	gcc -Ofast -Wall -lquadmath mandel_quad.c -o mandel_quad
+	gcc -O3 -Wall -lquadmath mandel_quad.c -o mandel_quad
-	gcc -Ofast -Wall -lmpfr mandel_mpfr.c -o mandel_mpfr
+	gcc -O3 -Wall -lmpfr -lm mandel_mpfr.c -o mandel_mpfr
-	gcc -Ofast -Wall -lm mandel.c -o mandel
+	gcc -O3 -Wall -lm mandel.c -o mandel
+	gcc -O3 -Wall -lm mandel_dd.c -o mandel_dd

doubledouble.h

@@ -0,0 +1,143 @@
+#include <math.h>
+
+typedef struct {
+    double hi;
+    double lo;
+} DoubleDouble;
+
+inline DoubleDouble dd_new(double hi, double lo) {
+    DoubleDouble num;
+    num.hi = hi;
+    num.lo = lo;
+    return num;
+}
+
+unsigned int dd_get_ui(DoubleDouble num) {
+    return num.hi;
+}
+
+inline DoubleDouble dd_sqrt(DoubleDouble num) {
+    double a, b, c, d, e;
+    d = 1 / sqrt(num.hi);
+    e = num.hi * d;
+    a = 0x08000001 * e;
+    a += e - a;
+    b = e - a;
+    c = e * e;
+    b = ((a * a - c) + a * b * 2) + b * b;
+    a = num.hi - c;
+    c = num.hi - a;
+    c = (a + ((((num.hi - (c + a)) + (c - c)) + num.lo) - b)) * d * 0.5;
+    a = e + c;
+    b = e - a;
+    return dd_new(a, (e - (b + a)) + (b + c));
+}
+
+inline DoubleDouble dd_div(DoubleDouble num1, DoubleDouble num2) {
+    double a, b, c, d, e, f, g;
+    f = num1.hi / num2.hi;
+    a = 0x08000001 * num2.hi;
+    a += num2.hi - a;
+    b = num2.hi - a;
+    c = 0x08000001 * f;
+    c += f - c;
+    d = f - c;
+    e = num2.hi * f;
+    c = (((a * c - e) + (a * d + b * c)) + b * d) + num2.lo * f;
+    b = num1.lo - c;
+    d = num1.lo - b;
+    a = num1.hi - e;
+    e = (num1.hi - ((num1.hi - a) + a)) + b;
+    g = a + e;
+    e += (a - g) + ((num1.lo - (d + b)) + (d - c));
+    a = g + e;
+    b = a / num2.hi;
+    f += (e + (g - a)) / num2.hi;
+    a = f + b;
+    return dd_new(a, b + (f - a));
+}
+
+inline DoubleDouble dd_ui_div(unsigned int num1, DoubleDouble num2) {
+    return dd_div(dd_new(num1, 0), num2);
+}
+
+inline DoubleDouble dd_div_ui(DoubleDouble num1, unsigned int num2) {
+    return dd_div(num1, dd_new(num2, 0));
+}
+
+inline DoubleDouble dd_si_div(signed int num1, DoubleDouble num2) {
+    return dd_div(dd_new(num1, 0), num2);
+}
+
+inline DoubleDouble dd_div_si(DoubleDouble num1, signed int num2) {
+    return dd_div(num1, dd_new(num2, 0));
+}
+
+inline DoubleDouble dd_add(DoubleDouble num1, DoubleDouble num2) {
+    double a, b, c, d, e, f;
+    e = num1.hi + num2.hi;
+    d = num1.hi - e;
+    a = num1.lo + num2.lo;
+    f = num1.lo - a;
+    d = ((num1.hi - (d + e)) + (d + num2.hi)) + a;
+    b = e + d;
+    c = ((num1.lo - (f + a)) + (f + num2.lo)) + (d + (e - b));
+    a = b + c;
+    return dd_new(a, c + (b - a));
+}
+
+inline DoubleDouble dd_mul(DoubleDouble num1, DoubleDouble num2) {
+    double a, b, c, d, e;
+    a = 0x08000001 * num1.hi;
+    a += num1.hi - a;
+    b = num1.hi - a;
+    c = 0x08000001 * num2.hi;
+    c += num2.hi - c;
+    d = num2.hi - c;
+    e = num1.hi * num2.hi;
+    c = (((a * c - e) + (a * d + b * c)) + b * d) + (num1.lo * num2.hi + num1.hi * num2.lo);
+    a = e + c;
+    return dd_new(a, c + (e - a));
+}
+
+inline DoubleDouble dd_mul2(DoubleDouble num1, DoubleDouble num2) {
+    double a, b, c, d, e;
+    a = 0x08000001 * num1.hi;
+    a += num1.hi - a;
+    b = num1.hi - a;
+    c = 0x08000001 * num2.hi;
+    c += num2.hi - c;
+    d = num2.hi - c;
+    e = num1.hi * num2.hi;
+    c = 2*((((a * c - e) + (a * d + b * c)) + b * d) + (num1.lo * num2.hi + num1.hi * num2.lo));
+    a = 2*e + c;
+    return dd_new(a, c + (2*e - a));
+}
+
+inline DoubleDouble dd_mul_ui(DoubleDouble num1, unsigned int num2) {
+    return dd_mul(num1, dd_new(num2, 0));
+}
+
+inline DoubleDouble dd_sub(DoubleDouble num1, DoubleDouble num2) {
+    double a, b, c, d, e, f, g;
+    g = num1.lo - num2.lo;
+    f = num1.lo - g;
+    e = num1.hi - num2.hi;
+    d = num1.hi - e;
+    d = ((num1.hi - (d + e)) + (d - num2.hi)) + g;
+    b = e + d;
+    c = (d + (e - b)) + ((num1.lo - (f + g)) + (f - num2.lo));
+    a = b + c;
+    return dd_new(a, c + (b - a));
+}
+
+inline DoubleDouble dd_sqr(DoubleDouble num) {
+    double a, b, c;
+    a = 0x08000001 * num.hi;
+    a += num.hi - a;
+    b = num.hi - a;
+    c = num.hi * num.hi;
+    b = ((((a * a - c) + a * b * 2) + b * b) + num.hi * num.lo * 2) + num.lo * num.lo;
+    a = b + c;
+    return dd_new(a, b + (c - a));
+}

mandel.c

@@ -4,20 +4,24 @@
 #include <stdlib.h>
 #include <stdbool.h>
 
-int main() {
+int main(int argc, char **argv) {
+    if (argc != 6) {
+        fprintf(stderr, "usage: %s width height centerx centery magnification", argv[0]);
+        return 1;
+    }
     long double eps = 1e-17;
     long double Q1LOG2 = 1.44269504088896340735992468100189213742664595415299;
     long double LOG2 = 0.69314718055994530941723212145817656807550013436026;
-    int width = 640;
+    int width = atoi(argv[1]);
-    int height = 480;
+    int height = atoi(argv[2]);
     char* image = malloc(width*height*3);
     int x, y;
-    long double centerx = -0.743643887037158704752191506114774;
+    long double centerx = strtold(argv[3], NULL);
-    long double centery = 0.131825904205311970493132056385139;
+    long double centery = strtold(argv[4], NULL);
     long double bailout = 128;
     long double logLogBailout = log(log(bailout));
     int foundperiods = 0;
-    long double magn = 1e18;
+    long double magn = strtold(argv[5], NULL);
     long maxiter = width * sqrt(magn);
     long double x0d = 4 / magn / width;
     long double x2 = -2 / magn + centerx;

mandel_dd.c

@@ -0,0 +1,139 @@
+#include "colors.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdbool.h>
+#include <math.h>
+#include <gmp.h>
+#include "doubledouble.h"
+
+int main(int argc, char **argv) {
+    if (argc != 6) {
+        fprintf(stderr, "usage: %s width height centerx centery magnification", argv[0]);
+        return 1;
+    }
+    DoubleDouble temp1;
+    double eps = 1e-17;
+    double Q1LOG2 = 1.44269504088896340735992468100189213742664595415299;
+    double LOG2 = 0.69314718055994530941723212145817656807550013436026;
+    unsigned int width = atoi(argv[1]);
+    unsigned int height = atoi(argv[2]);
+    char* image = malloc(width*height*3);
+    unsigned int x, y;
+    DoubleDouble centerx, centery;
+    centerx = dd_new(-0.7436438870371587, -3.628952515063387E-17);
+    centery = dd_new(0.13182590420531198, -1.2892807754956678E-17);
+    double bailout = 4; // the distance must not be greater than 2 (4 = 2*2)
+    double logLogBailout = log(log(bailout));
+    DoubleDouble magn = dd_new(strtod(argv[5], NULL), 0);
+    // maxiter = width * sqrt(magn);
+    temp1 = dd_sqrt(magn);
+    unsigned long maxiter = width * dd_get_ui(temp1);
+    DoubleDouble x2, y2, x0d, y1d;
+    // x0d = 4 / magn / width;
+    x0d = dd_ui_div(4, magn);
+    x0d = dd_div_ui(x0d, width);
+    // x2 = -2 / magn + centerx;
+    x2 = dd_si_div(-2, magn);
+    x2 = dd_add(x2, centerx);
+    // y1d = -4 / magn / width;
+    y1d = dd_si_div(-4, magn);
+    y1d = dd_div_ui(y1d, width);
+    // y2 = 2 / magn * height / width + centery;
+    y2 = dd_ui_div(2, magn);
+    temp1 = dd_new(height, 0);
+    temp1 = dd_div_ui(temp1, width);
+    y2 = dd_mul(y2, temp1);
+    y2 = dd_add(y2, centery);
+    DoubleDouble px, py, zx, zy, xx, yy;
+    unsigned int idx;
+    unsigned int imgidx = 0;
+    double hx, hy, d;
+    int foundperiods = 0;
+    /*fprintf(stderr, "centerx: %E\n", centerx.hi);
+    fprintf(stderr, "centery: %E\n", centery.hi);
+    fprintf(stderr, "magn:    %E\n", magn.hi);
+    fprintf(stderr, "x2:      %E\n", x2.hi);
+    fprintf(stderr, "y2:      %E\n", y2.hi);
+    fprintf(stderr, "x0d:     %E\n", x0d.hi);
+    fprintf(stderr, "y1d:     %E\n", y1d.hi);*/
+    for (y = 0; y < height; y++) {
+        for (x = 0; x < width; x++) {
+            fprintf(stderr, "\r%f %%", (float)imgidx/(width*height*3)*100);
+            //px = x*x0d + x2;
+            px = dd_mul_ui(x0d, x);
+            px = dd_add(px, x2);
+            //py = y*y1d + y2;
+            py = dd_mul_ui(y1d, y);
+            py = dd_add(py, y2);
+            // no Main bulb or Cardoid check to be faster
+            zx = dd_new(px.hi, px.lo);
+            zy = dd_new(py.hi, py.lo);
+            unsigned long i;
+            bool inside = true;
+            int check = 3;
+            int whenupdate = 10;
+            hx = 0;
+            hy = 0;
+            //for (i = 1; i <= maxiter; i++) {
+            for (i = 1; i <= 50000; i++) {
+                //xx = zx * zx;
+                xx = dd_sqr(zx);
+                //yy = zy * zy;
+                yy = dd_sqr(zy);
+                //if (xx + yy > bailout) {
+                if (xx.hi + yy.hi > bailout) {
+                    inside = false;
+                    break;
+                }
+                // iterate
+                //zy = 2 * zx * zy + py;
+                //zx = dd_mul_ui(zx, 2);
+                //zy = dd_mul(zx, zy);
+                zy = dd_add(dd_mul2(zx, zy), py);
+                //zx = xx - yy + px;
+                zx = dd_add(dd_sub(xx, yy), px);
+
+                // period checking
+                d = zx.hi - hx;
+                if (d > 0.0 ? d < eps : d > -eps) {
+                    d = zy.hi - hy;
+                    if (d > 0.0 ? d < eps : d > -eps) {
+                        // Period found.
+                        foundperiods++;
+                        break;
+                    }
+                }
+                if ((i & check) == 0) {
+                    if (--whenupdate == 0) {
+                        whenupdate = 10;
+                        check <<= 1;
+                        check++;
+                    }
+                    // period = 0;
+                    hx = zx.hi;
+                    hy = zy.hi;
+                }
+            }
+
+            if (inside) {
+                image[imgidx++] = 0;
+                image[imgidx++] = 0;
+                image[imgidx++] = 0;
+            } else {
+                double r = sqrt(zx.hi*zx.hi + zy.hi*zy.hi);
+                double c = i - 1.28 + (logLogBailout - log(log(r))) * Q1LOG2;
+                idx = fmod((log(c/64+1)/LOG2+0.45), 1)*GRADIENTLENGTH + 0.5;
+                image[imgidx++] = colors[idx][0];
+                image[imgidx++] = colors[idx][1];
+                image[imgidx++] = colors[idx][2];
+            }
+        }
+    }
+    // write out image
+    printf("P6 %d %d 255\n", width, height);
+    fwrite(image, 1, width*height*3, stdout);
+    //fprintf(stderr, "\r%d found periods", foundperiods);
+    fprintf(stderr, "\n");
+
+    return 0;
+}

run.sh

@@ -1,8 +1,8 @@
 #!/bin/sh
 
 start="1.0"
-end="1.0e+32"
+end="1.0e+30"
-frames="100"
+frames="14400"
 
 magn=`python -c "from math import pow; fac = pow($end / $start, 1.0 / ($frames - 1)); print \"\\n\".join([ str($start * pow(fac, i-1)) for i in range(1, $frames + 1)])"`
 
@@ -11,8 +11,15 @@ centery="0.131825904205311970493132056385139"
 
 i=0
 for mag in $magn; do
-	echo $mag
+	pngname=`printf "out_%05d_%s.png" $i $mag`
-	fname=`printf "out_%04d_%s.ppm" $i $mag`
+	if [ -s "$pngname" ]; then
-	./a.out 320 240 $centerx $centery $mag > "$fname"
+		i=$((i+1))
+		continue
+	fi
+	echo $i $mag
+	fname=`printf "out_%05d_%s.ppm" $i $mag`
+	/usr/bin/time -f "%e s" ./mandel_dd 320 240 $centerx $centery $mag > "$fname"
+	convert "$fname" -format png "$pngname"
+	rm "$fname"
 	i=$((i+1))
 done