Support complex powers of z

config.h

@@ -6,13 +6,13 @@
 #define MAXITER 1000
 
 // Pretty Julia set zoom
-#define SS 9
+//#define SS 9
-#define JULIA
+//#define JULIA
-#define JULIA_C_RE 0.285
+//#define JULIA_C_RE 0.285
-#define JULIA_C_IM 0.010
+//#define JULIA_C_IM 0.010
-#define CENTRE_RE -0.47353731
+//#define CENTRE_RE -0.47353731
-#define CENTRE_IM -0.18894516
+//#define CENTRE_IM -0.18894516
-#define ZOOMRATE 0.01
+//#define ZOOMRATE 0.01
 
 // Cool zoomy swirly animated Julia set
 //#define SS 9
@@ -43,3 +43,8 @@
 //#define JULIA_DTHETA -0.00002
 //#define POWRATE 0.0001
 //#define ZOOMRATE 0.0003
+
+// Complex power Mandelbrot
+#define SS 4
+#define CPOW_DTHETA 0.01
+#define MAXMAG 100.0

frag_shader.glsl

@@ -9,6 +9,10 @@
 
 #define NCOLS 8
 
+#ifndef CPOW_DTHETA
+#define MAXMAG 2.0
+#endif
+
 layout(location = 0) out vec4 out_colour;
 
 layout(push_constant) uniform Constants {
@@ -17,10 +21,13 @@ layout(push_constant) uniform Constants {
 	vec2 julia;
 #endif
 	vec2 centre;
-	float scale;
+#ifdef CPOW_DTHETA
+	vec2 cpow;
+#endif
 #ifdef POWRATE
 	float zpow;
 #endif
+	float scale;
 } params;
 
 vec3 incol = vec3(0.000000, 0.014444, 0.027321);
@@ -61,13 +68,37 @@ vec2 maptoz(vec2 xy)
 	return params.scale * z + params.centre;
 }
 
