I want to do some graphics on a PIC - so I need some "real" datatypes. But I'm not sure how many bits of fixed-point resolution I need. It's easy enough to code up fixed-point routines for any multiple of 8 bits, but I'd rather just do it once. So, I wrote some templated C++ code to simulate fixed-point types. Here's a 640x480 Mandelbrot set in 3.5 8-bit fixed point math (generated on a PC):

Not great. I knew I'd need more than 8-bits, though :-) Here's the same in 3.13 16-bits:

Much better. Of course, if you zoomed in, you'd need more bits of precision. You can test 24, 32, 40, 48, 56, or 64-bits with my trivial C++ class below. I started testing with fractals, because the algorithm is so simple: in my teens, I got it published as an Apple II BASIC one-liner in Nibble magazine. Anyway, templated over the datatypes, it looks like this:

template <typename iter_t, typename real_t>
iter_t
mandelbrot_test(real_t c, real_t d, iter_t max_iter, real_t max_mag)
{
  real_t a, b;
  a = real_t(0.);
  b = real_t(0.);
  iter_t i = 0;
  while (i < max_iter){
    real_t aa = a * a;
    real_t bb = b * b;
    if (aa + bb > max_mag){
      break;
    }
    b = real_t(2.) * a * b + d;
    a = aa - bb + c;
    i++;
  }
  return i;
}

This is a direct translation of the floating-point algorithm; it could be improved some to avoid issues with fixed-point overflows.

So far, I've only implemented addition, subtraction, and multiplication, which is sufficient for these fractals:

template <int width, int scale>
class fixed_pt
{
public:
  fixed_pt(int64_t val = 0)
    : w(width),
      s(scale),
      v(val)
  {
  }
  fixed_pt(double x)
    : w(width),
      s(scale)
  {
    v = int64_t((x * (int64_t(1) << w))) >> s;
    if (x >= 0){
      v &= 0xffffffffffffffffull >> (64-w);
    } else {
      v |= 0xffffffffffffffffull << (64-w);
    }
  }
  operator double()
  {
    return v / double(int64_t(1) << (w-s));
  }
  friend fixed_pt operator+ (fixed_pt<width, scale> a, fixed_pt<width, scale> b)
  {
    return fixed_pt<width, scale>(double(a) + double(b));
  }
  friend fixed_pt operator- (fixed_pt<width, scale> a, fixed_pt<width, scale> b)
  {
    return fixed_pt<width, scale>(double(a) - double(b));
  }
  friend fixed_pt operator* (fixed_pt<width, scale> a, fixed_pt<width, scale> b)
  {
    return fixed_pt<width, scale>(double(a) * double(b));
  }
private:
  int64_t v;
  int w;
  int s;
};

Yes, it's a hack - meant in the worst pejorative sense - but it was easy to implement, and will give me some quick results without having to dive into embedded issues yet. Once I decide on the number of bits I need, I can start coding up some routines - I know I'll need at least these operations.

I'll also need division if I want to get ray-tracing going. Square roots are probably also necessary, and those can be implemented in terms of the other operations - they might not be fast, but who cares?

I didn't implement any colors here: I'm saving that for the embedded code :-)

Now, I can try some ray-tracing. That will take a little more time...

EDIT

So, I really should have done this first, but I was pretty sure of the answer. I just compiled the most bare-bones calculation of the above fractal in C using Microchips XC8 compiler, a PIC16 target, and "float" datatypes: 1565 instructions = 2738.8 bytes. Yeah, you can't just take the easy way.