/*
(c) Copyright 2001 by Sergei Bespamyatnikh, All Rights Reserved 

The algorithm computes vertices of edgewise subdivision [1] of tetrahedron
with vertices (1,0,0,0), (0,1,0,0), (0,0,1,0), (0,0,0,1) in R^4.

[1] H. Edelsbrunner and D. R. Grayson, 
Edgewise subdivision of a simplex. DCG, 24:707--719, 2000.

*/

#include <stdio.h>      /* standard IO: FILE, NULL (*), EOF, ... */
#include <stdlib.h>     /* standard lib: abort(), system(), getenv(), ... */
#include <ctype.h>      /* char macros: isalpha(), ... */
#include <math.h>       /* math lib: sin(), ...; use with -lm */
#include <string.h>     

#define PI 3.141593

long int NV=0;

typedef struct quat
{
  float o, x, y, z;
} quat;

quat V[200000];

/*********************************************************************/
float dist (quat *a, quat *b)
{
  float o, x, y, z;
  o = a->o - b->o;
  x = a->x - b->x;
  y = a->y - b->y;
  z = a->z - b->z;
  return (o*o+x*x+y*y+z*z);
}

/*********************************************************************/
int add_vertex (quat *q)
{
  int i;
  float di;
  
  for (i=0; i<NV; i++)
  {
    di = dist (V+i, q);
    
    if (di<0.0000001)
    {
      return i;
    } 
  }
  printf ("\n %f %f %f %f", q->o, q->x, q->y, q->z);
  V[NV] = *q;
  NV++;
  return NV-1;
}


/*********************************************************************/
void subd2 (int *a, int l)
{
  int s, i, j, k, M[4][100], b[4], c[4], x;
  quat q, r, q0 = {0., 0., 0., 0.};
  int iq, jq;
  
  k = l;
  for (i=0; i<k; i++)
  {
    for (j=0; j<3; j++)
    {
      M[j][i] = 0;
    }
  }
  b[3] = -1;
 
  for (b[0]=0; b[0]<k; b[0]++)
  {
    for (b[1]=0; b[1]<k; b[1]++)
    {
      for (b[2]=0; b[2]<k; b[2]++)
      {
        s = 0;
        for (j=0; j<k; j++)
        {
          for (i=0; i<=3; i++)
          {
            M[i][j] = s;

            if (b[i]==j) 
            { 
              s++;
            }
            
          }
        }
        
        for (i=0; i<=3; i++)
        {
          q = q0;
          for (j=0; j<k; j++)
          {
            if (M[i][j]<0 || M[i][j]>3) 
            {
              printf("\nerror M[%d][%d]=%d\n",i,j,M[i][j]);
              printf("k=%d l=%d b=%d %d %d %d",k,l,b[0],b[1],b[2],b[3]);
              for (jq=0; jq<k; jq++)
              {
                printf("\n");
                for (iq=0; iq<=3; iq++)
                {
                  printf("%d",M[iq][jq]);
                }
              }
              exit(0);
            }
            if (a[M[i][j]]<0 || a[M[i][j]]>=NV)
            {
              printf("\nerror a %d %d\n", a[M[i][j]], M[i][j]);
              exit(0);
            }
            r = V[a[M[i][j]]];
            q.o += r.o;
            q.x += r.x;
            q.y += r.y;
            q.z += r.z;
          }
          q.o /= k;
          q.x /= k;
          q.y /= k;
          q.z /= k;
          x = NV;
          c[i] = add_vertex (&q);
        }
      }
    }
  }
}


/*********************************************************************/
void subd1 (int a0, int a1, int a2, int a3, int l)
{
  int a[4];
  
  a[0] = a0;
  a[1] = a1;
  a[2] = a2;
  a[3] = a3;
  
  subd2 (a, l);
}


/*********************************************************************/
void subd (int l)
{
  int i;
  quat q = {0.,0.,0.,0.};
  
  NV = 4;
  
  for (i=0; i<8; i++)
  {
    V[i] = q;
  }
  
  V[0].o = 1.;
  V[1].x = 1.;
  V[2].y = 1.;
  V[3].z = 1.;
  printf ("\n %f %f %f %f", V[0].o, V[0].x, V[0].y, V[0].z);
  printf ("\n %f %f %f %f", V[1].o, V[1].x, V[1].y, V[1].z);
  printf ("\n %f %f %f %f", V[2].o, V[2].x, V[2].y, V[2].z);
  printf ("\n %f %f %f %f", V[3].o, V[3].x, V[3].y, V[3].z);
  subd1 (0, 1, 2, 3, l);
}


/*********************************************************************/
int 
main ( int argc, char *argv[])
{
  int nl, i;
  

  if (argc != 2)
  {
  printf("\nSyntax: subd n");
  exit(0);
  }

  sscanf (argv[1], "%d", &nl);

  subd (nl);
  
}