#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

const unsigned DIM = 3;

vector<double> sa(const char point,
		  const vector<double> &a,
		  const vector<double> &b,
		  const vector<double> &c,
		  const vector<double> &d){ 
  vector<double> p(DIM, 0.0);

  switch(point){
  case 'a':
    return(a);
  case 'b':
    p[0] = a[0] + 0.5*(b[0]-a[0]);
    p[1] = a[1] + 0.5*(b[1]-a[1]);
    p[2] = a[2] + 0.5*(b[2]-a[2]);
    break;
  case 'c':
    p[0] = a[0] + 0.5*(c[0]-a[0]);
    p[1] = a[1] + 0.5*(c[1]-a[1]);
    p[2] = a[2] + 0.5*(c[2]-a[2]);
    break;
  case 'd':
    p[0] = a[0] + 0.5*(d[0]-a[0]);
    p[1] = a[1] + 0.5*(d[1]-a[1]);
    p[2] = a[2] + 0.5*(d[2]-a[2]);
    break;
  }

  return(p);
}

vector<double> sb(const char point,
		  const vector<double> &a,
		  const vector<double> &b,
		  const vector<double> &c,
		  const vector<double> &d){
  vector<double> p(DIM, 0.0);

  switch(point){
  case 'a':
    p[0] = b[0] + 0.5*(a[0]-b[0]);
    p[1] = b[1] + 0.5*(a[1]-b[1]);
    p[2] = b[2] + 0.5*(a[2]-b[2]);
    break;
  case 'b':
    return(b);
  case 'c':
    p[0] = b[0] + 0.5*(c[0]-b[0]);
    p[1] = b[1] + 0.5*(c[1]-b[1]);
    p[2] = b[2] + 0.5*(c[2]-b[2]);
    break;
  case 'd':
    p[0] = b[0] + 0.5*(d[0]-b[0]);
    p[1] = b[1] + 0.5*(d[1]-b[1]);
    p[2] = b[2] + 0.5*(d[2]-b[2]);
    break;
  }

  return(p);
}

vector<double> sc(const char point,
		  const vector<double> &a,
		  const vector<double> &b,
		  const vector<double> &c,
		  const vector<double> &d){
  vector<double> p(DIM, 0.0);

  switch(point){
  case 'a':
    p[0] = c[0] + 0.5*(a[0]-c[0]);
    p[1] = c[1] + 0.5*(a[1]-c[1]);
    p[2] = c[2] + 0.5*(a[2]-c[2]);
    break;
  case 'b':
    p[0] = c[0] + 0.5*(b[0]-c[0]);
    p[1] = c[1] + 0.5*(b[1]-c[1]);
    p[2] = c[2] + 0.5*(b[2]-c[2]);
    break;
  case 'c':
    return(c);
  case 'd':
    p[0] = c[0] + 0.5*(d[0]-c[0]);
    p[1] = c[1] + 0.5*(d[1]-c[1]);
    p[2] = c[2] + 0.5*(d[2]-c[2]);
    break;
  }
		
  return(p);
}

vector<double> sd(const char point,
		  const vector<double> &a,
		  const vector<double> &b,
		  const vector<double> &c,
		  const vector<double> &d){
  vector<double> p(DIM, 0.0);

  switch(point){
  case 'a':
    p[0] = d[0] + 0.5*(a[0]-d[0]);
    p[1] = d[1] + 0.5*(a[1]-d[1]);
    p[2] = d[2] + 0.5*(a[2]-d[2]);
    break;
  case 'b':
    p[0] = d[0] + 0.5*(b[0]-d[0]);
    p[1] = d[1] + 0.5*(b[1]-d[1]);
    p[2] = d[2] + 0.5*(b[2]-d[2]);
    break;
  case 'c':
    p[0] = d[0] + 0.5*(c[0]-d[0]);
    p[1] = d[1] + 0.5*(c[1]-d[1]);
    p[2] = d[2] + 0.5*(c[2]-d[2]);
    break;
  case 'd':
    return(d);
  }
		
  return(p);
}

// Sierpinski gasket function computes midpoints of each side,
// saves the midpoints in a vector, then recursively calls itself
// on the resulting three inner tetrahedra.
void sg(const int depth,
	const vector<double> &a,
	const vector<double> &b,
	const vector<double> &c,
	const vector<double> &d,
	vector< vector<double> > &points){
  
  if(depth == 0) return;
  
  vector<double> ab, ac, ad, bc, bd, cd;
  ab = sa('b',a,b,c,d);
  ac = sa('c',a,b,c,d);
  ad = sa('d',a,b,c,d);
  bc = sb('c',a,b,c,d);
  bd = sb('d',a,b,c,d);
  cd = sc('d',a,b,c,d);

  points.push_back(ab);
  points.push_back(ac);
  points.push_back(ad);
  points.push_back(bc);
  points.push_back(bd);
  points.push_back(cd);
  
  sg(depth-1,  a, ab, ac, ad, points);
  sg(depth-1, ab,  b, bc, bd, points);
  sg(depth-1, ac, bc,  c, cd, points);
  sg(depth-1, ad, bd, cd,  d, points);
}

int main(){
  
  vector<double> a(DIM), b(DIM), c(DIM), d(DIM);

  // Initialize the four points of the tetrahedron
  a[0]=0.5; a[1]=sqrt(3.0)/4.0; a[2]=3.0/4.0;
  b[0]=0.0; b[1]=0.0; b[2]=0.0;
  c[0]=0.5; c[1]=sqrt(3.0)/2.0; c[2]=0.0;
  d[0]=1.0; d[1]=0.0; d[2]=0.0;
  
  // An alternative initial tetrahedron
  //  a[0]=-1.0; a[1]=+1.0; a[2]=-1.0;
  //  b[0]=+1.0; b[1]=+1.0; b[2]=+1.0;
  //  c[0]=-1.0; c[1]=-1.0; c[2]=+1.0;
  //  d[0]=+1.0; d[1]=-1.0; d[2]=-1.0;

  int depth = 0;
  cerr << "Enter depth (integer > 0): ";
  cin >> depth;

  // This vector holds all data points generated
  // via recursive calls to sg()
  vector< vector<double> > points;

  points.push_back(a);
  points.push_back(b);
  points.push_back(c);
  points.push_back(d);

  // Recursively generate Sierpinski's gasket
  sg(depth-1, a, b, c, d, points);

  for(vector< vector<double> >::size_type i=0; i < points.size(); ++i){
    cout << points[i][0] << ' ' 
	 << points[i][1] << ' '
	 << points[i][2] << endl;
  }
  
  return 0;
}