My code

… was written in C++11 and can be compiled with G++, Clang++, Visual C++. You can download it, too. Or just jump to my GitHub repository.

       
#include <iostream>
       
#include <vector>
       
#include <set>
       
#include <map>
       
 
       
// a cuboid with 1+ added layers

       
struct Cuboid
       
{
       
  // cuboid size
       
  const unsigned short width;
       
  const unsigned short height;
       
  const unsigned short depth;
       
  // wrapped several times
       
  const unsigned short layers;
       
  // compute volume with calculated()
       
  const unsigned int   hullVolume;
       
 
       
  // create new cuboid
       
  Cuboid(unsigned short width_, unsigned short height_, unsigned short depth_, unsigned short layers_)
       
  : width(width_), height(height_), depth(depth_), layers(layers_),
       
    hullVolume(calculate())
       
  {}
       
 
       
  // calculate volume of the most recent / outermost layer
       
  unsigned int calculate() const
       
  {
       
    // innermost hull
       
    auto result = 2 * ((unsigned int)width  * height +
       
                       (unsigned int)height * depth  +
       
                       (unsigned int)depth  * width);
       
 
       
    // need just one layer ?
       
    if (layers == 1)
       
      return result;
       
 
       
    // add more layers
       
    auto increment = 4 * (width + height + depth);
       
    for (unsigned short current = 2; current <= layers; current++)
       
    {
       
      result    += increment;
       
      increment += 8;
       
    }
       
 
       
    return result;
       
  }
       
 
       
  // for std::set
       
  bool operator<(const Cuboid& other) const
       
  {
       
    // volume
       
    if (hullVolume < other.hullVolume) return true;
       
    if (hullVolume > other.hullVolume) return false;
       
    // the following lines are arbitrary => they only need to be consistent
       
    // sides
       
    if (width  < other.width)  return true;
       
    if (width  > other.width)  return false;
       
    if (height < other.height) return true;
       
    if (height > other.height) return false;
       
    if (depth  < other.depth)  return true;
       
    if (depth  > other.depth)  return false;
       
    // layers
       
    return layers < other.layers;
       
  }
       
};
       
 
       
// list of cuboids, sorted by volume of their outermost layer

       
std::set<Cuboid> todo;
       
// count how many cuboid layers have a certain volume (looking for 1000)

       
std::map<unsigned int, unsigned int> count; // could be replaced by std::vector if you can estimate the upper limit
       
// note: this information is implicitely stored in "todo" as well: count the consecutive cuboids with the same hullVolume

       
 
       
// add a cuboid with a certain amount of layers, return number of cuboids with the same hull volume

       
bool add(unsigned short width, unsigned short height, unsigned short depth, unsigned short layers)
       
{
       
  // accept only width >= height >= depth
       
  if (width  < height)
       
    return false;
       
  if (height < depth)
       
    return false;
       
 
       
  // create new cuboid
       
  Cuboid cuboid(width, height, depth, layers);
       
 
       
  // already queued ?
       
  if (todo.count(cuboid) != 0)
       
    return false;
       
 
       
  // add cuboid
       
  todo.insert(cuboid);
       
  count[cuboid.hullVolume]++;
       
 
       
  // successfully added
       
  return true;
       
}
       
 
       
const unsigned int Failed = 0;
       
// find the correct result if you know an upper bound maxVolume, if not found return Failed

       
unsigned int fast(unsigned int stopIf, unsigned int maxVolume)
       
{
       
  // allocate enough memory
       
  std::vector<unsigned int> count(maxVolume + 1, 0);
       
 
       
  // just four nested loops to generate all cuboids with hulls <= maxVolume
       
  // ensure that x >= y >= z to avoid enumerating some volume multiple times
       
  for (unsigned short x = 1; Cuboid(x,1,1,1).hullVolume <= maxVolume; x++)
       
    for (unsigned short y = 1; y <= x && Cuboid(x,y,1,1).hullVolume <= maxVolume; y++)
       
      for (unsigned short z = 1; z <= y && Cuboid(x,y,z,1).hullVolume <= maxVolume; z++)
       
        for (unsigned short layer = 1; ; layer++)
       
        {
       
          auto current = Cuboid(x,y,z,layer).hullVolume;
       
          if (current > maxVolume)
       
            break;
       
 
       
          count[current]++;
       
        }
       
 
       
  // find smallest hull volume that appears exactly 1000 times
       
  for (size_t i = 0; i < count.size(); i++)
       
    if (count[i] == stopIf)
       
      return i;
       
 
       
  // failed
       
  return Failed;
       
}
       
 
       
// forward declaration of my old helper code
       
//unsigned int naive(unsigned int x, unsigned int y, unsigned int z, unsigned int layers);

       
 
       
int main()
       
