text:
Extending OEIS: Counting Diamond Tilings
I promise, this will be my last challenge about diamong tilings (for a while, anyway). On the bright side, this challenge doesn't have anything to do with ASCII art, and is not a code golf either, so this is actually completely different.
So just as a reminder, every hexagon can be titled with three different diamonds:

An interesting question to ask is how many of these tilings exist for a given hexagon size. It seems these numbers have been studied fairly thoroughly and can be found in OEIS A008793.
However, the problem gets trickier if we ask how many tilings exist up to rotation and reflection. For instance, for side length N = 2, the following 20 tilings exist:
   ____     ____     ____     ____     ____     ____     ____     ____     ____     ____
  /\_\_\   /\_\_\   /\_\_\   /\_\_\   /_/\_\   /_/\_\   /\_\_\   /_/\_\   /_/\_\   /_/\_\
 /\/\_\_\ /\/_/\_\ /\/_/_/\ /\/_/\_\ /\_\/\_\ /\_\/_/\ /\/_/_/\ /\_\/\_\ /\_\/_/\ /_/\/\_\
 \/\/_/_/ \/\_\/_/ \/\_\_\/ \/_/\/_/ \/\_\/_/ \/\_\_\/ \/_/\_\/ \/_/\/_/ \/_/\_\/ \_\/\/_/
  \/_/_/   \/_/_/   \/_/_/   \_\/_/   \/_/_/   \/_/_/   \_\/_/   \_\/_/   \_\/_/   \_\/_/
   ____     ____     ____     ____     ____     ____     ____     ____     ____     ____
  /_/_/\   /\_\_\   /_/\_\   /_/_/\   /_/\_\   /_/\_\   /_/_/\   /_/_/\   /_/_/\   /_/_/\
 /\_\_\/\ /\/_/_/\ /_/\/_/\ /\_\_\/\ /\_\/_/\ /_/\/_/\ /_/\_\/\ /\_\_\/\ /_/\_\/\ /_/_/\/\
 \/\_\_\/ \/_/_/\/ \_\/\_\/ \/_/\_\/ \/_/_/\/ \_\/_/\/ \_\/\_\/ \/_/_/\/ \_\/_/\/ \_\_\/\/
  \/_/_/   \_\_\/   \_\/_/   \_\/_/   \_\_\/   \_\_\/   \_\/_/   \_\_\/   \_\_\/   \_\_\/

But many of these are identical under rotation and reflection. If we take these symmetries into account, only 6 distinct tilings remain:
   ____     ____     ____     ____     ____     ____
  /\_\_\   /\_\_\   /\_\_\   /_/\_\   /_/\_\   /_/\_\
 /\/\_\_\ /\/_/\_\ /\/_/_/\ /\_\/_/\ /\_\/_/\ /_/\/\_\
 \/\/_/_/ \/\_\/_/ \/\_\_\/ \/\_\_\/ \/_/\_\/ \_\/\/_/
  \/_/_/   \/_/_/   \/_/_/   \/_/_/   \_\/_/   \_\/_/

   2        2        6        6        1        3

where the numbers indicate the multiplicity of each tiling. Note that for larger hexagons there are also tilings with multiplicities 4 and 12.
It appears that the number of tilings up to symmetry has been studied less thoroughly. The OEIS entry A066931 only lists the five terms:
1, 1, 6, 113, 20174

where the first term is for side length N = 0 and the last term for side length N = 4.
I'm sure we can do better than that!
Your task is to compute the number of tilings for a given side length.
This is fastest-code. Your score will be the highest side-length N for which your code produces the correct result within 30 minutes on my machine. In case of a tie, I will accept the submission which produces the result for that N fastest.
As usual, you must not hardcode results you already know to win the tie-breaker. The algorithm that solve N = 3 should be identical to the one that solves N = 5.
Your submission must not use more than 4GB of memory. I will give some leeway on this if you're operating close to that limit, but if you're consistently above that limit, or if you spike significantly beyond it, I will not count that N for your submission.
I will test all submissions on my Windows 8 machine, so make sure your language of choice is freely available on Windows. The only exception to this is Mathematica (because I happen to have a licence for it). Please include instructions for how to compile/run your code.
Of course, feel free to compute more terms in your own time (for science, and for others to check their numbers against), but your answer's score will be determined in those 30 minutes.

Note that since N = 6 gives an output of more than 10^12, a non-constructive solution is almost certainly necessary to get that far.

@PeterTaylor I was hoping that would allow more room for improvement. Maybe a couple of simple constructive answers first that can do N = 5 to gain more insight into the problem, and then potentially hybrid approaches that don't need to construct all tilings but can extrapolate the total number from a few constructed ones... and then maybe something analytic if we're really lucky. :)

At the risk of stating the obvious, it seems to me that each such tiling corresponds to a projection of an assemblage of unit cubes as viewed from a distant point of view, e.g, from (100,-100,100).  I find that this lightens the burden of constructing tilings.

@DavidCarraher Indeed. More specifically, such an arrangement of unit cubes is a 3D Young diagram. (Maybe that helps someone.)

@DavidCarraher If you look hard enough at the big hexagon, you'll see there's 2 different ways to interpret it as a Young diagram. The obvious way (for me at least) is to see a flat area at the top and left with a 2x2x1 cuboid missing from the top left corner. But there's another way of seeing it: an empty zone in that area, with a 2x2x1 cuboid sitting in it. Tilting 60 degrees may help. It hurts my eyes but I think the two young diagrams fit together, possibly by reflection of one of them. OEIS A008793 is very careful with its wording: "number of plane partitions whose young diagrams..."

Yes, there is a figure-ground phenomenon with several of the figures.  This is one way to (manually) find some isomorphic figures.  I realize that the task is planar. But the 3D interpretation helps simplify the identification of valid 2D tilings.

To throw out a possible approach: since there's a formula for the total number of tilings, to compute the number of classes up to rotation and reflection, it suffices to count the tilings that satisfy each type of nontrivial symmetry. I expect only a small fraction of tilings to have a symmetry.

I'll be impressed with any solution submitted, but I'm sure some of the people on Project Euler would be able to produce good results.

I'm sorry to say that your statement regarding A066931 having 5 terms online is out of date ;) (Led here by your corresponding post in this post, so no surprises)

Algebra, graph theory, Möbius inversion, research, and Java
The symmetry group of the hexagon is the dihedral group of order 12, and is generated by a 60 degree rotation and a mirror flip across a diameter. It has 16 subgroups, but some of them are in non-trivial conjugacy groups (the ones which only have reflections have 3 choices of axis), so there are 10 fundamentally different symmetries which a tiling of the hexagon can have:

The number of diamond tilings of a subset of a triangular lattice can be calculated as a determinant, so my initial approach was to set up one determinant for each of the symmetries of the hexagon, to calculate the number of tilings which have at least those symmetries; and then use Möbius inversion in the incidence algebra of their poset (basically a generalisation of the inclusion-exclusion principle) to work out the number of tilings whose symmetry group is exactly each of the 10 cases. However, some of the symmetries have nasty edge conditions, so I was forced to sum over exponentially many determinants. Fortunately, the values obtained for n < 10 gave me enough data to be able to identify relevant sequences in OEIS and piece together a closed form (for some value of "closed" which allows finite products). There's a bit of discussion of the sequences, and references to proofs, in the formal writeup which I prepared to justify the OEIS sequence updates.
Once the double-counting is taken care of, it turns out that four of the ten values cancel out neatly, so we only have to calculate the remaining six and then do a weighted sum.
This code takes under 30 seconds for N=1000 on my machine.
import java.math.BigInteger;

public class OptimisedCounter {
    private static int[] minp = new int[2];

    public static void main(String[] args) {
        if (args.length > 0) {
            for (String arg : args) System.out.println(count(Integer.parseInt(arg)));
        }
        else {
            for (int n = 0; n < 16; n++) {
                System.out.format("%d\t%s\n", n, count(n));
            }
        }
    }

    private static BigInteger count(int n) {
        if (n == 0) return BigInteger.ONE;

        if (minp.length < 3*n) {
            int[] wider = new int[3*n];
            System.arraycopy(minp, 0, wider, 0, minp.length);
            for (int x = minp.length; x < wider.length; x++) {
                // Find the smallest prime which divides x
                for (wider[x] = 2; x % wider[x] != 0; wider[x]++) { /* Do nothing */ }
            }
            minp = wider;
        }

        BigInteger E = countE(n), R2 = countR2(n), F = countF(n), R3 = countR3(n), R = countR(n), FR = countFR(n);
        BigInteger sum = E.add(R3);
        sum = sum.add(R2.add(R).multiply(BigInteger.valueOf(2)));
        sum = sum.add(F.add(FR).multiply(BigInteger.valueOf(3)));
        return sum.divide(BigInteger.valueOf(12));
    }

    private static BigInteger countE(int n) {
        int[] w = new int[3*n];
        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j <= i + n; j++) w[j]--;
            for (int j = i + n + 1; j <= i + 2*n; j++) w[j]++;
        }
        return powerProd(w);
    }

    private static BigInteger countR2(int n) {
        int[] w = new int[3*n];
        for (int i = 0; i < n; i++) {
            w[3*i+2]++;
            for (int j = 3*i + 1; j <= 2*i + n + 1; j++) w[j]--;
            for (int j = 2*i + n + 1; j <= i + n + n; j++) w[j]++;
        }
        return powerProd(w);
    }

    private static BigInteger countF(int n) {
        int[] w = new int[3*n];
        for (int i = 0; i < n; i++) {
            for (int j = 2*i + 1; j <= 2*i + n; j++) w[j]--;
            for (int j = i + n + 1; j <= i + 2*n; j++) w[j]++;
        }
        return powerProd(w);
    }

    private static BigInteger countR3(int n) {
        if ((n & 1) == 1) return BigInteger.ZERO;
        return countE(n / 2).pow(2);
    }

    private static BigInteger countR(int n) {
        if ((n & 1) == 1) return BigInteger.ZERO;
        int m = n / 2;
        int[] w = new int[3*m-1];
        for (int i = 0; i < m; i++) {
            for (int j = 1; j <= 3*i+1; j++) w[j] += 2;
            for (int j = 1; j <= i + m; j++) w[j] -= 2;
        }
        return powerProd(w);
    }

    private static BigInteger countFR(int n) {
        if ((n & 1) == 1) return BigInteger.ZERO;
        int m = n / 2;
        int[] w = new int[3*n-2];
        for (int j = 1; j <= m; j++) w[j]--;
        for (int j = 2*m; j <= 3*m-1; j++) w[j]++;
        for (int i = 0; i <= 2*m-3; i++) {
            for (int j = i + 2*m + 1; j <= i + 4*m; j++) w[j]++;
            for (int j = 2*i + 3; j <= 2*i + 2*m + 2; j++) w[j]--;
        }
        return powerProd(w);
    }

    private static BigInteger powerProd(int[] w) {
        BigInteger result = BigInteger.ONE;
        for (int x = w.length - 1; x > 1; x--) {
            if (w[x] == 0) continue;

            int p = minp[x];
            if (p == x) result = result.multiply(BigInteger.valueOf(p).pow(w[p]));
            else {
                // Redistribute it. This should ensure we avoid negatives.
                w[p] += w[x];
                w[x / p] += w[x];
            }
        }

        return result;
    }
}

You are truly a god amongst mortals. I hope to see your solution published in a prestigious journal.

This is awesome. BTW my (currently not posted) code gives 22306956 for N=5: 22231176(12)+275(4)+75328(6)+352(2), a discrepancy of 1, which is odd. I've no idea you're doing here, is it suitable for a a breakdown by symmetries? For N=4 I'm 16 lower than you and https://oeis.org/A066931/a066931.txt From that reference it seems I have 16 too many of multiplicity 12, which I need to convert to 32 of multiplicity 6. I'm not too surprised, even N are more difficult for me. But I have no problems with odd N and I get right answers for 0<N<4. Will look for obvious problems and post my code tomorrow.

@steveverrill, if I understand the notation, for N=5 I make it 22231176(12) + 75328(6) + 275(4) + 176(2). I think you're failing to divide the index 2 ones by 2. (FWIW for odd numbers they all have an axis of symmetry passing through two vertices and rotational symmetry of order 3).

@steveverrill, and for N=4 your discrepancy seems to be a perfect fit for the number which have an axis of symmetry passing through the midpoints of two edges.

Impressive that you solved this.  I'm hoping you will eventually post an answer that non-mathematicians can follow.

Thanks, N=5 is solved. As I generate all tilings, I was using floating point to reduce to the unique tilings. I now check my divisions in a different way, answer revised.

Problem with N=4 seems to be what I expected. I knew some solids were self-complement and I assumed they were all antisymmetric with respect to (cube inversion)x(identity) which gives me (i think) 57 tilings with just a 2-fold rotation axis. I guess there are solids that are antisymmetric with respect to (cube inversion)x(some other symetry) these would be the ones I missed. I can't envisage the tilings with only a line of symmetry through the edges, is there an example anywhere? The title of OEIS A008793 seems to be clear that there's a 1:1 mapping between plane partitions and hexagon tilings

@steveverrill, 8 columns seems too few. There are 10 fundamentally different symmetries (including "no symmetry"), not 8. See e.g. http://dedekind.mit.edu/~rstan/pubs/pubfiles/65.pdf . For a tiling with only a line of symmetry through the edges, I think http://pastebin.com/raw.php?i=Y6F8bEEP is correct.

