Task
Produce a graphical or ASCII-art representation of a Sierpinski carpet of order   N.



For example, the Sierpinski carpet of order   3   should look like this:

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

The use of the   #   character is not rigidly required for ASCII art.

The important requirement is the placement of whitespace and non-whitespace characters.



Related task

   Sierpinski triangle


 
\\ Improved simple plotting using matrix mat (color and scaling added). 
\\ Matrix should be filled with 0/1. 7/6/16 aev
iPlotmat(mat,clr)={
  my(xz=#mat[1,],yz=#mat[,1],vx=List(),vy=vx,xmin,xmax,ymin,ymax,c=0.625);
  for(i=1,yz, for(j=1,xz, if(mat[i,j]==0, next, listput(vx,i); listput(vy,j))));
  xmin=listmin(vx); xmax=listmax(vx); ymin=listmin(vy); ymax=listmax(vy);
  plotinit(0); plotcolor(0,clr);
  plotscale(0, xmin,xmax,ymin,ymax);
  plotpoints(0, Vec(vx)*c,Vec(vy));
  plotdraw([0,xmin,ymin]);
  print(" *** matrix: ",xz,"x",yz,", ",#vy," DOTS");
}
\\ iPlotV2(): Improved plotting from a file written by the wrtmat(). (color added)
\\ Saving possibly huge generation time if re-plotting needed.
iPlotV2(fn, clr)={
  my(F,nf,vx=List(),vy=vx,Vr,xmin,xmax,ymin,ymax,c=0.625);
  F=readstr(fn); nf=#F;
  print(" *** Plotting from: ", fn, " - ", nf, " DOTS");
  for(i=1,nf, Vr=stok(F[i]," "); listput(vx,eval(Vr[1])); listput(vy,eval(Vr[2])));
  xmin=listmin(vx); xmax=listmax(vx); ymin=listmin(vy); ymax=listmax(vy);
  plotinit(0); plotcolor(0,clr);
  plotscale(0, xmin,xmax,ymin,ymax);  
  plotpoints(0, Vec(vx)*c,Vec(vy));
  plotdraw([0,xmin,ymin]);
}
\\ Are x,y inside Sierpinski carpet? (1-yes, 0-no) 6/10/16 aev  
inSC(x,y)={
  while(1, if(!x||!y,return(1)); 
    if(x%3==1&&y%3==1, return(0));
    x\=3; y\=3;);\\wend
}
 