Implementing series approximation/perturbation theory for WebGL2 based Mandelbrot viewer application

I have started building a mandelbrot viewer application with WebGL2 and JavaScript and am trying to implement series approximation at a basic level. I'm currently getting weird/distored/incomplete images depending on the number of iterations and the reference point I'm using (currently just the center of the viewport in the complex plane).

I don't have a heavy math background so it's entirely possible (apparent) that I've screwed up my calculations somewhere. I know there are separate concerns with choosing reference points, the number of iterations that can be done with the coefficients without distorting the image, etc., but I'm just trying to see the basic concept working. Currently I am only skipping 1 iteration with a reference point near the origin (-0.5, 0.0), and I don't think the result is supposed to be this distorted.

If someone is able to see and explain to me where the logic in the code has gone wrong to produce this result I would greatly appreciate it. I'll try to include all of the relevant code, including the fragment shader and JavaScript code so you should have all the information needed.

fragment.glsl

#version 300 es
precision highp float;

in vec2 uv;

uniform sampler2D palette;
uniform sampler2D orbit;

layout(std140) uniform blockData {
  vec2 center;
  float scale;
  float aspect;
  float radius;
  int iterations;
  int width;
  mat4x2 coefficients;
  vec2 zStart;
};

out vec4 fragColor;

ivec2 uvFromIndex(int i) {
  int col = i % width;
  int row = i / width;
  return ivec2(col, row);
}

