I just started implementing a simple mandelbrot zoom app based on the Processing example, and I think I'm hitting float precision limits.  It feels like I'm hitting the limits too early though, with numbers around e-7.

Can anybody comment on where in the app I'm losing precision, and what I can do to improve it?

Code:

int wth = 300;
int hgt = wth;

float centerX = -.5;
float centerY = 0;
float scale = 4;
float bailoutCircle = 10;

int maxIterations = 200;

color colorCore = color(0,0,0);

boolean bands = true;

color gradientColor1 = color(0,10,30);
color gradientColor2 = color(255,255,255);
float gradientShift = 1;
float gradientScale = 20;
float gradientNoise = .01;


void
setup() {
    size(wth, hgt);
    noLoop();
}


void
mousePressed() {
  centerX += (scale*float(mouseX)/wth-scale/2);
  centerY += (scale*float(mouseY)/hgt-scale/2);
  scale /= 2.0;
  println("X "+centerX+"   Y "+centerY+"   SC "+scale);
  redraw();
}

void
mouseDragged() {
  float dirX = mouseX - wth/2;
  float dirY = mouseY - hgt/2;
  float d = dist(0,0,dirX,dirY);
  dirX /= d;
  dirY /= d;
  centerX += dirX*scale*.01;
  centerY += dirY*scale*.01;
  scale /= 1.01;
  println("X "+centerX+"   Y "+centerY+"   SC "+scale);
  redraw();
}


void 
draw() {
    drawMandelbrot();
}
  
void 
drawMandelbrot() {
 
  loadPixels();
  
  float xmin = centerX - scale/2;
  float xmax = centerX + scale/2;
  float ymin = centerY - scale/2;
  float ymax = centerY + scale/2;
  
  // Calculate amount we increment x,y for each pixel
  float dx = (xmax - xmin)/wth;
  float dy = (ymax - ymin)/hgt;
  println("D "+dx);
  
  float aa = 0;
  float bb = 0;
  float y = ymin;
  for(int j = 0; j < height; j++) {
    float x = xmin;
    for(int i = 0;  i < width; i++) {
	
	// Now we test, as we iterate z = z^2 + cm does z tend towards infinity?
	float a = x;
	float b = y;
	int n = 0;
	while ( n < maxIterations ) {
	  float tmpa = a*a - b*b + x;
	  b = 2.0 * a * b + y;
	  a = tmpa;
	  aa = a*a;
	  bb = b*b;
	  if(aa + bb > bailoutCircle*bailoutCircle)  break;  // Bail
	  n++;
	}
	
	if (n == maxIterations) pixels[i+j*width] = colorCore;
	else {
	  if(bands) pixels[i+j*width] = color(abs(n*16 % 511 - 256)); 
	  else {
	    float gradient = n +  1.0f - log(log((aa+bb)))/log(2);
	    gradient = abs((gradientShift+gradient/gradientScale+random(gradientNoise))%2-1);
	    pixels[i+j*width] = gradientColor(gradientColor1,gradientColor2,gradient);
	  }
	}

	x += dx;
    }
    y += dy;
  }
  
  updatePixels();
  
}



color
gradientColor(color colorA, color colorB, float t) {
  t = constrain(t, 0, 1);
  float r = lerp(red(colorA),  red(colorB), t);
  float g = lerp(green(colorA), green(colorB), t);
  float b = lerp(blue(colorA),  blue(colorB), t);
  return color(r,g,b);
}

 

Isn't it possible to use double instead?

Btw: nice Mandelbrot

"The Java programming language provides two built-in classes for representing floating-point numbers: float, and double.

The "float" class takes 4 bytes of storage, and have 23 binary digits of precision. The "double" class takes 8 bytes of storage, and have 52 binary digits of precision."

I wasn't aware of the doubles - thanks!  

Reaching e-17 now - enough for some decent zoom fun on my current project.  Any advice or pointers about how to push this even further is still welcome though.

I also followed Ben's advice and did the bitshifting for the colors.


Code:

int wth = 1000;
int hgt = wth;
 
double centerX = -.5;
double centerY = 0;
double scale = 4;
double bailoutCircle = 10;


int maxIterations = 1000;
 
color colorCore = color(0,0,0);

// bands are a faster coloring method
boolean bands = false;

// if bands is false, a nice noised gradient will be used with following colors:
color gradientColor1 = color(0,10,30);
color gradientColor2 = color(255,255,255);
float gradientShift = 1;
float gradientScale = 20;
float gradientNoise = .01;
 
 
void
setup() {
    size(wth, hgt);
    noLoop();
}
 
 
void
mousePressed() {
  //centerX += (scale*double(mouseX)/wth-scale/2);
  centerX += (scale*(mouseX)/wth-scale/2);
  centerY += (scale*(mouseY)/hgt-scale/2);
  scale /= 2.0;
  println("X "+centerX+"   Y "+centerY+"   SC "+scale);
  redraw();
}
 
void
mouseDragged() {
  float dirX = mouseX - wth/2;
  float dirY = mouseY - hgt/2;
  float d = dist(0,0,dirX,dirY);
  dirX /= d;
  dirY /= d;
  centerX += dirX*scale*.01;
  centerY += dirY*scale*.01;
  scale /= 1.01;
  println("X "+centerX+"   Y "+centerY+"   SC "+scale);
  redraw();
}
 
 
void  
draw() {
    drawMandelbrot();
}
   
void  
drawMandelbrot() {
 
  loadPixels();
   
  double xmin = centerX - scale/2;
  double xmax = centerX + scale/2;
  double ymin = centerY - scale/2;
  double ymax = centerY + scale/2;
   
  // Calculate amount we increment x,y for each pixel
  double dx = (xmax - xmin)/wth;
  double dy = (ymax - ymin)/hgt;
  println("D "+dx);
   
  double aa = 0;
  double bb = 0;
  double y = ymin;
  for(int j = 0; j < height; j++) {
    double x = xmin;
    for(int i = 0;  i < width; i++) {
 
 // Now we test, as we iterate z = z^2 + cm does z tend towards infinity?
 double a = x;
 double b = y;
 int n = 0;
 while ( n < maxIterations ) {
   double tmpa = a*a - b*b + x;
   b = 2.0 * a * b + y;
   a = tmpa;
   aa = a*a;
   bb = b*b;
   if(aa + bb > bailoutCircle*bailoutCircle)  break;  // Bail
   n++;
 }
 
 if (n == maxIterations) pixels[i+j*width] = colorCore;
 else {
   if(bands) pixels[i+j*width] = color(abs(n*16 % 511 - 256));  
   else {
     float gradient = n +  1.0f - log(log((float)(aa+bb)))/log(2);
     gradient = abs((gradientShift+gradient/gradientScale+random(gradientNoise))%2-1);
     pixels[i+j*width] = gradientColor(gradientColor1,gradientColor2,gradient);
   }
 }
 
 x += dx;
    }
    y += dy;
  }
   
  updatePixels();
   
}
 
 
 
color
gradientColor(color colorA, color colorB, float t) {
  t = constrain(t, 0, 1);
  int r = int(lerp((colorA >> 16) & 0xFF,  (colorB >> 16) & 0xFF, t));
  int g = int(lerp((colorA >> 8 ) & 0xFF,  (colorB >> 8 ) & 0xFF, t));
  int b = int(lerp( colorA	  & 0xFF,   colorB	  & 0xFF, t));
//  println("T "+t+"   R "+r+"   G "+g+"   B "+b);
  int a = 255;
  return color((a<<24)|(r<<16)|(g<<8)|b);
}