Z-order curve and inverse Z-order for a cellular automaton

Question

I've written a program in Processing which simulates a cellular automaton in any number of dimensions (totalistic, nearest-neighbor, Moore neighborhood). To simplify things and make the program as general as possible, I store the states of the cells in a 1-dimensional integer array. To compute the n-dimensional neighbors of a cell, I treat the 1-dimensional array as a raster order (iterating the first coordinate from zero to the size of the space in question, then incrementing the next coordinate by one and iterating through the first one again, and so on.)

Since there are very few elegant ways to display n-dimensional data in the 2 dimensions of a computer screen (especially in this case, where n can in theory range from 2 to infinity and in practice can at least range from 2 to 5), I've decided to try using the Z-order, but everything I've read on the subject confuses me, and all the material seems to start with the assumption of Cartesian coordinates (x1,x2,x3...xn), which doesn't really apply in my case.

What I'd like is to be able to store the elements from the 1-dimensional cell array in another 1-dimensional array, sorted by n-dimensional Z-order, then convert that array to an array sorted by 2-dimensional Z-order, so that I can display it. I want to do this in an efficient manner (I'm trying to keep the program running at at least 0.5 frames per second, which is already tough in 5 dimensions, even for small spaces).

I'm aware of the bit-interleaving algorithm, but I'm not quite sure how to implement it in Processing without overflowing my 32-bit integers. Any help would be appreciated. If you want more details, I've included the relevant code below (and yes, I know it's messy; that's how I work):

void setup() {
  colorMode(HSB,255);
  background(0);
  frameRate(10);
  noStroke();
  ruledensity = random(1.1);
  celldensity = 0.1;
  dimension = 2;
  numstates = 5;
  scale = 20;
  rule = new int[int(pow(3.0,float(dimension)))*numstates];
  cells = new int[int(pow(float(scale),float(dimension)))];
  ncells = new int[int(pow(float(scale),float(dimension)))];
  zorder = new int[int(pow(float(scale),float(dimension)))];
  
  size(scale * 4,scale * 4);
  
  for (int a = 0; a < rule.length; a++) {
    if (random(1.0) < ruledensity) {
      rule[a] = 1 + floor(random(numstates - 1));
    } else {
      rule[a] = 0;
    }
  }
  
  for (int b = 0; b < cells.length; b++) {
    if (random(1.0) < celldensity) {
      cells[b] = 1 + floor(random(numstates - 1));
    } else {
      cells[b] = 0;
    }
    ncells[b] = cells[b];
  }
  
  
}

void draw() {
  background(0);
  for (int x = 0; x < cells.length; x++) {
    nsum = 0;
    for (int ex = 0; ex <= dimension; ex++) {
      for(int i = -1; i <= 1; i++) {
        nsum += cells[cyclize(x + i * int(pow(float(scale),float(ex-1))),cells.length)];
      }
    }
    ncells[x] = rule[nsum];
  }