Graeme Winter

Mandelbrot set, with integers

Started wondering if I can calculate fractals on an rp2040 but to make it efficient feels like a little more of an ask - certainly in the absence of a floating point unit, writing something quick will be a job of work. So, will be assembly then right?

How to assemble a Mandelbrot set from scratch?

Fixed point arithmetic

Untimately computers just work on 0, 1. More interesting numbers come from how you interpret these e.g. IEEE754 floating point values, where you have 23 bits for the fractional value and 8 bits for exponent. Since I know with a fractal that I won’t be working with values greater than about 16 (certainly if my “escape” is 4) then I can instead use fixed point where I do everything with signed integers, except the smallest bit is not equal to 1 but instead say 1>>24. Therefore, will use signed 7.24 format, which can be ignored except for the need to right shift by 24 bits on multiplication.

Code

Working Manedlbrot set calculation, writes out a stream of unsigned shorts very inefficiently to stdout but requires no storage (useful for assembly programming…) - however this is the C initial implementation:

#include <unistd.h>

// fixed point integer arithmetic formulation
//
// use s7.24 formulation for bits => have enough
// bits to give ~ 7 decimal places of accuracy
// i.e. similar to floats

#define mul(a, b) (int)(((long long)a) * ((long long)b) >> 24)

unsigned short iter(int cr, int ci) {
  unsigned short count = 0;
  int zr = 0;
  int zi = 0;
  int tmp;

  while (count < 4096) {
    int zr2 = mul(zr, zr);
    int zi2 = mul(zi, zi);
    if ((zr2 + zi2) > (4 << 24)) {
      break;
    }
    count++;
    tmp = zr;
    zr = zr2 - zi2 + cr;
    zi = 2 * mul(tmp, zi) + ci;
  }
  return count;
}

int main() {
  for (int j = 0; j < 1280; j++) {
    for (int i = 0; i < 1280; i++) {
      int cr = -(2 << 24) + 0x8000 * i + 0x4000;
      int ci = -(5 << 22) + 0x8000 * j + 0x4000;
      unsigned short c = iter(cr, ci);
      write(1, (void *)&c, 2);
    }
  }
  return 0;
}

Works, makes something which looks a lot like a Mandelbot set. Plotted as log2(counts) it certainly looks the part:

Next: do the actual assembly port and figure out how to extract the results…