Cubic Newton-Raphson iteration, in Processing3, explorable in 4d, within the limits of JRE float precision.

ClassicFractals

// Parameter space for Cubic Newton-Raphson fractals --- review Stuff, it's probably helpful

/************ section for basic complex functions ****/

float[] cOne = {1.0,0.0};
float[] cZro = {0.0,0.0};
float[] cmOne = { -1.0, 0.0 };

float[] cmul ( float[] c1 , float[] c2 ) {
    float ans[] = new float[2];
    ans[0] = c1[0] * c2[0] - c1[1] * c2[1];
    ans[1] = c1[0] * c2[1] + c1[1] * c2[0];
    return ans; 
}

float[] rmul ( float r , float[] c1){
  float ans[] = new float[2];
  ans[0] = c1[0] * r ;
  ans[1] = c1[1] * r ;
  return ans;
}

float[] cadd ( float[] c1, float[] c2 ) {
  float ans[] = new float[2];
  ans[0] = c1[0] + c2[0];
  ans[1] = c1[1] + c2[1];
  return ans;
}

float[] cdif( float[] c1, float[] c2) {
  float ans[] = new float[2];
  ans[0] = c1[0] - c2[0];
  ans[1] = c1[1] - c2[1];
  return ans;
}

float[] conj ( float[] c1 ) {
  float[] ans = new float[2];
  ans[0] = c1[0];
  ans[1] = -c1[1];
  return ans;
}

float cnorm( float[] c1 ) {
  return (c1[0] * c1[0] + c1[1] * c1[1]);
}

float[] cinv ( float[] c1){
  return rmul ( 1/cnorm(c1), conj(c1) );
}

float[] cdiv ( float[] c1, float[] c2) {
  float ans[] = new float[2];
  float R = cnorm(c2);
  ans[0] = (c1[0]*c2[0] + c1[1]*c2[1])/R ;
  ans[1] = (-c1[0]*c2[1]+c1[1]*c2[0])/R;
  return ans;
}

float[] cexp ( float[] c1 ){
  float ans[] = new float[2];
  float R = exp(c1[0]);
  ans[0] = R * cos(c1[1]);
  ans[1] = R * sin(c1[1]);
  return ans;
}

float[] clog( float[] c1){
  float ans[] = new float[2];
  float R = sqrt(cnorm(c1));
  ans[0] = log(R);
  ans[1] = 2 * atan(c1[1]/(R + c1[0]));
  return ans;
}
/************ section for coloring a pixel ***************/

/*** settling on
 the polynomials P(z) = (z-1)(z+1)(z-c) = (z^2-1)(z-c) = z^3-cz^2-z+c 
 satisfying P'(z) = 3 z^2 - 2 c z - 1;
 and because we know what the roots are, testing literal closeness to a root
 is "easy"
***/

float[] nrIter( float[] c , float[] z) {
  float[] numer = cmul(cdif(z,cOne),cmul(cadd(z,cOne),cdif(z,c)));
  float[] denom = cadd( cmul(cadd(z,cOne),cdif(z,c)) , cadd(cmul(cdif(z,cOne),cdif(z,c)) , cmul(cdif(z,cOne),cadd(z,cOne))));
  return cdiv ( numer , denom );
}

int nrColor( float[] c , float[] z, int maxIter ) {
  int countIter = 0, ans;
  float[] move;
  float p , q , r, s;
  for (countIter = 0 ; countIter < maxIter ; countIter++ ){
    move = nrIter(c,z);
    z = cdif( z , move );
    p = cnorm(cdif(z,c));
    q = cnorm(cdif(z,cOne));
    r = cnorm(cadd(z,cOne));
    ans = 0;
    if ( q < p ) { s = p ; p = q ; q = s; ans = 1; } // now p < q
    if ( r < p ) { s = p ; p = r ; r = s; ans = 2; } // now p < min(r,q);
    if ( p < min ( min(q,r)/5, 0.5 ) ) return ans;
  }
  return 5;
}

int nrparamColor( float[] c , int maxIter ){
  return nrColor( c , rmul(1/3.0, c), maxIter);
  /* ... we used to have the above iteration here, instead.
 */
}