{
       
  // start old code to analyze a single cuboid
       
  //naive(3,2,1, 5); return 0;
       
 
       
  // look for that number of cuboids with the same volume
       
  unsigned int stopIf = 1000;
       
  std::cin >> stopIf;
       
 
       
  // ---------- fast enumeration ----------
       
 
       
  // find solution if an upper limit of the result is known
       
#define FAST_ENUMERATION
       
#ifdef  FAST_ENUMERATION
       
  const unsigned int StepSize = 10000; // I deliberately choose not the optimal value !
       
  // assume that the result is below 10000, then 20000, ... until fast() returns > 0
       
  for (unsigned int maxVolume = StepSize; ; maxVolume += StepSize)
       
  {
       
    auto current = fast(stopIf, maxVolume); // returns 0 if not found
       
    // found ?
       
    if (current != Failed)
       
    {
       
      std::cout << current << std::endl;
       
      return 0;
       
    }
       
  }
       
#endif
       
 
       
  // ---------- slower enumeration ----------
       
 
       
  // note: I found the correct result with this approach but later created the faster enumeration you find above
       
 
       
  // initial cuboid: one unit wide, high and deep with one layer around it
       
  add(1,1,1,1);
       
  unsigned int volume = 0;
       
  while (count[volume] != stopIf)
       
  {
       
    // pick cuboid with smallest volume and remove it
       
    auto current = *(todo.begin());
       
    todo.erase(todo.begin());
       
 
       
    volume = current.hullVolume;
       
 
       
    // enlarge cuboid
       
    add(current.width + 1, current.height, current.depth, current.layers);
       
    add(current.width, current.height + 1, current.depth, current.layers);
       
    add(current.width, current.height, current.depth + 1, current.layers);
       
    // add one layer
       
    add(current.width, current.height, current.depth, current.layers + 1);
       
  }
       
 
       
  // display result
       
  std::cout << volume << std::endl;
       
  return 0;
       
}
       
 
       
// ---------- compute hull volume, slow approach ----------

       
 
       
// a single cube that is part of an cuboid (more or less a coordinate)

       
struct Cube
       
{
       
  // coordinates
       
  const int x, y, z;
       
 
       
  // create new cube
       
  Cube(int xx, int yy, int zz)
       
  : x(xx), y(yy), z(zz)
       
  {}
       
 
       
  // for std::set
       
  bool operator<(const Cube& other) const
       
  {
       
    if (x < other.x) return true;
       
    if (x > other.x) return false;
       
    if (y < other.y) return true;
       
    if (y > other.y) return false;
       
    return z < other.z;
       
  }
       
};
       
 
       
unsigned int naive(unsigned int x, unsigned int y, unsigned int z, unsigned int layers)
       
{
       
  // keep track of all cubes belonging to the shape
       
  std::set<Cube> cuboid;
       
  // cubes added in the most recent layer
       
  std::vector<Cube> hull;
       
 
       
  // initial shape
       
  for (int xx = 0; xx < (int)x; xx++)
       
    for (int yy = 0; yy < (int)y; yy++)
       
      for (int zz = 0; zz < (int)z; zz++)
       
      {
       
        Cube current(xx, yy, zz);
       
        cuboid.insert(current);
       
        hull.push_back(current);
       
      }
       
 
       
  // add desired numbers of layers
       
  for (unsigned int layer = 1; layer <= layers; layer++)
       
  {
       
    std::vector<Cube> next;
       
    for (auto current : hull)
       
    {
       
      // grow in each direction x,y,z
       
      Cube left (current.x - 1, current.y, current.z);
       
      Cube right(current.x + 1, current.y, current.z);
       
      Cube down (current.x, current.y - 1, current.z);
       
      Cube up   (current.x, current.y + 1, current.z);
       
      Cube back (current.x, current.y, current.z - 1);
       
      Cube front(current.x, current.y, current.z + 1);
       
 
       
      // add cube (if not already existing)
       
      if (cuboid.count(left)  == 0)
       
      {
       
        cuboid.insert(left);
       
        next.push_back(left);
       
      }
       
      if (cuboid.count(right) == 0)
       
      {
       
        cuboid.insert(right);
       
        next.push_back(right);
       
      }
       
      if (cuboid.count(down)  == 0)
       
      {
       
        cuboid.insert(down);
       
        next.push_back(down);
       
      }
       
      if (cuboid.count(up)    == 0)
       
      {
       
        cuboid.insert(up);
       
        next.push_back(up);
       
      }
       
      if (cuboid.count(back)  == 0)
       
      {
       
        cuboid.insert(back);
       
        next.push_back(back);
       
      }
       
      if (cuboid.count(front) == 0)
       
      {
       
        cuboid.insert(front);
       
        next.push_back(front);
       
      }
       
    }
       
 
       
    std::cout << "layer " << layer << ": "
       
              << "C(" << x << "," << y << "," << z << ")=" << next.size();
       
    if (layer > 1)
       
      std::cout << " delta=" << next.size() - hull.size();
       
    std::cout << std::endl;
       
 
       
    // next iteration
       
    hull = std::move(next);
       
  }
       
 
       
  return hull.size();
       
}

This solution contains 40 empty lines, 57 comments and 7 preprocessor commands.