ランダムに配置された粒子を格子で区切る

粒子をランダムに配置した後、空間を粒子が一つか0こしか含まないような格子に区切るプログラムを書きました。
クイックソートの要領で再帰を使って実装しました。

実はこのプログラムは別の数値計算に用いる予定で書いたのですが予想外に時間が掛かりそうなので単体で記録として残しておくことにしました。

#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
const int N = 1000;
const int x = 0, y = 1, z = 2;
typedef vector<double> Vec;
typedef vector<Vec> Mat;
typedef vector<Mat> Ten;
Mat X(N, Vec(3, 0));
Mat Dv;
void init() {
  for(int i = 0; i < N; i++) {
    for(int d = 0; d < 3; d++) {
      X[i][d] = drand48() * 20 - 10;
    }
  }
  return;
}
int how_many(double xmin, double xmax,
	     double ymin, double ymax,
	     double zmin, double zmax) {
  int count = 0;
  for(int i = 0; i < N; i++) {
    if((xmin < X[i][x] and X[i][x] < xmax) and
       (ymin < X[i][y] and X[i][y] < ymax) and
       (zmin < X[i][z] and X[i][z] < zmax)) {
      count ++;
    }
  }
  return count;
}
      
void dvide(double xmin, double xmax,
	   double ymin, double ymax,
	   double zmin, double zmax) {
  Vec tmp;
  tmp.push_back(xmin);
  tmp.push_back(xmax);
  tmp.push_back(ymin);
  tmp.push_back(ymax);
  tmp.push_back(zmin);
  tmp.push_back(zmax);
  Dv.push_back(tmp);
  if( 2 > how_many(xmin, xmax,
		   ymin, ymax,
		   zmin, zmax)) {
    return;
  } else {
    double dx = (xmax - xmin) / 2.0;
    double dy = (ymax - ymin) / 2.0;
    double dz = (zmax - zmin) / 2.0;
    dvide(xmin, xmin + dx,
	  ymin, ymin + dy,
	  zmin, zmin + dz);
    dvide(xmin, xmin + dx,
	  ymin, ymin + dy,
	  zmin + dx, xmax);
    dvide(xmin, xmin + dx,
	  ymin + dy, ymax,
	  zmin, zmin + dz);
    dvide(xmin, xmin + dx,
	  ymin + dy, ymax,
	  zmin + dz, zmax);
    dvide(xmin + dx, xmax,
	  ymin, ymin + dy,
	  zmin, zmin + dz);
    dvide(xmin + dx, xmax,
	  ymin, ymin + dy,
	  zmin + dz, zmax);
    dvide(xmin + dx, xmax,
	  ymin + dy, ymax,
	  zmin, zmin + dz);
    dvide(xmin + dx, xmax,
	  ymin + dy, ymax,
	  zmin + dz, zmax);
  }
}
void split() {
  Dv.clear();
  dvide(-10, 10,
	-10, 10,
	-10, 10);
}

int main(int argc, char **argv) {
  init();
  split();
  cout << Dv.size() << endl;
  return 0;
}