vec2 cxMul(vec2 a, vec2 b) {
  return vec2(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}

vec2 cxPow(vec2 a, float n) {
  float angle = atan(a.y, a.x);
  float r = length(a);
  float real = pow(r, n) * cos(n * angle);
  float im = pow(r, n) * sin(n * angle);
  return vec2(real, im);
}

float iterate() {
  vec2 dc = scale * vec2(aspect * (uv.x - 0.5), (uv.y - 0.5)) + center;
  vec2 z = zStart;
  vec2 dz = cxMul(coefficients[0], dc) + cxMul(coefficients[1], cxPow(dc, 2.)) + cxMul(coefficients[2], cxPow(dc, 3.)) + cxMul(coefficients[3], cxPow(dc, 4.));

  for(int i = 0; i < iterations; i++) {
    dz = cxMul(2. * z + dz, dz) + dc;
    if(length(dz) > radius)
      return float(i) + 1. - log(log(length(dz))) / log(2.);
    z = texelFetch(orbit, uvFromIndex(i), 0).rg;
  }
  return 0.;
}

void main() {
  float n = iterate();
  fragColor = texture(palette, vec2(n / float(iterations), 0.5));
}

app.ts

import PicoGL, { DrawCall } from "picogl";

import { Complex } from "./complex";
import palettes from "./palettes";
import shaders from "./shaders";
import UniformBufferWrapper from "./UniformBufferWrapper";
import { getReferenceOrbit, getSeriesCoefficients } from "./utils";

const digitCodeRegex = /^Digit(\d)$/;

const defaultValues = {
  center: [-0.5, 0] as [number, number],
  scale: 3,
  radius: 2,
  iterations: 256,
  width: 16,
};

export default (function main() {
  const canvas = document.querySelector("canvas") as HTMLCanvasElement;

  const app = PicoGL.createApp(canvas).resize(
    canvas.clientWidth,
    canvas.clientHeight
  );

  const quadPositions = app.createVertexBuffer(
    PicoGL.FLOAT,
    2,
    new Float32Array([
      -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0,
    ])
  );

  const quadArray = app
    .createVertexArray()
    .vertexAttributeBuffer(0, quadPositions);

  const uniformBuffer = new UniformBufferWrapper(
    app,
    [
      PicoGL.FLOAT_VEC2,
      PicoGL.FLOAT,
      PicoGL.FLOAT,
      PicoGL.FLOAT,
      PicoGL.INT,
      PicoGL.INT,
      PicoGL.FLOAT_MAT4x2,
      PicoGL.FLOAT_VEC2,
    ],
    [
      "center",
      "scale",
      "aspect",
      "radius",
      "iterations",
      "width",
      "coefficients",
      "zStart",
    ]
  );

  let drawCall: DrawCall;
  let image: HTMLImageElement;

  let z = new Complex(0);
  let coefficients = getSeriesCoefficients(
    z,
    new Complex(...defaultValues.center),
    // defaultValues.iterations
    1
  );
  console.log("coefficients:", coefficients);
  let zStart = z.toNumberArray();

  uniformBuffer
    .set("center", new Float32Array(defaultValues.center))
    .set("scale", new Float32Array([defaultValues.scale]))
    .set("aspect", new Float32Array([canvas.width / canvas.height]))
    .set("radius", new Float32Array([defaultValues.radius]))
    .set("iterations", new Int32Array([defaultValues.iterations]))
    .set("width", new Int32Array([defaultValues.width]))
    .set("coefficients", new Float32Array(coefficients))
    .set("zStart", new Float32Array(zStart))
    .update();

  let orbit = getReferenceOrbit(
    z,
    new Complex(...defaultValues.center),
    defaultValues.iterations
  );

  const orbitTexture = app.createTexture2D(
    orbit,
    defaultValues.width,
    defaultValues.width,
    {
      internalFormat: PicoGL.RG32F,
      wrapS: PicoGL.CLAMP_TO_EDGE,
      wrapT: PicoGL.CLAMP_TO_EDGE,
      magFilter: PicoGL.NEAREST,
      minFilter: PicoGL.NEAREST,
    }
  );

  app
    .createPrograms([shaders.vertex, shaders.fragment])
    .then(function ([program]) {
      image = new Image();

      image.onload = function () {
        const paletteTexture = app.createTexture2D(image, {
          magFilter: PicoGL.NEAREST,
          minFilter: PicoGL.NEAREST,
        });

        drawCall = app
          .createDrawCall(program, quadArray)
          .texture("palette", paletteTexture)
          .texture("orbit", orbitTexture)
          .uniformBlock("blockData", uniformBuffer.getBuffer());

        drawCall.draw();
      };

      image.src = palettes[0];
    });

  window.onresize = function () {
    app.resize(canvas.clientWidth, canvas.clientHeight);

    uniformBuffer
      .set("aspect", new Float32Array([canvas.width / canvas.height]))
      .update();

    drawCall.draw();
  };

  canvas.onclick = function (e) {
    updateScreenBounds(e, 0.8);
  };

  canvas.oncontextmenu = function (e) {
    e.preventDefault();
    updateScreenBounds(e, 1.25);
  };

  let result: RegExpExecArray | null;

  window.onkeydown = function (e) {
    if ((result = digitCodeRegex.exec(e.code)) !== null) {
      const index = +result[1];
      image.src = palettes[index];
    }
  };

  function updateScreenBounds(position: MousePosition, scaleFactor: number) {
    const [center, scale, aspect] = uniformBuffer.get(
      "center",
      "scale",
      "aspect"
    );

    const mouse = new Float32Array([
      aspect[0] * scale[0] * (position.x / canvas.width - 0.5),
      scale[0] * ((canvas.height - position.y) / canvas.height - 0.5),
    ]);

    const newCenter = [
      center[0] + mouse[0] * 0.5,
      center[1] + mouse[1] * 0.5,
    ] as [number, number];

    z = new Complex(0);
    coefficients = getSeriesCoefficients(
      z,
      new Complex(...newCenter),
      // defaultValues.iterations
      1
    );
    console.log("coefficients:", coefficients);
    zStart = z.toNumberArray();

    uniformBuffer
      .set("center", new Float32Array(newCenter))
      .set("scale", new Float32Array([scale[0] * scaleFactor]))
      .set("coefficients", new Float32Array(coefficients))
      .set("zStart", new Float32Array(zStart))
      .update();

    orbit = getReferenceOrbit(
      z,
      new Complex(...newCenter),
      defaultValues.iterations
    );

    orbitTexture.data(orbit);

    drawCall.draw();
  }
})();

type MousePosition = { x: number; y: number };

utils.ts

import { Complex } from "./complex";

export function getReferenceOrbit(
  z: Complex,
  c: Complex,
  iterations: number,
  radius = 2
) {
  const orbit = new Float32Array(iterations * 2);
  for (let i = 0; i < iterations * 2; i += 2) {
    iterate(z, c);
    if (z.abs().toNumber() > radius)
      console.warn(`${c} is not a valid reference point!`);
    orbit[i] = z.re.toNumber();
    orbit[i + 1] = z.im.toNumber();
  }
  return orbit;
}

export function getSeriesCoefficients(
  z: Complex,
  c: Complex,
  iterations: number,
  radius = 2
) {
  const coefficients = [
    new Complex(1),
    new Complex(0),
    new Complex(0),
    new Complex(0),
  ];
  for (let i = 0; i < iterations; i++) {
    iterate(z, c);
    if (z.abs().toNumber() > radius)
      console.warn(`${c} is not a valid reference point!`);
    const [A, B, C, D] = [...coefficients];
    coefficients[0] = A.mul(z.mul(2)).add(1);
    coefficients[1] = B.mul(z.mul(2)).add(A.mul(A));
    coefficients[2] = C.mul(z.mul(2)).add(B.mul(A.mul(2)));
    coefficients[3] = D.mul(z.mul(2))
      .add(C.mul(A.mul(2)))
      .add(B.mul(B));
  }
  return coefficients.flatMap((x) => x.toNumberArray());
}

function iterate(z: Complex, c: Complex) {
  const x = z.mul(z).add(c);
  z.re = x.re;
  z.im = x.im;
}

complex.ts

import Decimal from "decimal.js";

export namespace Complex {
  export type Constructor = typeof Complex;
  export type Value = [Decimal.Value, Decimal.Value?] | [Complex];
}

export class Complex {
  re: Decimal;
  im: Decimal;

  constructor(...n: Complex.Value) {
    if (n[0] instanceof Complex) {
      this.re = n[0].re;
      this.im = n[0].im;
    } else {
      this.re = new Decimal(n[0]);
      this.im = new Decimal(n[1] ?? 0);
    }
  }

  add(...n: Complex.Value) {
    const c = new Complex(...n);
    return new Complex(this.re.add(c.re), this.im.add(c.im));
  }

  sub(...n: Complex.Value) {
    const c = new Complex(...n);
    return new Complex(this.re.sub(c.re), this.im.sub(c.im));
  }

  mul(...n: Complex.Value) {
    const c = new Complex(...n);
    return new Complex(
      this.re.mul(c.re).sub(this.im.mul(c.im)),
      this.re.mul(c.im).add(this.im.mul(c.re))
    );
  }

  div(...n: Complex.Value) {
    const c = new Complex(...n);
    return new Complex(
      this.re
        .mul(c.re)
        .add(this.im.mul(c.im))
        .div(c.re.mul(c.re).add(c.im.mul(c.im))),
      this.im
        .mul(c.re)
        .sub(this.re.mul(c.im))
        .div(c.re.mul(c.re).add(c.im.mul(c.im)))
    );
  }

  abs() {
    return this.re.mul(this.re).add(this.im.mul(this.im)).sqrt();
  }

  toNumberArray() {
    return [this.re.toNumber(), this.im.toNumber()];
  }
}

I solved my own problem. The coefficients and the current "xn" (z) value were out of sync by 1 iteration. Re-reading this paper helped me realize that I was skipping an iteration of x which throws off the delta n calculations. I now have a better idea of how to apply the same logic for other reference points besides the center of the screen (should just have to add the distance of x0 from center when computing delta 0) and can begin exploring optimizations like choosing good reference points and detecting glitches, calculating how many times to iterate with the coefficients, etc.

The fractal looks normal with the x sub 0 coefficients plugged in at broad zoom levels, and when applied beyond that becomes increasingly accurate at deeper zoom levels like it's supposed to. Relevant fixed code:

fragment.glsl

// ...
float iterate() {
  vec2 d0 = scale * vec2(aspect * (uv.x - 0.5), (uv.y - 0.5)); // + center - x0;
  vec2 dn = cxMul(coefficients[0], d0) + cxMul(coefficients[1], cxPow(d0, 2.)) + cxMul(coefficients[2], cxPow(d0, 3.)) + cxMul(coefficients[3], cxPow(d0, 4.));
  vec2 xn;

  for(int i = 0; i < iterations; i++) {
    xn = texelFetch(orbit, uvFromIndex(i), 0).xy;
    dn = cxMul(2. * xn + dn, dn) + d0;
    if(length(dn) > radius)
      return float(i) + 1. - log(log(length(dn))) / log(2.);
  }
  return 0.;
}
// ...

app.ts

// ...
const defaultValues = {
  center: [-0.5, 0] as [number, number],
  scale: 2.5,
  radius: 2,
  iterations: 100,
  width: 10,
  coefficients: [1, 0, 0, 0, 0, 0, 0, 0],
};

(function main() {
// ...