# Computing Julia Set on Accelerator, using SYCL DPC++

### Created : September 19, 2021

This week I was exploring fractals - specifically computing fractals using Julia Set quadratic equation. I considered using SYCL DPC++ for implementation, which I've been working on for few weeks now. I was interested in using image interface of DPC++, which just same as buffer interface creates an abstraction over raw memory address based data read/ write ( read pointer arithmetic ), while providing some image specific functionality for data interpretation, reading & writing. I find buffer & image interfaces easier to work with rather than raw pointers, allocated using Unified Shared Memory, due to their elegance in design; accessor based data read/ write on device/ host, providing implicit data movement/ synchronization/ visibility, which in-turn lets me implicitly define task graph dependencies. So it's good, nice, powerful abstraction.

I'll now spend some time explaining how to sequentially compute Julia Set fractal.

z = z ** 2 + c, where z, c ∈ ℂ && c is a constant

I sequentially apply aforementioned equation on each pixel of target fractal image as long as certain conditions satisfy. Let us assume, I'm working with an image of dimension 4 x 4. This is how each cell ( read pixel ) is addressed in 4 x 4 image. Notice the origin of reference frame in use. These indices can easily be generated using two nested for loops.

(uint, uint) A; for i in 0..4 { for j in 0..4 { A[i][j] = (i, j); } } But for computing fractal, I've to reinterpret indices in a way where origin of reference frame should move to center of image. And those points on 2D cartesian plane to be interpreted as complex numbers ( read a + ib, where a, b are components of 2D point along X & Y-axis respectively ). It can be easily generated using meshgrid method.

\$ python3 >>> w = np.linspace(-2, 1, 4); w // along X-axis array([-2., -1., 0., 1.]) >>> h = np.linspace(2, -1, 4); h // along Y-axis array([ 2., 1., 0., -1.]) >>> a, b = np.meshgrid(w, h) >>> a // visualise 4 X-axes placed along rows of matrix array([[-2., -1., 0., 1.], [-2., -1., 0., 1.], [-2., -1., 0., 1.], [-2., -1., 0., 1.]]) >>> b // visualise 4 Y-axes placed along columns of matrix array([[ 2., 2., 2., 2.], [ 1., 1., 1., 1.], [ 0., 0., 0., 0.], [-1., -1., -1., -1.]]) >>> a + b*1j // complex numbers, required for fractal computation array([[-2.+2.j, -1.+2.j, 0.+2.j, 1.+2.j], [-2.+1.j, -1.+1.j, 0.+1.j, 1.+1.j], [-2.+0.j, -1.+0.j, 0.+0.j, 1.+0.j], [-2.-1.j, -1.-1.j, 0.-1.j, 1.-1.j]]) Notice along X-axis points are in range [-2, 2), while along Y-axis they're in range [2, -2). For scaling them to S, I first take them to (-1, 1) interval & then multiply by S. Now as I've these complex numbers in matrix form properly scaled, I'd like to go to actual fractal computation part.

>>> W, H = 4, 4 // image dimensions >>> factor = ((W/2) ** 2 + (H/2) ** 2) ** 0.5; factor // √((w/2) ** 2 + (h/2) ** 2) 2.8284271247461903 >>> a / factor // all values in interval (-1, 1) array([[-0.70710678, -0.35355339, 0. , 0.35355339], [-0.70710678, -0.35355339, 0. , 0.35355339], [-0.70710678, -0.35355339, 0. , 0.35355339], [-0.70710678, -0.35355339, 0. , 0.35355339]]) >>> # do same for b i.e. `b / factor` >>> a / factor * S + (b/factor * S) * 1j // S = scale, here 2.f array([[-1.41421356+1.41421356j, -0.70710678+1.41421356j, 0.+1.41421356j, 0.70710678+1.41421356j], [-1.41421356+0.70710678j, -0.70710678+0.70710678j, 0.+0.70710678j, 0.70710678+0.70710678j], [-1.41421356+0.j , -0.70710678+0.j , 0. +0.j, 0.70710678+0.j ], [-1.41421356-0.70710678j, -0.70710678-0.70710678j, 0.-0.70710678j, 0.70710678-0.70710678j]])

Quadratic iteration based Julia Set computation can be encapsulated in following function. As long as z.real()2 + z.imag()2 < S2, I keep applying z = z2 + c. There is another bound on number of iterations that can be performed.

const uint MAX_ITR = 1 << 10; // 1024 const float S = 2.f; // scale func length(complex z) -> float { return (z.real() ** 2 + z.imag() ** 2) ** .5f; // z = a + ib; √(a**2 + b**2) } func quadratic_iteration(complex z, c) {} uint i = 0; for ; i < MAX_ITR; i++ { z = z ** 2 + c; // z ** 2 = (a+ib) * (a+ib) = (a**2-b**2) + 2ab*i if (length(z) > S) { break; } } return i; }

I choose S ( = scale ) such that S2 - S >= abs(c), where c is a constant complex number. Say, I choose to work with c = (-0.8+0.156j), in that case setting S = 2.f, satisfies aforementioned constraint.