int jparamColor( float[] c, int maxIter){

/* just as a sanity-check; made several typos in the "nrXYZ",
wanted to make sure the problems weren't in the Complex-variables mess above */

  int countIter = 0;
  float[] z = cZro;
  for ( countIter = 0; countIter < maxIter; countIter++ ){
    z = cadd( cmul(z,z), c) ;
    if (cnorm(z) > 5) break;
  }
  return ( maxIter - countIter );
}

float theC[] = { 0 , sqrt(3) }; // had to initialize to Something; in the end, this doesn't get used anymore.
int theLastClickX = 200;
int theLastClickY = 200 - 173;
int doneState = 0;

/* the draw() function can (in theory) do three things: 
draw a classic Mandelbrot; draw a cubic NR picture; or,
draw NR cubic parameter space. (Grant Sanderson did a 
nice Video Two-Parter describing that... which is sort-of
why I wrote the present Thing.  That, and I can't find a 
good general fractals explorer packaged for archlinux... but nevermind.

middle-clicking switches between particular-NR and paramspace NR;
and in particular, the Third Root in particular-NR is the center of the current
paramspace NR, which is initially 0.0+0.0i;
left/right clicking zooms in
*/

int which = 2;

float moveX[] = { 0 , 0 , 0 };
float moveY[] = { 0 , 0 , 0 };
float zoom[] = { 1 , 1 , 1 };


void setup(){
  size(600,600);
  stroke(200,200,200);
  theLastClickX = width/2;
  theLastClickY = height/2 - floor(sqrt(3) * 100);
  noFill();
}

void draw(){
  int i, j, n;
  int iters = 320;
  PGraphics bgFrame = createGraphics(width,height);
  bgFrame.beginDraw();
  bgFrame.loadPixels();
  float[] myC = new float[2];
  if ( which == 2 ) { theC[0] = moveX[1] ; theC[1] = moveY[1] ; } 
  for ( i = 0 ; i < width ; i++ ){
      myC[0] = moveX[which] + ( i - width/2.0 ) / (zoom[which] * 100.0) ;
    for( j = 0 ; j < height ; j++) {
      myC[1] = moveY[which] + ( height/2.0 - j ) / (zoom[which] * 100.0);
      switch (which) {
        case 0 : 
            n = jparamColor(cdiv(cOne,myC),200);
            bgFrame.pixels[j * width + i] = color( n * 255/400 , n*n*255/16000, 128 * ( 1 + sin(5.0 * n/iters )));
            break;
        case 1 : // 
          switch ( n = nrparamColor(myC,iters)) {
            case 0 : bgFrame.pixels[j*width + i] = color(100, 10, 30); break;
            case 1 : bgFrame.pixels[j*width + i] = color(10, 30, 100); break;
            case 2 : bgFrame.pixels[j*width + i] = color(30, 100, 10); break;
            default : bgFrame.pixels[j*width + i] = color(0,0,0); break;
          } break;
        case 2: 
          switch ( n = nrColor(theC,myC,iters)) {
            case 0 : bgFrame.pixels[j*width + i] = color(100, 10, 30); break;
            case 1 : bgFrame.pixels[j*width + i] = color(10, 30, 100); break;
            case 2 : bgFrame.pixels[j*width + i] = color(30, 100, 10); break;
            default : bgFrame.pixels[j*width + i] = color( 0 ,0 ,0 ) ; break;
          } break;
      }
    }
  }
  bgFrame.updatePixels();
  image(bgFrame,0,0);
  bgFrame.endDraw();
}

void mouseReleased(){
  switch (mouseButton){
    case LEFT: // zoom in
      zoom[which] *= 2 ;
      moveX[which] = moveX[which] + ( mouseX - width/2.0 ) / (zoom[which] * 100.0) ; // these leave the clicked pixel fixed.
      moveY[which] = moveY[which] + ( height/2.0 - mouseY ) / (zoom[which] * 100.0); 
      break;
    case RIGHT: // zoom out
      moveX[which] = moveX[which] - ( mouseX - width/2.0 ) / (zoom[which] * 100.0) ; // and so do these
      moveY[which] = moveY[which] - ( height/2.0 - mouseY ) / (zoom[which] * 100.0);
      zoom[which] /= 2 ;
      break;
    case CENTER: // switch mode, somehow. 
      which = 3 - which;
      break;
  }
}