-vec2 zpow(vec2 z, float p) {
+#ifdef CPOW_DTHETA
+vec2 cmul(vec2 a, vec2 b)
+{
+	return vec2(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
+}
+
+vec2 cexp(vec2 z)
+{
+	return exp(z.x) * vec2(cos(z.y), sin(z.y));
+}
+
+vec2 cln(vec2 z)
+{
+	return vec2(log(length(z)), atan(z.y, z.x));
+}
+
+vec2 cpow(vec2 z, vec2 p)
+{
+	if (z.x == 0.0 && z.y == 0.0) return vec2(0.0, 0.0);
+	return cexp(cmul(p, cln(z)));
+}
+#else
+vec2 zpow(vec2 z, float p)
+{
     float r = length(z);
     float theta = atan(z.y, z.x);
     float mag = pow(r, p);
     float phi = p * theta;
     return mag * vec2(cos(phi), sin(phi));
 }
+#endif
 
 int fractal(vec2 c) {
 #ifdef JULIA
@@ -77,12 +108,14 @@ int fractal(vec2 c) {
 	vec2 z = vec2(0.0, 0.0);
 #endif
 	for (int i = 0; i < MAXITER; ++i) {
-#ifdef POWRATE
+#ifdef CPOW_DTHETA
+		z = cpow(z, params.cpow) + c;
+#elif defined(POWRATE)
 		z = zpow(z, params.zpow) + c;
 #else
 		z = vec2(z.x * z.x - z.y * z.y, 2.0 * z.x * z.y) + c;
 #endif
-		if (length(z) > 2.0)
+		if (length(z) > MAXMAG)
 			return i;
 	}
 	return -1;

main.c

@@ -10,6 +10,7 @@
 #include <GLFW/glfw3.h>
 #include <assert.h>
 #include <complex.h>
+#include <math.h>
 #include <stdbool.h>
 #include <stdio.h>
 #include <string.h>
@@ -26,6 +27,10 @@
 
 #define MAX_SHADER_SIZE (8 * 1024)
 
+#if defined(CPOW_DTHETA) && defined(POWRATE)
+#error "Cannot define CPOW_DTHETA and POWRATE simulataneously"
+#endif
+
 typedef struct {
 	float width, height;
 #ifdef JULIA
@@ -36,14 +41,23 @@ typedef struct {
 	struct {
 		float re, im;
 	} centre;
-	float scale;
+#ifdef CPOW_DTHETA
+	struct {
+		float re, im;
+	} cpow;
+#endif
 #ifdef POWRATE
 	float zpow;
 #endif
+	float scale;
 } params_t;
 
 static const char *layers[] = { "VK_LAYER_KHRONOS_validation" };
-static const char *swapchain_ext_name = VK_KHR_SWAPCHAIN_EXTENSION_NAME;
+static const char *extensions[] = {
+	VK_KHR_SWAPCHAIN_EXTENSION_NAME,
+};
+
+#define NEXTENSIONS (sizeof(extensions) / sizeof(extensions[0]))
 
 // Initial parameters
 static params_t params = {
@@ -67,10 +81,13 @@ static params_t params = {
 #endif
 	},
 #endif
-	.scale = 1.0,
+#ifdef CPOW_DTHETA
+	.cpow = { .re = 2.0, .im = 0.0 },
+#endif
 #ifdef POWRATE
 	.zpow = 2.0,
 #endif
+	.scale = 1.0,
 };
 
 static void
@@ -192,6 +209,28 @@ int main(void)
 	}
 	assert(found_queue_family);
 
+	// Check that the physical device supports all extensions
+	uint32_t ext_count = EXT_PROP_BUFFER_SIZE;
+	VkExtensionProperties ext_props[EXT_PROP_BUFFER_SIZE];
+	result = vkEnumerateDeviceExtensionProperties(
+	    physical_dev, NULL, &ext_count, ext_props);
+	assert(result == VK_SUCCESS);
+	bool supported[NEXTENSIONS] = { 0 };
+	for (unsigned i = 0; i < ext_count; ++i) {
+		for (unsigned j = 0; j < NEXTENSIONS; ++j) {
+			if (strcmp(extensions[j], ext_props[i].extensionName) == 0) {
+				supported[j] = true;
+				break;
+			}
+		}
+	}
+	for (unsigned j = 0; j < NEXTENSIONS; ++j) {
+		if (!supported[j]) {
+			fprintf(stderr, "Unsupported: %s\n", extensions[j]);
+			assert(false);
+		}
+	}
+
 	// Create the logical device and get the queue handle.
 	float queue_priority = 1.0;
 	const VkDeviceQueueCreateInfo queue_create_info = {
@@ -204,8 +243,8 @@ int main(void)
 		.sType = VK_STRUCTURE_TYPE_DEVICE_CREATE_INFO,
 		.pQueueCreateInfos = &queue_create_info,
 		.queueCreateInfoCount = 1,
-		.enabledExtensionCount = 1,
+		.enabledExtensionCount = NEXTENSIONS,
-		.ppEnabledExtensionNames = &swapchain_ext_name,
+		.ppEnabledExtensionNames = extensions,
 	};
 	VkDevice dev;
 	result = vkCreateDevice(physical_dev, &dev_create_info, NULL, &dev);
@@ -213,21 +252,6 @@ int main(void)
 	VkQueue queue;
 	vkGetDeviceQueue(dev, queue_family, 0, &queue);
 
-	// Check that the physical device has swap chain support
-	uint32_t ext_count = EXT_PROP_BUFFER_SIZE;
-	VkExtensionProperties ext_props[EXT_PROP_BUFFER_SIZE];
-	result = vkEnumerateDeviceExtensionProperties(
-	    physical_dev, NULL, &ext_count, ext_props);
-	assert(result == VK_SUCCESS);
-	bool swapchain_support = false;
-	for (unsigned i = 0; i < ext_count; ++i) {
-		if (strcmp(swapchain_ext_name, ext_props[i].extensionName) == 0) {
-			swapchain_support = true;
-			break;
-		}
-	}
-	assert(swapchain_support);
-
 	// Select a surface format, preferring R8G8B8A8_SRGB with
 	// SRGB_NONLINEAR colour space.
 	uint32_t surface_fmt_count = SURFACE_FMT_BUFFER_SIZE;
@@ -534,8 +558,10 @@ int main(void)
 	params.width = (float)swapchain_extent.width;
 	params.height = (float)swapchain_extent.height;
 #ifdef JULIA_DTHETA
-	const double dtheta = JULIA_DTHETA;
+	const double complex julia_rotz = cexp(I * JULIA_DTHETA);
-	const double complex rotz = cexp(I * dtheta);
+#endif
+#ifdef CPOW_DTHETA
+	const double complex cpow_rotz = cexp(I * CPOW_DTHETA);
 #endif
 
 	while (!glfwWindowShouldClose(window)) {
@@ -617,7 +643,7 @@ int main(void)
 
 #ifdef JULIA_DTHETA
 		double complex julia = params.julia.re + I * params.julia.im;
-		julia *= rotz;
+		julia *= julia_rotz;
 		params.julia.re = (float)creal(julia);
 		params.julia.im = (float)cimag(julia);
 #endif
@@ -626,6 +652,13 @@ int main(void)
 		params.scale *= (1 - ZOOMRATE);
 #endif
 
+#ifdef CPOW_DTHETA
+		double complex cpow = params.cpow.re + I * params.cpow.im;
+		cpow *= cpow_rotz;
+		params.cpow.re = (float)creal(cpow);
+		params.cpow.im = (float)cimag(cpow);
+#endif
+
 #ifdef POWRATE
 		params.zpow += POWRATE;
 #endif