>>> S = 2.; S // scale 2.0 >>> c = -0.8 + 0.156j; c // constant (-0.8+0.156j) >>> S ** 2 - S >= abs(c) True ✅

From quadratic_iteration(...), I obtain #-of evaluations it takes for z to determine whether it's in set with radius S or not. Depending upon #-of evaluations for certain z, I assign different color to respective image pixel. #-of iterations denote how many of it is required by z to get outside of Julia Set with radius S, which is alongside capped at maximum number of iterations ( = MAX_ITR ). After every evaluation of z = z2 + c, I check whether distance of updated z from center (0.+0j) is still less than radius of Julia Set defined on complex plane. This way I keep increasing z by (z2 + c) in each iteration & in some time z will reach some value in complex plane, which goes outside of set. I take that certain iteration number & use it to represent color.

Recently I discovered a nice way to assign beautiful colors to each pixel of image using cosine function. It simply requires me to evaluate color(t) = a + b * cos(2 * pi * (c * t + d)), where I can variate a, b, c, d ( 3-element floating point tuples ) to get different palettes. In this case t = quadratic_iteration(...) / MAX_ITR.
I recommend reading more about cosine palette based color generation technique here.

Today I'm using following method for colorizing each pixel of fractal.

[float; 4] colorize(uint itr) { [float; 3] d = {0.5, 0.5, 0.5}; [float; 3] e = {0.5, 0.5, 0.5}; [float; 3] f = {2.0, 1.0, 0.0}; [float; 3] g = {0.3, 0.2, 0.4}; // last 1.f dictates α channel value of RGBA image [float; 4] pixel_val = {d + e * cos(6.28318f * (f * ((float)itr / (float)(MAX_ITR)) + g)), 1.f}; return pixel_val; }

Now I've one sequential program which computes Julia Set Fractal of dimension 512x512, by evaluating quadratic equation z = z2 + c. Following is the computed fractal with c = (-0.8 + 0.156j). It's good time to move to parallel implementation of Julia Set computation on heterogeneous accelerators ( read CPU, GPGPU, FPGA etc. ). Being a good candidate for SIMD model, Julia Set computation fits well in data parallel programming domain, where each pixel of resulting image is processed independently & RGBA pixel value is written in respective memory location. The only challenge I see --- how to compute complex number (a + ib) located in pixel location (r, c) ?

In previous iteration, I was using meshgrid method to do so, this time I've to implement it manually. Each kernel invocation processes single ( unique ) pixel of image, so certain invocation has knowledge of which pixel it is responsible of processing. Let me assume certain kernel invocation knows, it has to process image pixel (r, c). Again my interest is to move reference frame to center of image.

r ∈ [0, H) && H = image height c ∈ [0, W) && W = image width

Notice each (x, y) point on 2D cartesian plane, is seperated by distance 1, so I can compute what should be (x, y) in current pixel location (r, c), with following formula.

w = h = 4 // image dimension func f(r, c) { step = 1.; // a + s * k return -(w/2) + step * c, (h/2) - step * r; } >>> assert f(0, 0) == (-2., 2.) ✔ // top left >>> assert f(0, 3) == (1., 2.) ✔ // top right >>> assert f(3, 0) == (-2., -1.) ✔ // bottom left >>> assert f(3, 3) == (1., -1.) ✔ // bottom right

Rest is same as I've already done while implementing sequential algorithm for fractal computation. I need to scale (x, y) point by factor S, so I rewrite aforementioned function as below, which I make use of for converting (r, c) image pixel to complex number (a+ib) which is to be operated on during fractal computation i.e. checking whether that complex number is part of Julia Set or not.

// check https://gist.github.com/itzmeanjan/0a13447c634afbc34f64d1f22751011d#file-julia_set-cpp-L16-L23 // for more info [float; 2] convert_to_complex_number(uint r, c) { float step = 1.; float factor = ((w/2) ** 2 + (h/2) ** 2) ** .5; // S = scale factor return {( ( -(w/2) + step * c ) / factor ) * S, ( ( (h/2) - step * r ) / factor ) * S}; // complex number (a+ib) as (a, b) }

I play around with complex constat c, with c = -0.123 + 0.745j, it generates following fractal. Here's a dendrite shaped fractal, with c = 0.f + 1.0j. With c = -0.75 + -0.2j, beautiful galaxy shaped fractal is generated. Good thing about cosine based color palette is just changing to following starts generating fractals with very different color.

{ // check https://gist.github.com/itzmeanjan/0a13447c634afbc34f64d1f22751011d#file-julia_set-cpp-L38-L41 // for more info [float; 3] d, e, f = /* unchanged */ ; [float; 3] g = {0.5, 0.2, 0.25}; // updated // ... } I keep whole DPC++ Julia Set computation implementation here for future reference. In coming days, I plan to compute Julia Set with sine/ cosine function, as this time I did with quadratic function. I'm also interested in Newton Fractal, which I'll consider implementing next week in DPC++.

Have a great time !