The paper looks like a tough read but by my own reckoning I'm missing 2: solids antisymmetric under (cube inversion)x(1/3 turn) only, and solids that are antisymmetric under (cube inversion)x(reflection through one of the C3v planes) only. Will see if my thoughts match the paper later. Thanks for the pastebin. Although it is a symmetric tiling, my eyes insist on interpreting it as a Young diagram that 1) contains exactly (n**3)/2 cells 2)is asymmetric and 3)is antisymmetric under some operation I can't identify (but definitely not inversion.) This makes it hard to conceive these tilings myself

@steveverrill, I've approached the problem in terms of tilings rather than plane partitions, and that might be conceptually simpler. It certainly made it fairly simple for me to reason about the symmetries.

C
Introduction
As commented by David Carraher, the simplest way of analysing the hexagon tiling seemed to be to take advantage of its isomorphism with the 3 dimensional Young Diagram, essentially an x,y square filled with integer height bars whose z heights must stay the same or increase as the z axis is approached.
I started with an algorithm for finding the totals which is more amenable to adaptation for symmetry counting than the published algorithm, which is based on a bias to one of the three cartesian axes.
Algorithm
I start by filling the cells of the x,y, and z planes with 1's, while the rest of the area contains zeroes. Once that is done, I build up the pattern layer by layer, with each layer containing the cells which have a common 3D manhattan distance from the origin. A cell can only contain a 1 if the three cells below it also contain a 1. if any of them contains a 0, then the cell must be a 0.
The advantage of building the pattern up in this way is that each layer is symmetrical about the x=y=z line. This means that each layer can be checked independently for symmetry.
Symmetry checking
The symmetries of the solid are as follows: 3 fold rotation about the x=y=z line --> 3 fold rotation about the hexagon centre; and 3 x reflections about the 3 planes containing the x=y=z line and each of the axes x,y,z --> reflection about the lines through the hexagon corners.
This adds up only to 6 fold symmetry. To get the full symmetry of the hexagon, another kind of symmetry must be considered. Each solid (built from 1's) has a complementary solid (built from 0's). Where N is odd, the complementary solid must be different from the original solid (because it is not possible for them to have the same number of cubes). Yet when the complementary solid is turned round, it will be found that its 2D representation as a diamond tiling is identical (except for a 2-fold symmetry operation) to the original solid. Where N is even, it is possible for the solid to be self-inverse.
This can be seen in the examples for N=2 in the question. If viewed from the left, the first hexagon looks like a solid cube with 8 small cubes, while the last hexagon looks like an empty shell with 0 small cubes. If viewed from the right, the reverse is true. The 3rd, 4th and 5th hexagons and the 16th, 17th and 18th hexagons look like they contain either 2 or 6 cubes, and thus they complement each other in 3 dimensions. They are related to each other in 2 dimensions by a 2-fold symmetry operation (2-fold rotation, or reflection about an axis through the hexagon's edges.) On the other hand the 9th,10th,11th and 12th hexagons show 3D patterns which are their own complements, and therefore have a higher symmetry (these are therefore the only patterns with odd multiplicity).
Note that having (N^3)/2 cubes is a necessary condition to be self complement, but in general it is not a sufficient condition if N>2. The result of all this is that for odd N, tilings always occur in pairs (N^3)/2 cubes must be carefully inspected.
Current code (generates the right total for N=1,2,3,5. Error as discussed for N=4.)
int n;                     //side length

char t[11][11][11];        //grid sized for N up to 10

int q[29][192], r[29];     //tables of coordinates for up to 10*3-2=28 layers

int c[9];                  //counts arrangements found by symmetry class. c[8] contains total.

//recursive layer counting function. m= manhattan distance, e= number of cells in previous layers, s=symmetry class.
void f(int m,int e,int s){

  int u[64], v[64], w[64]; //shortlists for x,y,z coordinates of cells in this layer
  int j=0;
  int x,y,z;

  for (int i=r[m]*3; i; i-=3){
    // get a set of coordinates for a cell in the current layer.
    x=q[m][i-3]; y= q[m][i-2]; z= q[m][i-1];
    // if the three cells in the previous layer are filled, add it to the shortlist u[],v[],w[]. j indicates the length of the shortlist.
    if (t[x][y][z-1] && t[x][y-1][z] && t[x-1][y][z]) u[j]=x, v[j]=y, w[j++]=z ;
  }

  // there are 1<<j possible arrangements for this layer.
  for (int i = 1 << j; i--;) {

    int d = 0;

    // for each value of i, set the 1's bits of t[] to the 1's bits of i. Count the number of 1's into d as we go.
    for (int k = j; k--;) d+=(t[u[k]][v[k]][w[k]]=(i>>k)&1);

    // we have no interest in i=0 as it is the empty layer and therefore the same as the previous recursion step.
    // Still we loop through it to ensure t[] is properly cleared.

    if(i>0){
      int s1=s;    //local copy of symmetry class. 1's bit for 3 fold rotation, 2's bit for reflection in y axis.
      int sc=0;    //symmetry of self-complement.

      //if previous layers were symmetrical, test if the symmetry has been reduced by the current layer
      if (s1) for (int k = j; k--;) s1 &= (t[u[k]][v[k]][w[k]]==t[w[k]][u[k]][v[k]]) | (t[u[k]][v[k]][w[k]]==t[w[k]][v[k]][u[k]])<<1;

      //if exactly half the cells are filled, test for self complement
      if ((e+d)*2==n*n*n){
        sc=1;
        for(int A=1; A<=(n>>1); A++)for(int B=1; B<=n; B++)for(int C=1; C<=n; C++) sc&=t[A][B][C]^t[n+1-A][n+1-B][n+1-C];
      }

      //increment counters for total and for symmetry class.
      c[8]++; c[s1+(sc<<2)]++;

      //uncomment for graphic display of each block stacking with metadata. not recommended for n>3.
      //printf("m=%d  j=%d  i=%d c1=%d-2*%d=%d c3=%d cy=%d(cs=%d) c3v=%d ctot=%d\n",m,j,i,c[0],c[2],c[0]-2*c[2],c[1],c[2],c[2]*3,c[3],c[8]);
      //printf("m=%d  j=%d  i=%d C1=%d-2*%d=%d C3=%d CY=%d(CS=%d) C3V=%d ctot=%d\n",m,j,i,c[4],c[6],c[4]-2*c[6],c[5],c[6],c[6]*3,c[7],c[8]);
      //for (int A = 0; A<4; A++, puts(""))for (int B = 0; B<4; B++, printf(" "))for (int C = 0; C<4; C++) printf("%c",34+t[A][B][C]);

      //recurse to next level.
      if(m<n*3-2)f(m + 1,e+d,s1);

    }
  }
}

main()
{
  scanf("%d",&n);

  int x,y,z;

  // Fill x,y and z planes of t[] with 1's
  for (int a=0; a<9; a++) for (int b=0; b<9; b++) t[a][b][0]= t[0][a][b]= t[b][0][a]= 1;

  // Build table of coordinates for each manhattan layer
  for (int m=1; m < n*3-1; m++){
    printf("m=%d : ",m);
    int j=0;
    for (x = 1; x <= n; x++) for (y = 1; y <= n; y++) {
      z=m+2-x-y;
      if (z>0 && z <= n) q[m][j++] = x, q[m][j++] = y, q[m][j++]=z, printf(" %d%d%d ",x,y,z);
      r[m]=j/3;
    }
    printf(" : r=%d\n",r[m]);
  }

  // Set count to 1 representing the empty box (symmetry c3v)
  c[8]=1; c[3]=1;

  // Start searching at f=1, with 0 cells occupied and symmetry 3=c3v
  f(1,0,3);

  // c[2 and 6] only contain reflections in y axis, therefore must be multiplied by 3.
  // Similarly the reflections in x and z axis must be subtracted from c[0] and c[4].
  c[0]-=c[2]*2; c[2]*=3;
  c[4]-=c[6]*2; c[6]*=3;

  int cr[9];cr[8]=0;
  printf("non self-complement                   self-complement\n");
  printf("c1  %9d/12=%9d           C1  %9d/6=%9d\n",   c[0], cr[0]=c[0]/12,     c[4], cr[4]=c[4]/6);
  if(cr[0]*12!=c[0])puts("c1 division error");if(cr[4]*6!=c[4])puts("C1 division error");

  printf("c3  %9d/4 =%9d           C3  %9d/2=%9d\n",   c[1], cr[1]=c[1]/4,      c[5], cr[5]=c[5]/2);
  if(cr[1]*4!=c[1])puts("c3 division error");if(cr[5]*2!=c[5])puts("C3 division error");

  printf("cs  %9d/6 =%9d           CS  %9d/3=%9d\n",   c[2], cr[2]=c[2]/6,      c[6], cr[6]=c[6]/3);
  if(cr[2]*6!=c[2])puts("cs division error");if(cr[6]*3!=c[6])puts("CS division error");

  printf("c3v %9d/2 =%9d           C3V %9d/1=%9d\n",   c[3], cr[3]=c[3]/2,      c[7], cr[7]=c[7]);
  if(cr[3]*2!=c[3])puts("c3v division error");

  for(int i=8;i--;)cr[8]+=cr[i];
  printf("total =%d unique =%d",c[8],cr[8]);
}

Output
The program generates an output table of 8 entries, in accordance with the 8 symmetries of the solid. The solid can have any of 4 symmetries as follows (Schoenflies notation)
c1: no symmetry
c3: 3-fold axis of rotation (produces 3-fold axis of rotation in hexagon tiling)
cs: plane of reflection (produces line of reflection in hexagon tiling)
c3v both of the above (produces 3-fold axis of rotation and three lines of reflection through the hexagon corners)

Additionally, when the solid has exactly half the cells with 1's and half with 0's, there exists the possibility of flipping all the 1's and 0's, then inverting the coordinates through the centre of the cube space. This is what I am calling self-complement, but a more mathematical term would be "antisymmetric with respect to a centre of inversion."
This symmetry operation gives a 2-fold axis of rotation in the hexagon tiling.
Patterns which have this symmetry are listed in a separate column. They only occur where N is even.
My count seems to be slightly off for N=4. In discussion with Peter Taylor it appears that I am not detecting tilings which only have a symmetry of a line through the hexagon edges. This is presumably because I have not tested for self complement (antisymmetry) for operations other than (inversion)x(identity.) Testing for self complement for the operatons (inversion)x(reflection) and (inversion)x(3-fold rotation) may uncover the missing symmetries. I would then expect the first line of the data for N=4 to look like this (16 less in c1 and 32 more in C1):
c1   224064/12=18672          C1  534/6=89

This would bring the totals in line with Peter's answer and https://oeis.org/A066931/a066931.txt
current output is as follows.
N=1
non self-complement     self-complement
c1      0/12= 0           C1  0/6= 0
c3      0/4 = 0           C3  0/2= 0
cs      0/6 = 0           CS  0/3= 0
c3v     2/2 = 1           C3V 0/1= 0
total =2 unique =1

non self-complement     self-complement
N=2
c1      0/12= 0           C1  0/6= 0
c3      0/4 = 0           C3  0/2= 0
cs     12/6 = 2           CS  3/3= 1
c3v     4/2 = 2           C3V 1/1= 1
total =20 unique =6

N=3
non self-complement     self-complement
c1    672/12=56           C1  0/6= 0
c3      4/4 = 1           C3  0/2= 0
cs    288/6 =48           CS  0/3= 0
c3v    16/2 = 8           C3V 0/1= 0
total =980 unique =113

N=4 (errors as discussed)
non self-complement     self-complement
c1   224256/12=18688          C1  342/6=57
c3       64/4 =16             C3  2/2= 1
cs     8064/6 =1344           CS  54/3=18
c3v      64/2 =32             C3V 2/1= 2
total =232848 unique =20158

N=5
non self-complement     self-complement
c1  266774112/12=22231176        C1  0/6= 0
c3       1100/4 =275             C3  0/2= 0
cs     451968/6 =75328           CS  0/3= 0
c3v       352/2 =176             C3V 0/1= 0
total =267227532 unique =22306955

To-do list (updated)

Tidy up current code.

Done, more or less

Implement symmetry checking for current layer, and pass a parameter for the symmetry of the previous layer (no point in checking if the last layer was asymmetrical.)

Done, results for odd N agree with published data

Add an option to suppress counting asymmetrical figures (should run much faster)

This can be done by adding another condition to the recursion call: if(s1 && m<n*3-2)f(m + 1,e+d,s1) It reduces the run time for N=5 from 5 minutes down to about a second. As a result the first line of the output becomes total garbage (as do the overall totals) but if the total is already known from OEIS the number of asymmetric tilings can be reconstituted, at least for odd N.
But for even N, the number of asymmetric (according to c3v symmetries) solids that are self-complement would be lost. For this case, a separate program dedicated to solids with exactly (N**3)/2 cells with a 1 may be useful. With this available (and counting correctly) it may be possible to try N=6, but it will take a long time to run.

Implement counting of cells to reduce search to up to (N^3)/2 cubes.

Not done, savings expected to be marginal

Implement symmetry (complementary solid) checking for patterns containing exactly (N^3)/2 cubes.

Done, but seems to have omissions, see N=4.

Find a way to pick the lexically lowest figure from an asymetrical one.

Savings are not expected to be that great. Supressing asymmetrical figures eliminates most of this. The only reflection that is checked is the plane through the y axis (x and z are calculated later by multiplying by 3.) Figures with only rotational symmetry are counted in both their enantiomeric forms. Perhaps it would run nearly twice as fast if only one was counted.

To facilitate this, possibly improve the way the coordinates in each layer are listed (they form degenerate groups of 6 or 3, with possibly a group of 1 in the exact centre of the layer.)

Interesting but probably there are other questions on the site to explore.

Rust, 0.03 s for n=0..=15
Port of @Peter Taylor's Java answer in Rust.
Run it on RustPlayground!

src/main.rs
use std::collections::HashMap;
use num_bigint::BigUint;
use num_traits::{Zero, One, Pow};

fn main() {
    for n in 0..16 {
        println!("{}\t{}", n, count(n));
    }
}

fn count(n: usize) -> BigUint {
    if n == 0 { return BigUint::one() }

    let mut minp = vec![0; 3 * n];

    for index in 2..3 * n {
        minp[index] = (2..=index).find(|p| index % p == 0).unwrap_or(2);
    }

    let e = count_e(n, &minp);
    // println!("n={:?}  e={:?}", n, e);
    let r2 = count_r2(n, &minp);
    // println!("n={:?}  r2={:?}", n, r2);
    let f = count_f(n, &minp);
    // println!("n={:?}  f={:?}", n, f);
    let r3 = count_r3(n, &minp);
    // println!("n={:?}  r3={:?}", n, r3);
    let r = count_r(n, &minp);
    // println!("n={:?}  r={:?}", n, r);
    let fr = count_fr(n, &minp);
    // println!("n={:?}  fr={:?}", n, fr);

    let sum = e
        + r3
        + (r2 + r.clone()) * BigUint::from(2u8)
        + (f + fr) * BigUint::from(3u8);

    sum / BigUint::from(12u8)
}

fn count_e(n: usize, minp: &[usize]) -> BigUint {
    let mut w = vec![0; 3 * n];

    for i in 0..n {
        for j in (i + 1)..=(i + n) { w[j] -= 1 }
        for j in (i + n + 1)..=(i + 2 * n) { w[j] += 1 }
    }

    // println!("w={:?}  minp={:?}",w, minp);

    power_prod(&mut w, minp)
}

fn count_r2(n: usize, minp: &[usize]) -> BigUint {
    let mut w = vec![0; 3 * n];

    for i in 0..n {
        w[3 * i + 2] += 1;
        for j in (3 * i + 1)..=(2 * i + n + 1) { w[j] -= 1 }
        for j in (2 * i + n + 1)..=(i + 2 * n) { w[j] += 1 }
    }

    power_prod(&mut w, minp)
}

fn count_f(n: usize, minp: &[usize]) -> BigUint {
    let mut w = vec![0; 3 * n];

    for i in 0..n {
        for j in (2 * i + 1)..=(2 * i + n) { w[j] -= 1 }
        for j in (i + n + 1)..=(i + 2 * n) { w[j] += 1 }
    }

    power_prod(&mut w, minp)
}

fn count_r3(n: usize, minp: &[usize]) -> BigUint {
    if n & 1 == 1 { return BigUint::zero() }
    let res = count_e(n / 2, minp);
    res.pow(2u8)
}

fn count_r(n: usize, minp: &[usize]) -> BigUint {
    if n & 1 == 1 { return BigUint::zero() }
    let m = n / 2;
    // let mut w = vec![0; 3 * m - 1];
    let mut w = vec![0; 3 * m - 1];

    for i in 0..m {
        for j in 1..=(3 * i + 1) { w[j] += 2 }
        for j in 1..=(i + m) { w[j] -= 2 }
    }

    // power_prod(&w, minp)
    power_prod(&mut w, minp)
}

fn count_fr(n: usize, minp: &[usize]) -> BigUint {
    if n & 1 == 1 { return BigUint::zero() }
    let m:i32 = n as i32 / 2;
    let mut w = vec![0; 3 * n - 2];
    // println!("n={:?}   m={:?}", n, m);

    for j in 1..=m { w[j as usize] -= 1 }
    for j in (2 * m)..(3 * m) { w[j as usize] += 1 }
    for i in 0..=(2 * m - 3) {
        // println!("start{:?}    end{:?}", (i + 2 * m + 1),  (i + 4 * m));
        for j in (i + 2 * m + 1)..=std::cmp::min(w.len() as i32 - 1, i + 4 * m) { w[j as usize] += 1 }
        // println!("start{:?}    end{:?}", (2 * i + 3),  (2 * i + 2 * m + 2));
        for j in (2 * i + 3)..=std::cmp::min(w.len() as i32 - 1, (2 * i + 2 * m + 2)) { w[j as usize] -= 1 }
        // println!("i={:?}",i);
    }

    power_prod(&mut w, minp)
}

fn power_prod(w: &mut [i32], minp: &[usize]) -> BigUint {
    let mut result = BigUint::one();
    // println!("left={:?}   right={:?}", 2, w.len());
    for x in (2..w.len()).rev() {
        if w[x] != 0 {
            let p = minp[x];
            if p == x {
                result *= BigUint::from(p as u8).pow(w[p] as u32);
            } else {
                w[p] += w[x];
                w[x / p] += w[x];
            }
        }
        // println!("x={:?}", x);
    }
    result
}

Cargo.toml
[package]
name = "rust_hello"
version = "0.1.0"
edition = "2021"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[dependencies]
num-bigint = "0.3"
num-traits = "0.2"

[profile.release]
lto = "y"
panic = "abort"

Build and run
$ cargo build --release
$ cargo run --release

    Finished release [optimized] target(s) in 0.03s
     Running `target\release\rust_hello.exe`
0       1
1       1
2       6
3       113
4       20174
5       22306955
6       123222909271
7       3283834214485890
8       421263391026827547540
9       260028731850596651411721718
10      772086476515163830856527013278243
11      11025620741283840573496993339545350520150
12      757129347300072898736973484532998417574513923224
13      249988147151001673705019407088899436686198854480028470000
14      396838668136562149340506906750722486726178589993407032745345886384
15      3028477458789741699983305448780147795991274413836547758038383842702471625120

Scala
A port of @Peter Taylor's Java answer in Scala.
TIO
import scala.collection.mutable.ArrayBuffer
import scala.math._
import java.math.BigInteger

object OptimisedCounter extends App {
  private var minp = ArrayBuffer.fill(2)(0)

  if (args.length > 0) {
    args.map(arg => println(count(arg.toInt)))
  } else {
    (0 until 16).foreach(n => println(s"$n\t${count(n)}"))
  }

  private def count(n: Int): BigInteger = {
    if (n == 0) return BigInteger.ONE

    if (minp.length < 3 * n) {
      val wider = ArrayBuffer.fill(3 * n)(0)
      for (index <- 0 until (minp.length) by 1) wider(index) = minp(index)
            //println(minp.mkString(","))
      (minp.length until wider.length).foreach { x =>
        // Find the smallest prime which divides x
        wider(x) = (2 to x).find(x % _ == 0).getOrElse(2)
      }
      minp = wider

            //println(minp.mkString(","))
    }

    val E = countE(n)
    val R2 = countR2(n)
    val F = countF(n)
    val R3 = countR3(n)
    val R = countR(n)
    val FR = countFR(n)
    val sum = E
      .add(R3)
      .add(R2.add(R).multiply(BigInteger.valueOf(2)))
      .add(F.add(FR).multiply(BigInteger.valueOf(3)))
    sum.divide(BigInteger.valueOf(12))
  }

  private def countE(n: Int): BigInteger = {
    val w = ArrayBuffer.fill(3 * n)(0)
    (0 until n).foreach { i =>
      ((i + 1) to (i + n)).foreach(j => w(j) -= 1)
      ((i + n + 1) to (i + 2 * n)).foreach(j => w(j) += 1)
    }
    powerProd(w.toArray)
  }

  private def countR2(n: Int): BigInteger = {
    val w = ArrayBuffer.fill(3 * n)(0)
    (0 until n).foreach { i =>
      w(3 * i + 2) += 1
      ((3 * i + 1) to (2 * i + n + 1)).foreach(j => w(j) -= 1)
      ((2 * i + n + 1) to (i + 2 * n)).foreach(j => w(j) += 1)
    }
    powerProd(w.toArray)
  }

  private def countF(n: Int): BigInteger = {
    val w = ArrayBuffer.fill(3 * n)(0)
    (0 until n).foreach { i =>
      ((2 * i + 1) to (2 * i + n)).foreach(j => w(j) -= 1)
      ((i + n + 1) to (i + 2 * n)).foreach(j => w(j) += 1)
    }
    powerProd(w.toArray)
  }

  private def countR3(n: Int): BigInteger = {
    if ((n & 1) == 1) return BigInteger.ZERO
    countE(n / 2).pow(2)
  }

  private def countR(n: Int): BigInteger = {
    if ((n & 1) == 1) return BigInteger.ZERO
    val m = n / 2
    val w = ArrayBuffer.fill(3 * m - 1)(0)
    (0 until m).foreach { i =>
      (1 to (3 * i + 1)).foreach(j => w(j) += 2)
      (1 to (i + m)).foreach(j => w(j) -= 2)
    }
    powerProd(w.toArray)
  }

  private def countFR(n: Int): BigInteger = {
    if ((n & 1) == 1) return BigInteger.ZERO
    val m = n / 2
    val w = ArrayBuffer.fill(3 * n - 2)(0)
    (1 to m).foreach(j => w(j) -= 1)
    (2 * m until 3 * m).foreach(j => w(j) += 1)
    (0 to 2 * m - 3).foreach { i =>
      ((i + 2 * m + 1) to (i + 4 * m)).foreach(j => w(j) += 1)
      ((2 * i + 3) to (2 * i + 2 * m + 2)).foreach(j => w(j) -= 1)
    }
    powerProd(w.toArray)
  }

  private def powerProd(w: Array[Int]): BigInteger = {
    var result = BigInteger.ONE
    (2 until w.length).reverse.foreach { x =>
      if (w(x) == 0) ()
      else {
        val p = minp(x)
        if (p == x) result = result.multiply(BigInteger.valueOf(p).pow(w(p)))
        else {
          // Redistribute it. This should ensure we avoid negatives.
          w(p) += w(x)
          w(x / p) += w(x)
        }
      }
    }
    result
  }
}


text:
xkcd-Style Page Numbering
Randall Munroe's book "xkcd, volume 0" uses a rather odd number system for the page numbers. The first few page numbers are
1, 2, 10, 11, 12, 20, 100, 101, 102, 110, 111, 112, 120, 200, 1000, 1001, ...

This looks a bit like ternary, but notice that he skips from 20 straight to 100, from 120 to 200 and from 200 to 1000. One way to define this sequence is to say that it enumerates all ternary numbers that contain at most one 2 and no 1 after that 2. You can find this on OEIS in entry A169683. This number system is known as skew binary.
Your task is to find the representation of a given positive integer N in this number system.
You may write a program or function, taking input via STDIN (or closest alternative), command-line argument or function argument and outputting the result via STDOUT (or closest alternative), function return value or function (out) parameter.
Output may be a string, a number with a decimal representation equal to the skew binary representation, or a list of digits (as integers or characters/strings). You must not return leading zeroes.
This is code golf, so the shortest answer (in bytes) wins.
Fun fact: There is actually some merit to this number system. When incrementing a number, you will always change at most two adjacent digits - you'll never have to carry the change through the entire number. With the right representation that allows incrementing in O(1).
Test Cases
1 => 1
2 => 2
3 => 10
6 => 20
7 => 100
50 => 11011
100 => 110020
200 => 1100110
1000 => 111110120
10000 => 1001110001012
100000 => 1100001101010020
1000000 => 1111010000100100100
1048576 => 10000000000000000001

1000000000000000000 => 11011110000010110110101100111010011101100100000000000001102

I will give a bounty to the shortest answer which can solve the last test case (and any other input of similar magnitude, so don't think about hardcoding it) in less than a second.
Leaderboards
Here is a Stack Snippet to generate both a regular leaderboard and an overview of winners by language.
To make sure that your answer shows up, please start your answer with a headline, using the following Markdown template:
# Language Name, N bytes

where N is the size of your submission. If you improve your score, you can keep old scores in the headline, by striking them through. For instance:
# Ruby, <s>104</s> <s>101</s> 96 bytes

I've had that book since it first came out and never noticed the page numbering.

@AlexA. I have it on my Kindle where you can't do any interesting page numbering.

@LegionMammal978: Tragic.

Is a newline separated list of digits permitted?

@isaacg no that seems like a bit of a stretch.

Is it permissible to take input in binary?

@SuperJedi224 Consensus seems to be no, so sorry (I'd make an exception from consensus if that was the only sort of input your language could handle, but that doesn't seem to be the case).

This reminds me of how I used to count when I was 3. I reasoned: "there's no fifty-ten, so after one-hundred-nine must be two-hundred." I only realized my mistake upon seeing the difference between 59->60 and 109->110, with the extra 0.

Pyth, 17 bytes
Ljb3ye.f!-P-yZ01Q

This is somewhat ridiculously slow - O(output^log_2(3)). It's exponential in the length of the input, but not doubly exponential, like some of the answers on the page. Some ideas taken from @Dennis's answer, here.
Demonstration.
It makes use of .f, Pyth's "loop until n matches have been found" function.

CJam, 24 23 22 21 19 bytes
ri_)2b,,W%{2\#(md}/

This is an O(log n) approach, where n is the input, that completes the last test case instantly. It converts n directly to skew binary, using modular division by the values of the digit 1.
This code finishes with an error that goes to STDERR with the Java interpreter, which is allowed according to the consensus on Meta.
If you try this code in the CJam interpreter, just ignore everything but the last line of output.
The error can be eliminated, at the cost of 2 bytes, by prepending 2> to W%.
Thanks to @MartinBüttner for golfing off one byte.
Background
The skew binary representation ak...a0 corresponds to the integer n = (2k+1-1)ak+...+(21-1)a0.
Since both (2k-1) + ... + (21-1) = 2k+1 - (k + 2) and (2k-1) + ... + 2(2j-1) = 2k+1 - (2j+1 - 2j + k + 1) are less than 2k+1-1, the values of ak to a0 can be recovered by successive modular division by 2k+1-1, 2k-1, etc.
To begin, we have to find the value of 2k+1-1 first. Since n is at most 2(2k+1-1), the integer n + 1 must be strictly smaller than 2k+2.
Thus, taking the integer part of the binary logarithm of n + 1 yields k + 1.
Finally, we observe that the integer n + 1 has ⌊log2(n + 1)⌋ digits in base 2.
How it works
ri    e# Read an integer N from STDIN.
2b,   e# Push the number of binary digits of N + 1, i.e, K + 2 = log(N + 1) + 1.
,W%   e# Push the range [K+1 ... 0].
{     e# For each J in the range:
  2\# e#   J -> 2**J
  (   e#     -> 2**J - 1
  md  e#   Perform modular division of the integer on the stack (initially N)
      e#   by 2**J - 1: N, 2**J - 1 -> N/(2**J - 1), N%(2**J - 1)
}/    e#

In the last two iterations, we perform modular division by 1 and 0. The first pushes an unwanted 0 on the stack. The last attempts to execute 0 0 md, which pops both unwanted 0s from the stack, exits immediately instead of pushing anything and dumps the stack to STDOUT.

Python 2, 67 bytes
def f(n):x=len(bin(n+1))-3;y=2**x-1;return n and n/y*10**~-x+f(n%y)

Seems to work for the given test cases. If I've got this right, this should be about O(place values set in output), so it does the last case with ease.
Call like f(100). Returns a decimal representation equal to the skew binary.
Python 3, 65 bytes
def g(n,x=1):*a,b=n>x*2and g(n,x-~x)or[n];return a+[b//x,b%x][:x]

Slightly less efficient but still logarithmic, so the last case is near-instant.
Call like g(100). Returns a list of digits.

does 2and compile in 3? I'm in 2 and 2and2 throws a syntax error

@TankorSmash 2and2 wouldn't work because it'd get parsed as 2 and2 - try 2and 2, which should work if your Python version is new enough (tested in Python 2.7.10)

Oh nice, you're right. Even on 2.7.3 it works.

CJam, 22 21 20 bytes
ri_me,3fb{0-W<1-!},=

This is an O(enn) approach, where n is the input. It lists the first ⌊en⌋ non-negative integers in base 3, eliminates those that have 2s or 1s after the first 2 (if any) and selects the n+1th.
Try it online in the CJam interpreter.
How it works
ri    e# Read an integer N from STDIN.
_me,  e# Push [0 ... floor(exp(N))-1].
3fb   e# Replace each integer in the range by the array of its digits in base 3.
{     e# Filter; for each array A:
  0-  e#   Remove all 0's.
  W<  e#   Remove the last element.
  1-  e#   Remove all 1's.
  !   e#   Logical NOT. Pushes 1 iff the array is empty.
},    e#   If ! pushed 1, keep the array.
=     e# Select the (N+1)th element of the filtered array.

Pyth, 20 bytes
Jt^2hslhQ#p/=%QJ=/J2

Runs in O(log(input())), well under a second for the final test case. Based around a run until error loop. No trailing newline.
Demonstration.
Explanation:
Jt^2hslhQ#p/=%QJ=/J2
                        Implicit: Q is the input.
      lhQ                          log(Q+1,2)
     slhQ                    floor(log(Q+1,2))
    hslhQ                    floor(log(Q+1,2))+1
  ^2hslhQ                 2^(floor(log(Q+1,2))+1)
 t^2hslhQ                 2^(floor(log(Q+1,2))+1)-1
Jt^2hslhQ               J=2^(floor(log(Q+1,2))+1)-1
         #              until an error is thrown:
            =%QJ        Q=Q%J
                =/J2    J=J/2
           /            The value Q/J, with the new values of Q and J.
          p             print that charcter, with no trailing newline.

J is initialized to the value of the smallest skew binary digit position that is larger than the input. Then, each time through the loop, we do the following:

Remove each digit of value J from Q with =%QJ. For instance, if Q=10 and J=7, Q becomes 3, which corresponds to the skew binary changing from 110 to 10. This has no effect in the first iteration.
Change J to the next smaller skew binary base value with =/J2. This is floored division by 2, changing J=7 to J=3, for instance. Because this happens before the first digit is output, J is intialized one digit position higher than needed.
Find the actual digit value with /QJ (effectively).
Print that value with p, instead of Pyth`s default printing, to avoid the trailing newline.

This loop will repeat until J becomes zero, at which point a divide by zero error will be thrown, and the loop will end.

ES6, 105 bytes
f=n=>{for(o=0;n--;c?o+=Math.pow(3,s.length-c):o++)s=t(o),c=s.search(2)+1;return t(o);},t=a=>a.toString(3)

Usage is:
f(1048576) => `"10000000000000000001"
Test the last argument at your own peril. I gave up after 5 seconds.
And pretty print with comments!
f=n=>{ //define function f with input of n (iteration)
    for(o=0; //define o (output value in decimal)
        n--; //decrement n (goes towards falsy 0) each loop until 0
        c?o+=Math.pow(3,s.length-c):o++) //if search finds a 2, increment output by 3^place (basically moves the 2 to the left and sets the place to 0), else ++
        s=t(o), //convert output to base 3
        c=s.search(2)+1; //find the location of 2, +1, so a not-found becomes falsy 0.
    return t(o); //return the output in base 3
},

t=a=>a.toString(3);  //convert input a to base 3

Btw, unnamed functions are perfectly acceptable, so you don't need f=.

-16 bytes: f=n=>{for(o=0;~n--;o+=c?Math.pow(3,s.length+c):1)s=o.toString(3),c=~s.search(2);return s}

@nderscore Pretty cool :D

-7 bytes if you use ES7: replace Math.pow(3,s.length+c) with 3**(s.length+c).

@GustavoRodrigues I haven't even finished learning ES6! @_@

Those technologies are still being developed, so it's normal. You may like this table, there are many new features from ES6 and ES7, links to the specifications and also where those are supported.

I think they should wait until ES6 is more widely supported first.

Retina, 55 bytes
^
0a
(+`12(.*a)1
20$1
0?2(.*a)1
10$1
0a1
1a
)`1a1
2a
a
<empty line>

Takes input in unary.
Each line should go to its own file but you can run the code as one file with the -s flag. E.g.:
> echo 11111|retina -s skew
12

Method: Executes incrementation on a string input number times starting from the string 0.
We use the following incrementation rules:

if contains 2: ^2 -> ^12; 02 -> 12; 12 -> 20
if doesn't contain 2: 0$ -> 1$; 1$ -> 2$

(There can be at most one 2 in the string; ^ and $ marks the start and end of the string in the rules.)
More information about Retina.

Java, 154 148
n->{String s="0";for(;n-->0;)s=s.contains("2")?s.replaceAll("(^|0)2","10").replace("12","20"):s.replaceAll("1$","2").replaceAll("0$","1");return s;}

This answer take the form of a single anonymous function that takes an integer argument and returns the answer as a String. Below is a complete class for testing this solution.
import java.util.function.Function;
public class Skew {
    public static void main(String[] args){
        Function<Integer,String> skew = n->{String s="0";for(;n-->0;)s=s.contains("2")?s.replaceAll("(^|0)2","10").replace("12","20"):s.replaceAll("1$","2").replaceAll("0$","1");return s;};

        for(String s:args){
            System.out.println(skew.apply(Integer.parseInt(s)));
        }
    }
}

Java, 337 335 253 246 244 bytes
A method that takes in the index as a long and returns the result as a string
Uses a long for the index so can theoretically handle the last test case, but I really wouldn't suggest it.
String f(long i){List<Long>l=new ArrayList<>();l.add(0L);for(;i-->0;){int j=l.indexOf(2);if(j!=-1){l.set(j,0L);if(j==0){l.add(0,1L);}else{l.set(j-1,l.get(j-1)+1);}}else{j=l.size()-1;l.set(j,l.get(j)+1);}}String s="";for(long q:l)s+=q;return s;}

You could shorten this quite a bit by making it a function instead of a full program (and taking input as an argument instead of from a Scanner).

You don't need the braces in the if(j == 0)  statement (four bytes). You don't need to declare k; you can just use j again (four more). You can use generic type inference (in Java 7) in your List declaration (new ArrayList<>();) (another 4)

Bash + coreutils, 52
dc -e3o0[r1+prdx]dx|grep -v 2.\*[12]|sed -n $1{p\;q}

This is a brute-forcer, so its rather slow for larger numbers.
Output:
$ ./xkcdnum.sh 1000
111110120
$

CJam, 24 bytes
Q{0\+2a/())+a\+0a*}ri*si

This is an O(n log n) approach, where n is the input. It starts with the skew binary representation of 0 and increments the corresponding integer n times.
Try it online in the CJam interpreter.
Background
Incrementing a number in skew binary can be done by following two easy steps:

Replace an eventual 2 by a 0.
If a 2 has been replaced, increment the digit to its left.
Otherwise, increment the last digit.

How it works
Q     e# Push an empty array.
{     e# Define an anonymous function:
  0\+ e#   Prepend a 0 to the array.
  2a/ e#   Split the array at 2's.
  (   e#   Shift out the first chunk of this array.
  )   e#   Pop the last digit.
  )+  e#   Increment it and append it to the array.
  a\+ e#   Prepend the chunk to the array of chunks.
  0a* e#   Join the chunks, using [0] as separator.
      e#   If there was a 2, it will get replaced with a 0. Otherewise, there's
      e#   only one chunk and joining the array dumps the chunk on the stack.
}     e#
ri*   e# Call the function int(input()) times.
si    e# Cast to string, then to integer. This eliminates leading 0's.

Haskell, 73 72
Thanks to @nimi for 1 byte!
This solution won't be winning any bounty, the last couple tests take an inordinate amount of time to run, but, I think I've golfed it fairly well.
i(2:b)=1:0:b
i[b]=[b+1]
i(b:2:c)=b+1:0:c
i(b:c)=b:i c
s=(iterate i[0]!!)

This solution is a rather naïve approach that calculates the skew binary number n by incrementing 0 n times.

VBA, 209 147 142 bytes
Sub p(k)
For i=1To k
a=StrReverse(q)
If Left(Replace(a,"0",""),1)=2Then:q=q-2*10^(InStr(a,2)-1)+10^InStr(a,2):Else:q=q+1
Next
msgbox q
End Sub

My math is inefficient and my golfing could use work. But it's my first attempt at PoG and thought I'd give this one a try. Sort of a Brute Force way though.
It's just counting by 1 unless the last digit is a 2 then steps back by 2 and jumps forward 10. Plus the trailing 0's.
This stops working at 65534 because VBA insists on giving output in scientific notation, but the logic should work fine for even higher numbers
Looking forward to golfing suggestions, VBA is not very golf-friendly but it's not often represented, and I think it can beat Java for length.
Edit1:
Thanks Manu for helping shave off 62 bytes
Edit2: Swapped debug.print for msgbox as output. Saved 5 bytes

You can remove the brackets from Debug.Print (q). Also, you can remove most of the whitespace (the editor will put them back, but they aren't necessary). You don't need to declare k as Long, simply write k. It will be a variable of the type Variant and the code will still work. With these tips you should get down to ~165 bytes.

Some more thoughts: You can leave out the first and last argument of InStr, they are optional. Trim() isn't necessary, since you don't have whitespace. Applied correctly, I get to 147 bytes.

Thanks for the Help Manu, One quick question, Output should be the Standard Output. I'm unsure what that would be in VBA. debug.print q would be the Standard output? msgbox q is shorter but Again seems like its not quite the standard output.  Sheet1.cells(1,1) seems like the Typical output but assumes its running in excel. I'm just not quite sure on how strict code-golf is about these sort of things.

Congratulations, you beat the java answer ;) I don't know either... Just use MsgBox, if someone complains, you can still change it back.

Javascript ES6, 99 86 78 76 72 chars
f=n=>{for(s="1";--n;s=s.replace(/.?2|.$/,m=>[1,2,10][+m]||20));return s}

// Old version, 76 chars:
f=n=>{for(s="1";--n;s=s.replace(/02|12|2|.$/,m=>[1,2,10][+m]||20));return s}

// Old version, 86 chars:
f=n=>{for(s="1";--n;s=s.replace(/(02|12|2|.$)/,m=>[1,2,10,,,,,,,,,,20][+m]));return s}

// Old version, 99 chars:
f=n=>{for(s="1";--n;s=s.replace(/(^2|02|12|20|.$)/,m=>({0:1,1:2,2:10,12:20,20:100}[+m])));return s}

Test:
;[1,2,3,6,7,50,100,200,1000,10000,100000,1000000,1048576].map(f) == "1,2,10,20,100,11011,110020,1100110,111110120,1001110001012,1100001101010020,1111010000100100100,10000000000000000001"

Fun fact: There is actually some merit to this number system. When incrementing a number, you will always change at most two adjacent digits - you'll never have to carry the change through the entire number. With the right representation that allows incrementing in O(1).

Thank you for the fact - it is the base of my solution :)

How could I leave unnecessary braces in the regex? o_O

Octave, 107 101 bytes
Should be O(log n) if I'm figuring this right...
function r=s(n)r="";for(a=2.^(uint64(fix(log2(n+1))):-1:1)-1)x=idivide(n,a);r=[r x+48];n-=x*a;end;end

Pretty-print:
function r=s(n)
  r="";
  for(a=2.^(uint64(fix(log2(n+1))):-1:1)-1)
    x=idivide(n,a);
    r=[r x+48];
    n-=x*a;
  end
end

I was somewhat stymied addressing the final challenge, since Octave by default treats everything as floating point numbers and I didn't have the necessary precision to calculate the last one. I got around that by spending precious bytes to force everything to be an unsigned integer. The result of the last one was too big to treat as a number, so the result is a string.
Output (I'm including 1e18 - 1 to show I can do that accurately, and the last set of outputs shows how long it takes to calculate that value):
octave:83> s(uint64(1e18))
ans = 11011110000010110110101100111010011101100100000000000001102

octave:84> s(uint64(1e18)-1)
ans = 11011110000010110110101100111010011101100100000000000001101

octave:85> tic();s(uint64(1e18)-1);toc()
Elapsed time is 0.0270021 seconds.

Turing Machine Code, 333 293 bytes
I'm using an encoding as used here.
This machine uses 9 states and 11 colors.
If binary input is permissible, this can be reduced to a mere 4 colors, saving a few tens of bytes in the process.
0 _ _ l 1
0 * * r 0
1 9 8 l 2
1 8 7 l 2
1 7 6 l 2
1 6 5 l 2
1 5 4 l 2
1 4 3 l 2
1 3 2 l 2
1 2 1 l 2
1 1 0 l 2
1 0 9 l 1
1 _ _ r 8
2 _ _ l 3
2 * * l 2
3 _ 1 r 4
3 * * l 5
4 _ _ r 0
4 * * r 4
5 * * l 5
5 _ _ r 6
6 _ _ l 7
6 2 0 l 7
6 * * r 6
7 _ 1 r 4
7 0 1 r 4
7 1 2 r 4
8 _ _ * halt
8 * _ r 8

If the above link isn't working (sometimes it works for me, other times the page refuses to load) you may also test this using this java implementation.

T-SQL, 221 189 177 bytes
EDIT: The original versions of this code would produce incorrect output for some numbers, this has been corrected.
With every query here, just add the number to calculate before the first comma.
Everyone know that T-SQL is the best golfing language. Here is a version that will calculate even the last test case. On the machine I tested it on, it ran in under a second, I'd be interested to see how it runs for everyone else.
DECLARE @ BIGINT=,@T VARCHAR(MAX)='';WITH M AS(SELECT CAST(2AS BIGINT)I UNION ALL SELECT I*2FROM M WHERE I<@)SELECT @T += STR(@/(I-1),1),@%=(I-1)FROM M ORDER BY I DESC SELECT @T

And here it is again, but readable:
DECLARE
    @ BIGINT=,
    @T VARCHAR(MAX)='';

WITH M AS
(
    SELECT
        CAST(2 AS BIGINT) I

    UNION ALL

    SELECT I * 2
    FROM M
    WHERE I < @
)

SELECT
    @T+=STR(@/(I-1),1),
    @%=(I-1)
FROM M
ORDER BY I DESC

SELECT @T

If I only use ints, this can be a bit shorter, coming out at 157 bytes:
DECLARE @ INT=,@T VARCHAR(MAX)='';WITH M AS(SELECT 2I UNION ALL SELECT I*2FROM M WHERE I<@)SELECT @T+=STR(@/(I-1),1),@%=(I-1)FROM M ORDER BY I DESC SELECT @T

And once more again, more readable:
DECLARE
    @ INT=,
    @T VARCHAR(MAX)='';

WITH M AS
(
    SELECT
        2I

    UNION ALL

    SELECT
        I * 2
    FROM M
    WHERE I < @
)

SELECT
    @T+=STR(@/(I-1),1),
    @%=(I-1)
FROM M
ORDER BY I DESC

SELECT @T

Remember @ is a valid identifier in sql and you most likely can get away with Char(8000)  which is still cheaper than nvarchar(max). You can also convert to char instead of varchar, or use the str function.

@MichaelB Oh, I thought that I had used @, silly me. The CHAR(8000) advice is pretty good, I'll try that. I always seem to forget about the existence of STR() thanks for the heads up.

doesn't actually help. However, you can rewrite the part after the CTE into::
select @t=concat(@t,@/i) that should be smaller. Requires sql2012 though.

@MichaelB ah. CONCAT, I'm on 2008. So I can't test it without using an SQL fiddle right now. Good call though.

Perl, 66 bytes
Number is to be input through STDIN.
$_=1;$c=<>;s/(.*)(.?)2(.*)/$1.$2+1 .$3.0/e||$_++while(--$c);print;

Could you explain how your solution works? I don't see how you need (.?) in $2 since (.*) in $1 should be greedy and get that character first. But if it gets removed the code doesn't produce the correct results anymore! By the way, you don't need the final ;.

@CJDennis Thanks for that byte saved. Anyway, the .? will get the digit coming right before the two unless there is no digit there (e.g. 20). In cases such as 120 or 10020, the regex groups like this: ()(1)2(0) and (10)(0)2(0). Then the first group is simply ignored, the second group (which is always either one digit if possible or blank) is incremented, and the third group (always consisting of zeroes) is ignored and a zero is added. I simply used the OEIS entry as my guide for this regex.

I got your code down to 53 bytes: $c=<>;s/(.*)2(.*)/$1+1 .$2.0/e||$_++while($c--);print. I was right, the (.?) never captured anything.

Can be optimized to $_=1;$c=<>;s/(.?)2/1+$1.0/e||$_++while(--$c);print, which is 50 bytes. .* at the beginning or end can be optimized away, if you just replace it with the original text. Also there is no reason to add the 0 at the end, as there are always only zeros in the original $3.

Perl, 84 70 67 bytes
$n=<>;$d*=2while($d++<$n);$_.=int($n/$d)while($n%=$d--,$d/=2);print

Not very golfy, Getting better but it does work very quickly!
Dennis's suggestion gets it down to 51 (50 bytes + -p switch)
$d*=2while$d++<$_;$\.=$_/$d|0while$_%=$d--,$d/=2}{

It must be called like perl -p skew_binary.pl num_list.txt where num_list.txt contains a single line with the number to encode on it.

@frederick we're on the edge of our seats waiting for 2.

@RobertGrant He golfed his comment from two things to one thing!

Two things:

Instead of using $ARGV[0], use <> as input. It will take input from stdin unless there are any files as arguments.
For odd counting schemes like this, use regular expressions. Instead of doing odd math operations, you can replace digits in a number as if it is a string. The best part is that you can use math operations (like increment) on it at the same time, provided that it is a string consisting completely of digits. Check the documentation for the regular expression operators, as they can be very useful in so many occasions.

Sorry, I pressed enter and it saved before I finished the comment.

@frederick don't be so modest. We know what happened!

@frederick I don't understand your solution at all (which uses regex) and I got fed up waiting for the longest one to finish. I am very good with regex. How would you suggest I turn maths into regular expressions?

@CJDennis I used the /e switch on the regex. However, my solution is incapable of efficiently calculating high numbers such as the last test case.

With the -p switch (1 byte), you can use $d*=2while($d++<$_);$\.=$_/$d|0while($_%=$d--,$d/=2)}{.

@Dennis Every time I type a new number it appends it to the previous one. Only the first one gives the correct output.

I'm not sure what you mean. The program is supposed to take a single number as input, no? (While were at it, the parentheses can be eliminated as well.)

@Dennis I was just looking up exactly what the -p switch does. It wraps the program in an endless loop, so to get the right results you have to break out of it every time and rerun it.

The trailing }{ in my proposed code mean that s/0// is not needed (I've tried explaining here how this works out). Eliminating the parentheses as well, $d*=2while$d++<$_;$\.=$_/$d|0while$_%=$d--,$d/=2}{ gets you down to 51 bytes.

OK, I've finally worked out how to run it! perl -p skew_binary.pl num_list.txt

Pyth, 19 bytes
m/=%Qtydtd^L2_SslhQ

Logarithmic complexity. Easily finishes in the necessary amount of time. Output in the form of a list of digits.
Demonstration.

JavaScript (Node.js), 40 bytes
v=>(g=n=>n>v?'':g(n-~n)+(v-(v%=n))/n)(1)

Try it online!
This solution does not work on last testcase. It just cause stack overflow on it.

Mathematica, 65
Should be quite fast enough, although I have to admit I peeked at the other submissions before making this.
f = (n = #;
     l = 0;
     While[n > 0,
      m = Floor[Log2[1 + n]];
      l += 10^(m - 1);
      n -= 2^m - 1
     ]; l)&

Usage:
f[1000000000000000000]

Output:
11011110000010110110101100111010011101100100000000000001102

Starts giving MaxExtraPrecision error messages somewhere past 10^228 (for which it calculates the result in .03 seconds on my machine)
After removing the MaxExtraPrecision limit, it will handle numbers up to around 10^8000 in a second.
Input:
Timing[Block[{$MaxExtraPrecision = Infinity}, f[10^8000]];]

Output:
{1.060807, Null}

C, 95 bytes
void f(unsigned long i,int*b){for(unsigned long a=~0,m=0;a;a/=2,b+=!!m)m|=*b=i/a,i-=a**b;*b=3;}

This accepts an integer and a buffer in which to return the digits.  The results are stored in b, terminated with value 3 (which can't occur in the output).  We don't have to handle input of 0 (as the question specified only positive integers), so there's no special-casing to avoid empty output.
Expanded code
void f(unsigned long i,int*b)
{
    for (unsigned long a=~0, m=0;  a;  a/=2, b+=(m!=0)) {
        *b = i/a;               /* rounds down */
        i -= *b * a;
        m = m | *b;             /* m != 0 after leading zeros */
    }
    *b=3;                       /* add terminator */
}

We operate by successive subtraction, starting with the most significant digit.  The only complication is that we use the variable m to avoid printing leading zeros.  A natural extension to unsigned long long may be made if desired, at the cost of 10 bytes.
Test program
Pass numbers to be converted as command arguments.  It converts the int array buffer to a printable digit string.  Runtime is less than one millisecond for input 1000000000000000000.
#include <stdio.h>
#include <stdlib.h>
int main(int argc, char**argv)
{
    while (*++argv) {
        unsigned long i = strtoul(*argv, NULL, 10);
        int result[1024];
        f(i,result);

        /* convert to string */
        char s[1024];
        {char*d=s;int*p=result;while(*p!=3)*d++=*p+++'0';*d=0;}
        printf("%lu = %s\n", i, s);
    }

    return EXIT_SUCCESS;
}

Test results
$ ./51517 $(seq 20)
1 = 1
2 = 2
3 = 10
4 = 11
5 = 12
6 = 20
7 = 100
8 = 101
9 = 102
10 = 110
11 = 111
12 = 112
13 = 120
14 = 200
15 = 1000
16 = 1001
17 = 1002
18 = 1010
19 = 1011
20 = 1012

I guess a C++ version is similar, but can use auto a=~0ull for a slight advantage...

CoffeeScript, 92 69 bytes
Based of Qwertiy's answer and updates:
f=(n)->s='1';(s=s.replace /(.?2|.$)/,(m)->[1,2,10][+m]||20)while--n;s

# older version, 92 bytes
f=(n)->s='1';(s=s.replace /(^2|02|12|20|.$)/,(m)->{0:1,1:2,2:10,12:20,20:100}[+m])while--n;s

Converting to other language without even simple optimization by removing unnecessary braces in the regex doesn't seem cool for me...

@Qwertiy I've provided attribution to your answer, and with both of the answers I can't produce the same results without the brackets in the regex

The whole occurrence is replaced. Why do you need it to be inside of the group? JS version works in Firefox without brackets.

Japt, 31 bytes
_r/?2|.$/g_÷C ç20 ª°Zs3}}gU['0]

Try it online!
Almost direct port of this JS solution. No idea if there's better way.
Unpacked & How it works
X{Xr/?2|.$/gZ{Z÷C ç20 ||++Zs3}}gU['0]

X{     Declare a function...
Xr       Accept a string, replace the regex...
/?2|.$/g   /.?2|.$/   (g is needed to match *only once*, opposite of JS)
Z{       ...with the function... (matched string will be 0,1,2,02 or 12)
Z÷C        Implicitly cast the matched string into number, divide by 12
ç20        Repeat "20" that many times (discard fractions)
||         If the above results in "", use the next one instead
++Z        Increment the value
s3         Convert to base-3 string (so 0,1,2 becomes 1,2,10)
}}
gU['0] Repeatedly apply the function on "0", U (input) times

Stax, 16 bytes
üxëàè£öΦGΩ│Je5█ò

Run and debug it
I'm not sure exactly what the formal complexity class is, but it's fast enough to do all the test cases in a tenth of a second on this machine.
Unpacked, ungolfed, and commented, it looks like this.  In this program, the x register originally contains the input.
z       push []
{       start block for while loop
 |X:2N  -log2(++x)
 {^}&   increment array at index (pad with 0s if necessary)
 xc:G-X unset high bit of x; write back to x register
w       while; loop until x is falsy (0)
$       convert to string

Run this one

Husk, 16 bytes
!mdfö¬f≠1hfImB3N

Try it online! or Verify first 100 values
Same method as Dennis' CJam answer.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Golf Your Favorite Minecraft Tool
In Minecraft, the default item textures are all reasonably simple 16×16 pixel images, which makes them seem ideal for golfing.
Below are simplified textures of the five "core" diamond tools in Minecraft: pickaxe, shovel, axe, sword, and hoe.
The images shown are enlarged to show their detail. Click on an image to view its correctly sized 16×16 pixel version.

To make golfing easier, I've modified each of them from the originals to only use the five same 24-bit RGB colors:

R=75 G=82 B=73 for the background.
R=51 G=235 B=203 for the diamond tool heads.
R=14 G=63 B=54 for the diamond outlines.
R=137 G=103 B=39 for the wooden handle core.
R=40 G=30 B=11 for the wooden handle outlines.

Choose your favorite tool out of the five and write a program that outputs its simplified 16×16 pixel texture in any common lossless truecolor image format (such as bpm, png, ppm, etc.).
So, for example, if you chose the axe, you would write a program that outputs this image:
No input should be taken and a web connection should not be required. The image can be output as a file with the name of your choice, or the raw image file data can be output to stdout, or you can simply display the image.
You only need to choose one of the five images. The program that outputs any one of the five images in the fewest number of bytes is the winner.
You may write programs for more than one of the images, but only the one with the minimum number of bytes counts towards your score. If there's a tie, the highest voted post wins.

If you enjoy PPCG and play Minecraft, I invite you to come join our trial Minecraft server. Just ask in the dedicated chatroom.

"Enlarged to show detail." I can only think of cereal boxes.

I think there may be some compression issues or you have some typos in the colours. Mathematica claims, these five colours are used (at least for the shovel): [[76 82 73] [26 63 54] [106 234 204] [39 30 13] [133 102 45]]

Similar: http://codegolf.stackexchange.com/questions/39915/pixel-art-episode-1-display-super-mario

@MartinBüttner Are you sure? I re-downloaded the shovel image and made sure that the 5 colors are correct. I also made sure there were exactly 5 colors. imgur may be compressing the enlarged images but those are not the ones you should be using.

@Calvin'sHobbies Is input allowed?

@user41805 Input is never allowed if it isn't explicitly mentioned in the challenge.

CJam, 119 bytes
'P3NGSGN255N]o67T"…vîþáy$<OW¥ÓNZ"256b9b~99T]2/e~W%"LRI?6jêÌ'
…f-":i3/f=F,_W%:)+{)/(\:~}%{G/({)S*S+oNo}%1>\:~+}G*

Try it online!
I've chosen the shovel.
This program prints a PPM file to STDOUT.
The basic idea is to unroll the image along antidiagonals and then use run-length encoding. With this technique, the shovel contains the fewest runs. For reference the number of runs per image (in the order given in the challenge) is:
{60, 26, 38, 43, 37}

This is 129 bytes, not 119.

@Makonede I don't think so?

Well, TIO also claims Pyth uses an SBCS when in reality it does not. If CJam really does use an SBCS, could you link to it?

CJam doesn't so much use an SBCS as allow any byte value in the program @Makonede. These are often shown in older CJam answers, such as this, using extended ASCII characters to avoid using unprintable characters.

Minecraft 18w11a (.mcfunction), 757 bytes
fill ~ ~ ~ ~15 ~ ~15 ice
fill ~13 ~ ~13 ~7 ~ ~11 cyan_wool
fill ~12 ~ ~14 ~10 ~ ~8 cyan_wool
fill ~12 ~ ~13 ~10 ~ ~11 diamond_block
fill ~11 ~ ~12 ~9 ~ ~10 diamond_block
fill ~10 ~ ~11 ~8 ~ ~9 diamond_block
fill ~3 ~ ~4 ~1 ~ ~2 dirt
setblock ~3 ~ ~4 oak_planks
setblock ~2 ~ ~3 oak_planks
clone ~3 ~ ~4 ~1 ~ ~2 ~4 ~ ~5
setblock ~4 ~ ~5 oak_planks
setblock ~4 ~ ~7 ice
setblock ~6 ~ ~5 ice
clone ~6 ~ ~5 ~4 ~ ~7 ~7 ~ ~8
setblock ~9 ~ ~10 diamond_block
setblock ~4 ~ ~4 dirt
setblock ~3 ~ ~5 dirt
setblock ~7 ~ ~7 dirt
setblock ~6 ~ ~8 dirt
setblock ~1 ~ ~2 ice
fill ~12 ~ ~9 ~12 ~ ~8 ice
setblock ~11 ~ ~8 ice
fill ~8 ~ ~13 ~7 ~ ~13 ice
setblock ~7 ~ ~12 ice
fill ~ ~ ~ ~15 ~ ~15 light_gray_concrete replace ice
fill ~ ~ ~ ~9 ~ ~10 dark_oak_bark replace dirt

Of course someone had to answer the question with Minecraft. Place answer inside of a datapack, and run with /function <packname>:<filename>. The shovel is drawn relative to you in the +X and +Z direction. Colors are wrong but I'll count that as a language limitation ;)
But the shovel is actually made out of wood and diamonds!!!!
Output

Answers in languages that postdate the question no longer need to be marked as non-competing

oh hey, that's cool!

I really like the idea of using fill ... replace to save bytes. +1

You can use TNT or air instead of dirt to save 6 bytes.

You can save a lot of bytes by using absolute coordinates.

@xigoi Wow! That's a 76 byte save. Combined with my tip, you save 82 bytes!

@Makonede air can use keep shorter?

@l4m2 Oh yeah! That saves 7 more bytes. 89 bytes saved so far

JavaScript ES6, 353 bytes
document.write(`
(p=parseInt)(e[0],17)+`px ${p(e[1],17)}px 0 #${['33EBCB','0E3F36','896727','281E0B'][e[2]]},`)}9px 9px 0 8px #4B5249"`)

This heavily abuses CSS3 box-shadows to create a pixelized version of the image, in this case the Minecraft hoe. The Stack Snippet below uses ES5 for easy testing and is somewhat ungolfed (You'll have to zoom in to see it well).

s='931a31b31841940a40b40c41951a51b50c50d51e53f52b61c60d60e62f63c73d70e70f71b83c82d83e81a93b92c939a3aa2ba38b39b2ab37c38c29c36d37d28d35e36e27e34f35f26f34g35g3'.replace(/.../g,function(e){
  return parseInt(e[0],17)+'px '+parseInt(e[1],17)+'px 0 #'+['33EBCB','0E3F36','896727','281E0B'][e[2]]+','
})
document.write('<p style="width:1px;height:1px;box-shadow:'+s+'9px 9px 0 8px #4B5249"')

Could this be shortened by using base64 encoding instead of hex?

The largest source image is 297 bytes. Kinda sad to be over that :)

@JB It's [kolmogorov-complexity] for a reason ;)

@ThomasKwa Doesn't base64 actually make it longer? (It's actually base-17 because there's gs near the end.)

I don't know how expensive interpreting base64 encoding is in Javascript, but the string literal itself should be 2/3rds the size in base64 as in hex.

05AB1E, 113 97 bytes
-16 by using RLE (run-length encoding) for the color palette indices.
„P3žvD₅•P{¡æjž²+JGÀ‰íÔñ•Ƶ‡в3ô•N…ý~îc|\мï¶„Ä^€¿^•5в¦•RθÅª`<≠cλS²bÖ"₁Þ@„ÉI¬5øëY4ÀΘØ¥P@û h•29вÅΓè˜»»

Try it online! Beats all other answers. Outputs a plain PPM of a diamond hoe.
Note: The first » can be a ) (or a ¶, ð, or õ) to output the image data separated by spaces: Try it online! However, in the official PPM Format Specification, under Plain PPM, it states:

No line should be longer than 70 characters.

It still renders properly, but is not recommended.
The first » can also be a ` with no change in functionality: Try it online!
Pure GIF compression, 109 107 bytes
•+Ùι1Ù³œ¤∞4ûи„ĆQ´ηÏzεÝ₅¼úΣ¨ü6Ó…Ï,È–ð:<û∞¢Þ§α»«¥<çp:Ù“ζÞ¢ó‡lÿнδßiÌ$}ªiƒ*₅xΛxJ=yh‹×м’’®']!JÍ¦Z¹JQc›Vàÿ¸Ìм•h2ô

Try it online! Outputs as a list of hexadecimal codepoints.
•...•h2ô  # trimmed program
•...•     # 121686424391392225557581524495943776461362082974476220323034665137613008241067075561955081975776948192923502794622971144661356950543760600319115998403445050769039958830830234731371315671091547736084642764295720991029129326629299863698388004372539...
     h    # in hexadecimal...
       ô  # split into pieces of length...
      2   # literal
          # implicit output

Node.JS using jimp, 393 377 376 bytes
377, 376 thanks to Makonede
new(require('jimp'))(16,16,1263684095,(e,i)=>{i.z=i.setPixelColor;i.z(x=673057791,2,12);i.z(x,4,14);i.z(d=871091199,13,3);for(j=0;j<8;j++){i.z(w=2305239039,3+j,13-j);i.z(x,2+j,13-j);i.z(x,3+j,14-j)}for(j=0;j<3;j++){i.z(e=239023871,11+j,2);i.z(e,14,3+j);i.z(e,8+j,5-j);i.z(d,9+j,5-j);i.z(d,10+j,5-j);i.z(d,10+j,6-j);i.z(d,11+j,6-j);i.z(d,11+j,7-j);i.z(e,11+j,8-j)}i.write('')})

I'm not sure if this one is allowed as it uses a separate library for processing the image.
Uses hex colors in decimal, sets i.s to i.setPixelColor to save bytes, sets color variables as it's first used

as long as the question does not disallow libraries, and the library is explicitly stated, it is fine.

However, it is a good idea to tell us which tool you've made.

-49 by removing some semicolons, using an empty file name, using s instead of i.s, and only assigning variables when first needed: new(require('jimp'))(16,16,1263684095,(e,i)=>{(s=i.setPixelColor)(x=673057791,2,12);s(x,4,14);s(d=871091199,13,3);for(j=0;j<8;j++){s(w=2305239039,3+j,13-j);s(x,2+j,13-j);s(x,3+j,14-j)}for(j=0;j<3;j++){s(e=239023871,11+j,2);s(e,14,3+j);s(e,8+j,5-j);s(d,9+j,5-j);s(d,10+j,5-j);s(d,10+j,6-j);s(d,11+j,6-j);s(d,11+j,7-j);s(e,11+j,8-j)}i.write('')}).

@Makonede that's actually really clever assigning variables like that. however i tried using s instead of i.s and jimp threw an error for this being undefined.

actually it seems it's if it's s not i.s or if it's using the assignment in the use or both

also a blank filename shouldn't/doesn't work, however single letter names work, and .png with no name

@mekb While it is true that most operating systems don't support blank filenames, it's OK to assume that a blank filename should work, so it's allowed. :)

ah okay ;;;;;;;

Python 3, 483 bytes
I chose to make the sword
from PIL import Image as IG, ImageColor as IC
s=IG.new('RGB',(16,16))
w='#6b6727'
b='#4b5249'
d='#33ebcb'
a='#0e3f36'
n='#281e0b'
t=b*13+a*3+b*12+a+d*2+a+b*11+a+d*3+a+b*10+a+d*3+a+b+b*9+a+d*3+a+b*2+b*8+a+d*3+a+b*3+b*2+a*2+b*3+a+d*3+a+b*4+b*2+a+d+a+b+a+d*3+a+b*5+b*3+a+d+a+d*3+a+b*6+b*3+a+d+a+d*2+a+b*7+b*4+a+d+a*2+b*8+b*3+n+w+a+d*2+a+b*7+b*2+n+w+n+b+a*2+d+a+b*6+a*2+w+n+b*4+a*2+b*6+a+d+a+b*13+a*3+b*13
s.putdata([IC.getrgb(t[i:i+7]) for i in range(0,len(t),7)])
s.save('s.png','PNG')

here is the output:

I created a string for each color, and combined them to get a string of hexadecimal numbers. Then I used the python image library to convert that string into an image.

You can save 5 bytes by replacing 's.png' with ''. It assumes that the system the program is being run on allows empty file names. This does not apply to any popular OS (AFAIK) but it's OK to make assumptions here.

Actually, you can save 7 bytes by using s.save('.png'). The format argument isn't necessary if a corresponding extension is given, and s.save('.png') is 2 bytes shorter than s.save('','PNG').

You can save another 9 on top of that by changing Image as IG, ImageColor as IC to Image,ImageColor as c, changing IG to Image, changing IC to c, removing the space between (t[i:i+7]) and for, and replacing len(t) with 1792 (the guaranteed length of t).

And yet another 18 saved by 1. fixing some inefficiencies in the definition of t (for example, b+b*9 -> b*10), 2. making t a tuple instead of a string, 3. not using a variable (this method uses t only once) and 4. using all that to allow for a much shorter putdata method: s.putdata(map(c.getrgb,[*[b]*13,a,a,a,*[b]*12,a,d,d,a,*[b]*11,a,d,d,d,a,*[*[b]*10,a,d,d,d,a]*2,*[b]*5,a,a,b,b,b,a,d,d,d,a,*[b]*6,a,d,a,b,a,d,d,d,a,*[b]*8,a,d,a,d,d,d,a,*[b]*9,a,d,a,d,d,a,*[b]*11,a,d,a,a,*[b]*11,n,w,a,d,d,a,*[b]*9,n,w,n,b,a,a,d,a,*[b]*6,a,a,w,n,*[b]*4,a,a,*[b]*6,a,d,a,*[b]*13,a,a,a,*[b]*13])).

(continued from previous comment because of length limitations) A total of 34 bytes saved!

Actually, s.save('') should work instead of s.save('.png'), saving you another 4 for a total of 38 bytes saved so far.

You can save yet another 10 by using import PIL. Current saved bytes: 48! Also, here's a pastebin with all of my golfs so far applied.

-15 more by using Python 3.8's walrus operator and decimal color values: pastebin

-39 more with a ton more *s and a ton more :=s. Saved 102 bytes so far!! Pastebin

-3 more with a few less *s and a few more +s. Current saved bytes: 105. Pastebin

-2 more by showing the image instead of saveing it. Current saved bytes: 107. Pastebin

-2 more by using a few more *s and :=s. Current saved bytes: 109. Pastebin

GIF, 102 bytes
Hexdump:
00000000: 4749 4638 3961 1000 1000 7000 002c 0000  GIF89a....p..,..
00000010: 0000 1000 1000 824b 5249 0e3f 3633 ebcb  .......KRI.?63..
00000020: 281e 0b89 6727 0000 0000 0000 0000 0003  (...g'..........
00000030: 3308 babc f12d 8220 4490 8cd6 8b67 b59d  3....-. D....g..
00000040: b785 1a98 71a2 9981 25da 506d 078f a4e5  ....q...%.Pm....
00000050: 4a03 a1dd 4c3e 4cb0 c883 f073 f082 21c5  J...L>L....s..!.
00000060: 8327 4900 003b                           .'I..;

Renders as:

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Generate boxes!
Your task is to generate boxes using any one ASCII character with respect to the inputs given.
Test Cases
1 1   --> =====
          =   =
          =====

1 2   --> =========
          =   =   =
          =========

2 1   --> =====
          =   =
          =====
          =   =
          =====

2 2   --> =========
          =   =   =
          =========
          =   =   =
          =========

2 5   --> =====================
          =   =   =   =   =   =
          =====================
          =   =   =   =   =   =
          =====================

Input

Input can be taken from one of the following

stdin
Command-line arguments
Function arguments (2 arguments, one for each number)

Input, if taken from stdin or command line arguments, will contain two positive integers, seperated by a space.

The two numbers denote the number of boxes in each column and row

Output

Boxes must be outputted in stdout (or closest equivalent)
Each box should have three horizontal spaces in them

Rules

Both the numbers will be greater than 0, but will not go beyond 1000
Any visible character can be used for outputting the boxes. (as long as they aren't too harsh on the eye!)
You are permitted to write a full program or a function.
There should be no unnecessary characters except an optional trailing newline character.

Scoring
This is code-golf, so the shortest submission (in bytes) wins.

Pyth, 23 bytes
Mjbmsm@"= "&%k4dh*4HhyG

Defines the function g which works as desired.
Demonstration.

Retina, 57 bytes
1(?=.* (1*))
$1#$1#

1(?=(1*#1*#)*$)
=
1
====
#
=#

Takes input in unary with a trailing newline.
Each line should go to its own file and # should be changed to newline in the files. This is impractical but you can run the code as is as one file with the -s flag, keeping the # markers (and changing the trailing newline to # in the input too). You can change the #'s back to newlines in the output for readability if you wish. E.g.:
> echo -n 11 111#|retina -s boxes|tr # '\n'
=============
=   =   =   =
=============
=   =   =   =
=============

Method: 5 single substitution steps. The first two (first 4 lines) creates an 2*N+1 by M grid of ones and the next 3 steps (6 lines) format it into the desired output.
The intermediate strings (delimited by -'s):
11 111
------------------
111
111
111
111
 111
------------------
111
111
111
111
111
------------------
111
=   =   =
111
=   =   =
111
------------------
============
=   =   =
============
=   =   =
============
------------------
=============
=   =   =   =
=============
=   =   =   =
=============

One file per line always seems like a very odd design choice.

@curiousdannii It was made this way so the newline character would be usable in the regex expressions without escaping but it is planned to be changed into a more reasonable 1-file format.

JavaScript (ES6), 83
A function with parameters rows and columns. Using template strings, there are 2 embedded newlines that are significant and counted.
Output via alert (popup).
As alert use a proportional font, we get a better result using a letter with a width similar to the space instead of =.
Test in Firefox using the console to have the real alert output, or run the snippet below.

f=(r,c)=>{for(x=y="|";c--;)x+="||||",y+="   |";for(o=x;r--;)o+=`
${y}
`+x;alert(o)}

// TEST

// redefine alert to avoid that annoying popup during test
alert=x=>O.innerHTML=x

test=_=>([a,b]=I.value.match(/\d+/g),f(a,b))

test()
<input id=I value='2 3'><button onclick="test()">-></button>

Run code snippet does not work for me, just see a box with '2 3' inside and an arrow that seems to suggest that you should click it to get the result. Using the full page button or tryhing a different browser does not change anything.

@DennisJaheruddin , The same happens for me when using chrome (doesn't support ES6). Try firefox.

@DennisJaheruddin Test in Firefox

JavaScript (ES6), 75 73
Using copious repeat calls, we repeat |, then repeat | with trailing spaces, and repeat both of those repeats to make rows:
f=(y,x)=>alert(((s="|"[r="repeat"](x*4)+`|
`)+"|   "[r](x)+`|
`)[r](y)+s)

(Newlines are significant, per edc65's suggestion to use template strings.)
Snippet:

<input id="x" type="range" max="10" min="1" step="1" value="3"><input id="y" type="range" max="10" min="1" step="1" value="2">

Without template strings, without method-name reuse, and with added whitespace:
f=(y,x)=>alert(
    (
      (s="|".repeat(x*4)+"|\n") +
      "|   ".repeat(x)+"|\n"
    ).repeat(y)+s
)

(Using edc65's recommendation to use | for better visual spacing.)

Pip, 29 24 = 23 + 1 bytes
Requires the -n flag. Takes the height and width as command-line args and builds boxes of of 1s.
([1X4m]XbRLa+1)@<v.1R0s

Explanation:
                         a,b are cmdline args; m is 1000; v is -1; s is " " (implicit)
 [1X4m]                  List containing 1111 and 1000
       Xb                String repetition of each element b times
         RLa+1           Repeat the list a+1 times
(             )@<v       Now we have one row too many at the end, so take everything
                         but the last item (equiv to Python x[:-1])
                  .1     Concatenate a 1 to the end of each row
                    R0s  Replace 0 with space
                         Print, joining list on newlines (implicit, -n flag)

This program takes heavy advantage of the fact that strings are numbers and numbers are strings in Pip. (And the three spaces in those boxes happened to be just right to take advantage of the built-in m variable!)
Here's how the list gets built with the input 2 3:
[1111;1000]
[111111111111;100010001000]
[111111111111;100010001000;111111111111;100010001000;111111111111;100010001000]
[111111111111;100010001000;111111111111;100010001000;111111111111]
[1111111111111;1000100010001;1111111111111;1000100010001;1111111111111]
[1111111111111;"1   1   1   1";1111111111111;"1   1   1   1";1111111111111]

And the final output:
C:\> pip.py 2 3 -ne "([1X4m]XbRLa+1)@<v.1R0s"
1111111111111
1   1   1   1
1111111111111
1   1   1   1
1111111111111

More on Pip

Perl, 72 63 52 50 bytes
My shortest version yet uses $/ to get a newline char more compactly:
$ perl -e 'print((($,="="."   ="x pop.$/)=~s/./=/gr)x(1+pop))' 2 5
=====================
=   =   =   =   =   =
=====================
=   =   =   =   =   =
=====================

The previous update puts the mostly empty lines in the output record separator $,, and prints a list of continuous lines.
$ perl -e 'print((($,="="."   ="x pop."\n")=~s/./=/gr)x(1+pop))' 2 5

The previous version might be a bit clearer for the casual reader:
$ perl -E 'say($y=($x="="."   ="x pop)=~s/./=/gr);for(1..pop){say$x;say$y}' 2 5

The first version used @ARGV instead of pop:
$ perl -E 'say($y=($x="="."   ="x$ARGV[1])=~s/./=/gr);for(1..$ARGV[0]){say$x;say$y}' 2 5

Python 2, 58 57 Bytes
Fairly straightforward implementation.
def f(x,y):a="="*(4*y+1);print(a+"\n="+"   ="*y+"\n")*x+a

Check it out here.
Thanks to Sp3000 for saving one byte.

GNU sed -r, 160
Sigh, I thought this would be smaller, but here it is anyway.  Unfortunately sed regexes don't have any lookaround capability.
:
s/(.*)1$/   =\1/;t
s/([= ]+)/\1\n\1/
:b
s/   (.*\n)/===\1/;tb
s/(1*)1 $/\n\1/
:c
s/([^\n]*\n[^\n]*\n)(1*)1$/\1\1\2/;tc
s/(=+)(.*)/\1\2\1/
s/(^|\n)(.)/\1=\2/g

Taking input as unary from STDIN:
$ sed -rf boxes.sed <<< "11 111"
=============
=   =   =   =
=============
=   =   =   =
=============
$

CJam, 25 bytes
q~)S3*'=5*+4/f*3f>\)*1>N*

Try it online in the CJam interpreter.
How it works
q~                        e# Read H and W from STDIN.
   S3*'=5*+               e# Push "   =====".
           4/             e# Chop into ["   =" "===="].
  )          f*           e# Repeat each string W+1 times.
               3f>        e# Cut off the first three characters.
                  \)*     e# Repeat the resulting array H+1 times.
                     1>   e# Remove the first string.
                       N* e# Join the lines.

Marbelous, 168 bytes
This answer only works up to 255x255, not 1000x1000, due to limitations of the Marbelous language. I'm working on a 32-bit library, but it's not going to be ready any time soon.
Input is provided as two command line parameters or function parameters to the main board.
@2@3}1@0
SLEL//>0\/
@3@1}0--
}1&0@0/\&0
@1/\@2}1\/
:SL
..}0@0
&0/\>0&0
EN..--\/
{0@0/\ES
:EL
..@0
..>0EN
}0--\/
@0/\EX
:EX
}0
\/3D3D3D3D
:ES
}0
\/3D202020
:EN
}0
{03D0A

Pseudocode:
MB(H,W):
    EL(W)
    for 1..H:
        SL(W)
        EL(W)
EL(W):
    for 1..W:
        EX()
    EN()
SL(W):
    for 1..W:
        ES()
    EN()
EX():
    print "===="
ES():
    print "=   "
EN():
    print "=\n"

Python 3 2, 160 87 85 79 bytes
I know this can be golfed a lot more, I would like some suggestions as this is my first try at golfing.
def d(x,y):
    for i in[1]*x:print'='*(4*y+1)+'\n'+'=   '*(y+1)
    print'='*(4*y+1)

Thanks to @Cool Guy and @Sp3000's tips, I narrowed the size down to just above below half.
Eg: d(3,3)
=============
=   =   =   =
=============
=   =   =   =
=============
=   =   =   =
=============

Try it out here.
Excuse the trailing whitespaces.

Indentation level can be reduced.

You don't need to build a list then join. Strings can be multiplied '='*(4*y+1)

Setting w=4*y+1 saves 3 bytes.

@Cool Guy I'm using a tab, not 4 spaces.

Oh ok. Didn't notice that.

CJam 52 51 46 41 bytes
l~:B;:A;'=:U{{U4*}B*N}:V~{U{SZ*U}B*NUV}A*

Thanks to Sp3000 for -5 chars
Thanks to Martin Büttner♦ for another 5 chars
Try it

c function, 81
x,y;f(h,w){for(;y<=h*2;y++)for(x=0;x<w*4+2;x++)putchar(x>w*4?10:x&3&&y&1?32:61);}

Test program:
x,y;f(h,w){for(;y<=h*2;y++)for(x=0;x<w*4+2;x++)putchar(x>w*4?10:x&3&&y&1?32:61);}

int main (int argc, char **argv)
{
  f(2, 3);
}

I dropped a few characters to treat it as a single line instead of double for: x;f(h,w){for(w=w*4+2;x<=w*h*2+w;x++)putchar(x%w?x/w%2?x%w%4!=1?32:61:61:10);} -- 78

Should have looked at the other answers first =/, my comment is a longer version of Reto Koradi's answer.

yes, I tried quite hard (and failed) to get this into a single (shorter) loop

PHP4.1, 76 71 69 bytes
This is as golfed as I can get.
$R=str_repeat;echo$R(($T=$R('-',4*$H+1))."
|{$R('   |',$H)}
",$V),$T;

This expects the key H to have the number of lines and V the number of boxes per line.
I'm using - and | just so the boxes actually look like boxes. If required, I can change it to =. It doesn't matter the char that is used.
Leaving - and | also helps to understand the mess that's going on.

Old method, 76 bytes:
for($R=str_repeat;$H--;)echo$T=$R('-',($V*4)+1),"
|{$R('   |',$V)}
";echo$T;

Example of output:
http://localhost/file.php?H=2&V=3

-------------
|   |   |   |
-------------
|   |   |   |
-------------

Julia, 59 bytes
(n,m)->(l="="^(4m+1)*"\n";print(l*("="*"   ="^m*"\n"*l)^n))

This creates an unnamed function that accepts two integers as input and prints the result to STDOUT. To call it, give it a name, e.g. f=(n,m)->....
Ungolfed + explanation:
function f(n, m)
    # Store the solid row
    l = "="^(4m + 1) * "\n"

    # Print all rows by repeating a pattern n times
    print(l * ("=" * "   ="^m * "\n" * l)^n)
end

Examples:
julia> f(2, 3)
=============
=   =   =   =
=============
=   =   =   =
=============

julia> f(1, 5)
=====================
=   =   =   =   =   =
=====================

Any suggestions are welcome.

R, 83 81
As an unnamed function taking h and w as parameters.
Builds the 1st and second lines into a vector for each row and replicates it h times.  Appends a vector for the bottom line and cats out the vector using fill to restrict the length of the lines.  Now takes advantage of the any visible character rule.
function(h,w)cat(rep(c(A<-rep(1,w*4+2),rep('   1',w)),h),A[-1],fill=w*4+1,sep='')

Test run
> f=function(h,w)cat(rep(c(A<-rep(1,w*4+2),rep('   1',w)),h),A[-1],fill=w*4+1,sep='')
> f(4,2)
111111111
1   1   1
111111111
1   1   1
111111111
1   1   1
111111111
1   1   1
111111111
>

bash + coreutils, 57
dc -e"2do$2 4*1+^1-pd8*F/1+$1si[fli1-dsi0<c]dscx"|tr 0 \

This uses dc to print binary numbers that have 1s for the box edges and 0s for the spaces.

the all-ones number X is calculated by 2 ^ (width * 4 + 1) - 1, then pushed and printed
the 10001...0001 number is calculated by X* 8 / 15 + 1, then pushed
the stack is then dumped height times

The tr then converts the 0s to space characters.
Output
$ ./boxes.sh 2 4
11111111111111111
1   1   1   1   1
11111111111111111
1   1   1   1   1
11111111111111111
$

JavaScript, 92 85 bytes
I had hoped this would be shorter than the other JS answer (nice work as always, edc65), but oh well. I have a feeling the math here can be further golfed.
f=(r,c)=>(x=>{for(i=s='';(y=i++/x)<r-~r;)s+=i%x?' *'[-~y%2|!(-~i%4)]:'\n'})(4*c+2)||s

Sorry, can't help with the math - my head is spinning ... but here is some microop: f=(r,c)=>(x=>{for(i=s='';(y=i++/x)<r-~r;)s+=i%x?' *'[-~y%2|!(-~i%4)]:'\n'})(4*c+2)||s --> 85

@edc65 that's great, thanks! Things like 2*r+1 => r-~r are exactly what I meant by golfing the math, and that particular one is genius. :)

KDB(Q), 37 bytes
{(3+2*x-1)#(5+4*y-1)#'(4#"=";"=   ")}

Explanation
                      (4#"=";"=   ")     / box shape
           (5+4*y-1)#'                   / extend y axis
 (3+2*x-1)#                              / extend x axis
{                                   }    / lambda

Test
q){(3+2*x-1)#(5+4*y-1)#'(4#"=";"=   ")}[2;5]
"====================="
"=   =   =   =   =   ="
"====================="
"=   =   =   =   =   ="
"====================="

CJam, 30 29 bytes
New version with redundant + at the end removed (thanks, Dennis):
l~_4*)'=*N+:F\'=S3*+*'=+N++*F

I know that Dennis already posted a CJam solution that beats this by miles. But this is my very first attempt at CJam, so it's a miracle that it works at all. :)
Fairly brute force. Builds the first line from 4 * H + 1 = signs, then the second line from = repeated H times, with another = added. Then concatenates the two lines, repeats the whole thing V times, and then adds another copy of the first line.
It feels like I have way too many stack manipulations, and even ended up storing the first line in a variable because I had to shuffle stuff around on the stack even more otherwise.
Not very elegant overall, but you have to start somewhere... and I wanted to try a simple problem first.

You don't need the + at the end. CJam prints the entire stack.

Octave, 69 65 64
y=ones([2,4].*input()+1);y(1:2:end,:)=y(:,1:4:end)=4;char(y+31)

Thanks to DLosc for pointing out issues that led to -1
Takes input as [1 1] and outputs:

#####
#   #
#####

You can also just input '1' and get 1x1. If the input really needs to be 1 1, the size goes up to 88 85 84:
y=ones([2,4].*eval(['[',input(0,'s'),']'])+1);y(1:2:end,:)=y(:,1:4:end)=4;char(y+31)

Note: Matlab doesn't allow Octave's chaining or input(integer), but here is the Matlab version (67):
y=ones([2,4].*input('')+1);y(1:2:end,:)=4;y(:,1:4:end)=4;char(y+31)

C, 76 bytes
w,n;f(r,c){for(w=c*4+2,n=w*r*2+w;n--;)putchar(n%w?n/w%2&&n%w%4-1?32:61:10);}

Invoked as a function with number of rows and columns as arguments. For example:
f(5, 2)

CJam, 23
q~[F8]f{2b*1+' f+N}*_0=

Try it online
Explanation:
q~        read and evaluate the input (height width)
[F8]      make an array [15 8] - 1111 and 1000 in base 2
f{…}      for width and each of [15 8]
  2b      convert the number to base 2
  *       repeat the digits "width" times
  1+      append a 1 to the array of digits (right edge)
  ' f+    add the space character to each digit (0->' ', 1->'!')
  N       push a newline
*         repeat the resulting array "height" times
_0=       copy the first line (bottom edge)

Brainfuck, 157 156 152 (123 120 118 117 with unary)
This works with Alex Pankratov's bff (brainfuck interpreter used on SPOJ and ideone) and Thomas Cort's BFI (used on Anarchy Golf).
The input must contain a single pair of integers (in decimal) separated by a space, without a trailing newline.
Demonstration on ideone.
,+
[
  >++++++[<-[------->>] >[->>] <<-]
  <
  [
    -
    <[>++++++++++<-]
    <
  ]
  >>[<+>-]
  ,+
]
+++++++++++[>+++>+++>+<<<-]
<<<+
[
  >>
  [<+> >>>.<...<<-]
  >+>+>.>-.+
  [<]
  <
  [
    <->
    >>>--<,<
  ]
  <[>+<-]
  <
]

Brief explanation:
At the end of first loop we are left with
r 0 c [0]
We then change this to
[r+1] 0 c 0 33 33 11 0
and begin the main loop. In that loop, the tricky part is alternating between 33 (exclamation mark) and 32 (space), and alternating between decrementing and not decrementing the number of rows. This is done with a flag using the cell to the right of c.
Update:
I had not considered unary input until reading some other answers. Here is a unary version:
+
[
  ->++++[<-------->-]
  <[<+<]
  >>,+
]
+++++++++++[>+++>+++>+<<<-]
<<<
[
  >>
  [<+> >>>.<...<<-]
  >+>+>.>-.+
  [<]
  <
  [
    <->
    >>>--<,<
  ]
  <[>+<-]
  <
]

Lua, 95 93 85 bytes.
Works with command-line arguments.
Gotta love those string methods
a=arg[2]c="="l=c:rep(a*2+1)print((l.."\n"..(("= ")):rep(a)..c.."\n"):rep(arg[1])..l)

You can replace c..(" "):rep(3) by ("=   ")

Python 2.7, 66 bytes.
I know there are already better solutions in python, but that's the best I came up with.
def f(x,y):a,b,c="="*4,"   =","\n=";print "="+a*y+(c+b*y+c+a*y)*x

Example:
f(3,4)

=================
=   =   =   =   =
=================
=   =   =   =   =
=================
=   =   =   =   =
=================

Sorry, but this produces wrong output for test cases 2,3 and 5. You mixed up the columns and rows and got it the opposite way.

SAS, 117 119
macro a array a[0:1]$4('#   ' '####');do x=1 to 2+2*&x-1;t=repeat(a[mod(x,2)],&y-1);put t $char%%eval(&y*3). '#';end;%

Example:
%let x=4;
%let y=4;
data;a;run;

Result:
#############
#   #   #   #
#############
#   #   #   #
#############
#   #   #   #
#############
#   #   #   #
#############

Is there any online compiler where I could test this? BTW, as per your result, your program produces wrong output. See output2 : Each box should have three horizontal spaces in them

@CoolGuy, you're right, I did not catch that, updated my post.  You can try SAS On Demand, would be the least troublesome way to access an online compiler without setting up your own AWS Instance

Your first link doesn't work for me. Chrome gives DNS_PROBE_FINISHED_NXDOMAIN

Try this one, although the first link works fine for me SAS On Demand or follow the link on this page to the 'Control Center' here

I don't know why, but the going to the control center causes the same error as mentioned in my earlier comment :/

Java, 123 119 bytes
void p(int h,int w){String b="=",d="\n=";for(;w-->0;d+="   =")b+="====";for(d+="\n"+b;h-->0;b+=d);System.out.print(b);}

Abusing the input parameters as counters greatly helped in decreasing code size.
Thanks to Cool Guy for suggesting the abuse of for syntax.

Golf it more by using a for loop instead of a while loop.

Unfortunately for(;w-->0;) is the same length as while(w-->0)

I meant while(w-->0){d+="   =";b+="====";} --> for(;w-->0;b+="====")d+="   ="; which saves three bytes

You are right. I actually managed to squeeze 4 bytes out of this thanks to your suggestion.

GolfScript, 48 bytes
2 5:b;:a;0`:z;{4b*1+z*n+}:r~{{z' '3*}b*z+n+r}a*

Takes the size as the first two inputs (here the first two digits).
Demonstration
000000000000000000000
0   0   0   0   0   0
000000000000000000000
0   0   0   0   0   0
000000000000000000000

Java, 181
I hope that according to

You are permitted to write a full program or a function.

it is compliant to the rules to count the bytes of the function, which is 181 in this case
import static java.lang.System.*;
public class Boxes
{
    public static void main(String[] args)
    {
        Boxes b=new Boxes();
        System.out.println("1,1:");
        b.b(1,1);
        System.out.println("1,2:");
        b.b(1,2);
        System.out.println("2,1:");
        b.b(2,1);
        System.out.println("2,2:");
        b.b(2,2);
        System.out.println("2,5:");
        b.b(2,5);
    }

    void b(int R, int C){String s="",e=s,x,y,z=s,a="====",n="=\n";int r,c;for(r=R;r-->0;){x=y=e;for(c=C;c-->0;){x+=a;y+="=   ";}s+=x+n+y+n;}for(c=C;c-->0;){z+=a;}s+=z+n;out.println(s);}
}

The output is wrong. See output no 2 : "Each box should have three horizontal spaces in them". Your code outputs two spaces, not three

@CoolGuy Miscounted this - now 2 bytes more, but that doesn't change much...

ok. +1. Save two bytes by changing for(r=R;r-->0;){x=y=e;for(c=C;c-->0;){x+=a;y+="=   ";}s+=x+n+y+n;} to for(r=R;r-->0;s+=x+n+y+n){x=y=e;for(c=C;c-->0;y+="=   "){x+=a;}}

C#, 153 151 150
This can't really compete but here it is just for fun:
(h,w)=>{string s="=",e="\n",a="====",b="   =",m=a,o;int i=0,j;for(;i++<=h*2;){o=s;for(j=0;j++<w+1;o=m)System.Console.Write(o);m=m==a?b:a;s=e+s;e="";}}

How to run:
public class Program
{
    public static void Main()
    {
        new System.Action<int, int>((h,w)=>{string s="=",e="\n",a="====",b="   =",m=a,o;int i=0,j;for(;i++<=h*2;){o=s;for(j=0;j++<w+1;o=m)System.Console.Write(o);m=m==a?b:a;s=e+s;e="";}})(3, 4);
    }
}

Open for improvements.

replace string with var.

Unfortunately that's not allowed. See this

Declaring those ints outside the loop can save a byte.

Dart, 57
b(y,x,{z:'='}){z+=z*4*x;print('$z\n=${'   ='*x}\n'*y+z);}

Try it at: https://dartpad.dartlang.org/36ed632613395303ef51

Ruby, 57 48 45
f=->h,w{l=?=*w*4+?=;(l+$/+'=   '*w+"=
")*h+l}

Usage:
print f[2,5]

Thanks to manatwork for saving 3 bytes.

Two small improvement possibilities: '=' → ?= and "<LF>" → $/.

Another small one: ?=*(w*4+1) → ?=+?=*w*4

Haskell, 153
At least it's beating Java!
s r c
 |r==1=a c
 |otherwise=s (r-1)c++(take 2$tail(a c))
 where
 a 1=[f,"=   =",f]
 a n=map(\x->x++((take 4.tail)x))(a(n-1))
 f="====="

p r=unlines.s r

Works in ghci, although ghci doesn't take a new line on \n.
Here is the ungolfed version:
import System.Environment

xsquare :: Int -> [String]
xsquare 1 = ["=====","=   =","====="]
xsquare n = map (\x -> x++((take 4 . tail) x)) (xsquare (n-1))

ysquare :: Int -> [String] -> [String]
ysquare 1 x = x
ysquare n x = ysquare (n-1) (x ++ (take 2 $ tail x))

squares :: Int -> Int -> [String]
squares x y = ysquare y (xsquare x)

main = do
  a <- getArgs
  let [x,y] = (map read . take 2) a
  putStrLn $ unlines $ squares y x

It calculates the horizontal shape first, then finds the vertical shape by adding the last two lines of the horizontal shape the appropriate number of times.
The main function is so I could test it properly with command line arguments and proper newline behaviour.