#include <stdio.h>
#include <conio.h>
#include <dos.h>
#include <mem.h>
#include <stdlib.h>
#include <math.h>
#include <values.h>

#define MAX_DIM  7    /* Don't forget to change MAX_EDGE accordingly */
#define MAX_EDGE 448   /* Depends on MAX_DIM */
#define X_MAG   60.0
#define Y_MAG   50.0
#define X_OFFS  160.0
#define Y_OFFS  100.0
#define X_RES   320
#define Y_RES   200
#define VIDEO 0xA000
#define put_pixel(x,y,color) canvas[x+(y<<8)+(y<<6)] = color
#define get_pixel(x,y) canvas[x+(y<<8)+(y<<6)]
#define RAD(x) (float)(x)*M_PI/180

#define VERT_SYNC1  asm push dx;     \
		    asm push ax;     \
		    asm mov  dx,3dah;\
		VRT1:                \
		    asm in   al,dx;  \
		    asm test al,8;   \
		    asm jnz  VRT1;   \
		NoVRT1:              \
		    asm in   al,dx;  \
		    asm test al,8;   \
		    asm jz   NoVRT1; \
		    asm pop  ax;     \
		    asm pop  dx;     \

#define VERT_SYNC2  asm push dx;     \
		    asm push ax;     \
		    asm mov  dx,3dah;\
		VRT2:                \
		    asm in   al,dx;  \
		    asm test al,8;   \
		    asm jnz  VRT2;   \
		NoVRT2:              \
		    asm in   al,dx;  \
		    asm test al,8;   \
		    asm jz   NoVRT2; \
		    asm pop  ax;     \
		    asm pop  dx;     \

#define VERT_SYNC3  asm push dx;     \
		    asm push ax;     \
		    asm mov  dx,3dah;\
		VRT3:                \
		    asm in   al,dx;  \
		    asm test al,8;   \
		    asm jnz  VRT3;   \
		NoVRT3:              \
		    asm in   al,dx;  \
		    asm test al,8;   \
		    asm jz   NoVRT3; \
		    asm pop  ax;     \
		    asm pop  dx;     \

unsigned char far *SCREEN=MK_FP(VIDEO,0);
unsigned char far *canvas;

struct {
    float vert1[MAX_DIM];
    float vert2[MAX_DIM];
  } ncube[MAX_EDGE];
int prev_cube_draw[MAX_EDGE][4];
float rot[MAX_DIM][MAX_DIM];
float rotm[MAX_DIM][MAX_DIM];
int n;   /* dimension of vector space */
int e;
int cube_size;
int mat_size;
int colpal; //color palettes, user index into palcol
struct {
  float a;
  float b;
  } palcol[6]; //subpalette half-widths and centers
int defvgapal[256][3]; //to save the default vga palette
float ninc_targ_bridge_length; //incrementing n destination bridge length
float n_inc_ref[MAX_DIM]; //incrementing n reference vector for bridge expansion animation
int vertmchs[2*MAX_EDGE], vertmatch[2*MAX_EDGE]; //vertmatch search results, see vertmatcher()
float pushes[2*MAX_EDGE][MAX_DIM]; //2e is the max number of different verts in ncube,
				   //might make do with less here,
			  //,or rely on the smallness of the animation step and do without a push array,
			  //pushing the verts while doing the pushes' calc loop
  //, or, if using float pushes[2*MAX_EDGE][MAX_DIM], why not just copy the whole ncube instead?
int pushvi[2*MAX_EDGE]; //1 or -1 flag for whether there is or isn't an active vector in push[2e+0or1]

void vga_mode(int mode)
{
  union REGS regs;

  regs.x.ax = mode;
  int86(0x10, &regs, &regs);
}

void clear_screen(void)
{
  _fmemset(canvas, 0, 64000);
}

void update_screen(void)
{
}

void old_draw_line(int x1, int y1, int x2, int y2, int color)
{
  register int dx, dy, incr1, incr2, d, x, y, xend, yend, yinc, xinc;

  dx = abs(x2-x1); dy=abs(y2-y1);
  if (dx >= dy)
  {
    if (x1 >x2)
    {
      x=x2;y=y2;xend=x1;
      if(dy==0)	yinc = 0;
      else
      {
	if (y2>y1) yinc = -1;
	else yinc = 1;
      }
    }
    else
    {
      x = x1; y = y1; xend = x2;
      if (dy==0) yinc=0;
      else
      {
	if (y2>y1) yinc = 1;
	else yinc = -1;
      }
    }

    incr1 = 2*dy; d = incr1-dx; incr2 = 2*(dy-dx);
    if (x >= 0 && x < X_RES && y >= 0 && y < Y_RES)
      put_pixel(x,y,color);
    while (x<xend)
    {
      x++;
      if (d<0) d+=incr1;
      else
      {
	y += yinc;
	d = d+incr2;
      }
      if (x >= 0 && x < X_RES && y >= 0 && y < Y_RES)
	put_pixel(x,y,color);
    }
  }
  else
  {
    if (y1 >y2)
    {
      x=x2;y=y2;yend=y1;
      if(dx==0)	xinc = 0;
      else
      {
	if (x2>x1) xinc = -1;
	else xinc = 1;
      }
    }
    else
    {
      x = x1; y = y1; yend = y2;
      if (dx==0) xinc=0;
      else
      {
	if (x2>x1) xinc = 1;
	else xinc = -1;
      }
    }

    incr1 = 2*dx; d = incr1-dy; incr2 = 2*(dx-dy);
    if (x >= 0 && x < X_RES && y >= 0 && y < Y_RES)
      put_pixel(x,y,color);
    while (y<yend)
    {
      y++;
      if (d<0) d+=incr1;
      else
      {
	x += xinc;
	d = d+incr2;
      }
      if (x >= 0 && x < X_RES && y >= 0 && y < Y_RES)
	put_pixel(x,y,color);
    }
  }
}

void draw_line(int x1, int y1, int x2, int y2, int col1, int col2)
{
  register int dx, dy, incr1, incr2, d, x, y, xend, yend, yinc, xinc;
  int color, cs; float colf;

  dx = abs(x2-x1); dy=abs(y2-y1);
  if (dx >= dy)
  {
    if (x1 >x2)
    {
      x=x2;y=y2;xend=x1; color=col2; cs = col1-col2;
      if(dy==0)	yinc = 0;
      else
      {
	if (y2>y1) yinc = -1;
	else yinc = 1;
      }
    }
    else
    {
      x = x1; y = y1; xend = x2; color=col1; cs = col2-col1;
      if (dy==0) yinc=0;
      else
      {
	if (y2>y1) yinc = 1;
	else yinc = -1;
      }
    }

    incr1 = 2*dy; d = incr1-dx; incr2 = 2*(dy-dx);
    if (x >= 0 && x < X_RES && y >= 0 && y < Y_RES)
      put_pixel(x,y,color);
    colf=(float)color;
    while (x<xend)
    {
      x++;
      colf = colf + (float)cs/(float)dx;
      color = (int)colf;
      if (d<0) d+=incr1;
      else
      {
	y += yinc;
	d = d+incr2;
      }
      if (x >= 0 && x < X_RES && y >= 0 && y < Y_RES)
	{
	put_pixel(x,y,color);
	}
    }
  }
  else
  {
    if (y1 >y2)
    {
      x=x2;y=y2;yend=y1; color=col2; cs = col1-col2;
      if(dx==0)	xinc = 0;
      else
      {
	if (x2>x1) xinc = -1;
	else xinc = 1;
      }
    }
    else
    {
      x = x1; y = y1; yend = y2; color=col1; cs = col2-col1;
      if (dx==0) xinc=0;
      else
      {
	if (x2>x1) xinc = 1;
	else xinc = -1;
      }
    }

    incr1 = 2*dx; d = incr1-dy; incr2 = 2*(dx-dy);
    if (x >= 0 && x < X_RES && y >= 0 && y < Y_RES)
      put_pixel(x,y,color);
    colf=(float)color;
    while (y<yend)
    {
      y++;
      colf = colf + (float)cs/(float)dy;
      color = (int)colf;
      if (d<0) d+=incr1;
      else
      {
	x += xinc;
	d = d+incr2;
      }
      if (x >= 0 && x < X_RES && y >= 0 && y < Y_RES)
	put_pixel(x,y,color);
    }
  }
}

void draw_edge(int edg, int col)
{
  int x1, y1, x2, y2;
  float ompr1=0, omcs1=0, ompr2=0, omcs2=0, sqn;
  int col1=col, col2=col, i;

  x1 = X_MAG * (ncube[edg].vert1[0]) + X_OFFS;
  y1 = Y_MAG * (ncube[edg].vert1[1]) + Y_OFFS;
  x2 = X_MAG * (ncube[edg].vert2[0]) + X_OFFS;
  y2 = Y_MAG * (ncube[edg].vert2[1]) + Y_OFFS;

// calculate the pythagorean radical for the truncated component and
// shade pixels in the edge line in a gradient,
// the sense of it depending on the sum of the omitted coordinates

// omitted coordinates pythagorean radicals will range from 0 to root n
// omitted coordinates sums
//   negative omcs shaded bright, positive omcs shaded dimmer

  if (n>2 && col!=0)
  {
  for (i=2;i<n;i++)
  {
    omcs1 = omcs1 + ncube[edg].vert1[i];
    omcs2 = omcs2 + ncube[edg].vert2[i];
    ompr1 = ompr1 + (ncube[edg].vert1[i])*(ncube[edg].vert1[i]);
    ompr2 = ompr2 + (ncube[edg].vert2[i])*(ncube[edg].vert2[i]);
  }
// later maybe try ompr's as abs sums--no squares, no sqrts--xor no sqrts
  ompr1 = sqrt(ompr1);
  ompr2 = sqrt(ompr2);
  sqn = sqrt(n);
  if (omcs1<0) ompr1=(-1)*ompr1;
  if (omcs2<0) ompr2=(-1)*ompr2;
// colors will be 16 grays from color# 16 to 31 (black to white),
// until improved with custom palette
  col1 = (ompr1*palcol[colpal].a/sqn)+palcol[colpal].b;
  col2 = (ompr2*palcol[colpal].a/sqn)+palcol[colpal].b;
// look into int roundings here and in draw_line
  }
  draw_line(x1, y1, x2, y2, col1, col2);
  prev_cube_draw[edg][0]=x1;
  prev_cube_draw[edg][1]=y1;
  prev_cube_draw[edg][2]=x2;
  prev_cube_draw[edg][3]=y2;
}

void wipe_prev_cube(void)
{
  int i;
  int x1, y1, x2, y2, col=0;
  for (i=0;i<e;i++)
  {
    x1=prev_cube_draw[i][0];
    y1=prev_cube_draw[i][1];
    x2=prev_cube_draw[i][2];
    y2=prev_cube_draw[i][3];
    draw_line(x1, y1, x2, y2, col, col);
  }
}

void old_draw_edge(int edg, int col)
{
  float x1, y1, x2, y2;

  x1 = X_MAG * (ncube[edg].vert1[0]) + X_OFFS;
  y1 = Y_MAG * (ncube[edg].vert1[1]) + Y_OFFS;
  x2 = X_MAG * (ncube[edg].vert2[0]) + X_OFFS;
  y2 = Y_MAG * (ncube[edg].vert2[1]) + Y_OFFS;
  old_draw_line(x1, y1, x2, y2, col);
// new erase
  prev_cube_draw[edg][0]=x1;
  prev_cube_draw[edg][1]=y1;
  prev_cube_draw[edg][2]=x2;
  prev_cube_draw[edg][3]=y2;
}

void draw_ncube(int col)
{
  int i;

  for (i=0;i<e;i++)
    draw_edge(i, col);
}

void old_draw_ncube(int col)
{
  int i;

  for (i=0;i<e;i++)
    old_draw_edge(i, col);
}

void add_edge(float *v1,float *v2)
{
  int i, j, found;

  /* check if edge has already been added */
  found = 0;
  i = 0;
  while (!found && i<cube_size)
  {
    found = 1;
    for (j=0;j<n;j++)
      if (ncube[i].vert1[j] != v1[j] ||
	  ncube[i].vert2[j] != v2[j])
	found = 0;
    i++;
  }
  if (found)
    return;
  i=cube_size;
  for (j=0;j<n;j++)
  {
    ncube[i].vert1[j] = v1[j];
    ncube[i].vert2[j] = v2[j];
  }
  cube_size++;
}

void recurse_ncube(float *v)
{
  int i, j;
  float new[MAX_DIM];

  for (i=0;i<n;i++)
    if (v[i] == -1)
      {
	for (j=0;j<n;j++)
	  new[j] = v[j];
	new[i] = 1;
	add_edge(v, new);
	recurse_ncube(new);
      }
}

int edge_count()
{
  int i, edges=1;
  for(i=1; i<n; i++)
    edges=2*edges+pow(2,i);
  return edges;
}

void init_ncube()
{
  int i;
  float zero[MAX_DIM];

  e = edge_count();
  cube_size = 0;
  for (i=0;i<n;i++)
    zero[i] = -1;
  recurse_ncube(zero);
  if (cube_size != e)
    {
      vga_mode(0x03);
      printf("Trouble creating ncube!\n");
      exit(1);
    }
}

void set_rot_mat(float a, int b1, int b2)
{
  int i;

  memset(rot,0,mat_size);
  for (i=0;i<n;i++)
    rot[i][i] = 1;
  rot[b1][b1] = cos(a);
  rot[b1][b2] = sin(a);
  rot[b2][b1] = -sin(a);
  rot[b2][b2] = cos(a);
}

void rotate_point(float *p)
{
  int i, j;
  float result, newp[MAX_DIM];

  for (i=0;i<n;i++)	/* row of m */
  {
    result = 0;
    for (j=0;j<n;j++)   /* col of m and row of p */
      result += rot[i][j] * p[j];
    newp[i] = result;
  }
  for (j=0;j<n;j++)
    p[j] = newp[j];
}

void rotate_edge(int i2)
{
  rotate_point(ncube[i2].vert1);
  rotate_point(ncube[i2].vert2);
}

void rotate(float a, int b1, int b2)
{
  int i;

  set_rot_mat(a, b1, b2);
  for (i=0;i<e;i++)
    rotate_edge(i);
  //also rotate the reference vector for n_inc
  rotate_point(n_inc_ref);
}

void init_dimension(int d)
{
  n=d;
  mat_size=MAX_DIM*MAX_DIM*sizeof(float);
  init_ncube();
}

void pauser()
{
       while (!kbhit());
       getch();
}

void breakdist()
{
  int i,j;

  for (i=0;i<e;i++)
  {
    for (j=0;j<n;j++)
    {
    ncube[i].vert1[j]+=(float)rand()/(RAND_MAX*5.0)-0.1;
    ncube[i].vert2[j]+=(float)rand()/(RAND_MAX*5.0)-0.1;
    }
  }
}

void shrink()
{
  int i,j;

  for (i=0;i<e;i++)
  {
    for (j=0;j<n;j++)
    {
    ncube[i].vert1[j]*=.92;
    ncube[i].vert2[j]*=.92;
    }
  }
}

void grow()
{
  int i,j;

  for (i=0;i<e;i++)
  {
    for (j=0;j<n;j++)
    {
    ncube[i].vert1[j]*=1.087;
    ncube[i].vert2[j]*=1.087;
    }
  }
}

void distort()
{
// Distorts vertices in the NCube by up to +/-.1
//
// Too bad the NCube structure is an array of edges with
// duplicated+ vertices,
// instead of being an array of vertices and an array of double
// indexes to vertices for edges.
// /*
  int i,j,k, match[2];
  int vertmem[MAX_EDGE][2];
  float bump[MAX_DIM];

  for (i=0;i<e;i++)
  {
    vertmem[i][0]=0;
    vertmem[i][1]=0;
  }
  for (i=0;(i<(e-1));i++)
  {
    if (vertmem[i][0]==0) //vert1's pass, and matches (vert1 or 2)
    {
      for (k=0;k<n;k++)
	bump[k]=(float)rand()/(RAND_MAX*5.0)-0.1;
      for (j=i+1;j<e;j++) //search ncube struct for identical vertices
      {
	match[0]=1; match[1]=1;
	if (vertmem[j][0]==0)
	{
	  for (k=0;k<n;k++)
	    if (fabs(ncube[j].vert1[k]-ncube[i].vert1[k])>0.01)
	    match[0] = 0;
	  if (match[0]==1)
	  {
	    vertmem[j][0]=1;
	    for (k=0;k<n;k++)
	      ncube[j].vert1[k]+=bump[k];
	  }
	}
	if (vertmem[j][1]==0)
	{
	  for (k=0;k<n;k++)
	    if (fabs(ncube[j].vert2[k]-ncube[i].vert1[k])>0.01)
	    match[1] = 0;
	  if (match[1]==1)
	  {
	    vertmem[j][1]=1;
	    for (k=0;k<n;k++)
	      ncube[j].vert2[k]+=bump[k];
	  }
	}
      } // matches have been bumped
      for (k=0;k<n;k++)
	ncube[i].vert1[k]+=bump[k]; //vertex bumped
    }
// /*
    if (vertmem[i][1]==0) //vert2's pass, and matches (vert1 or 2)
    {
      for (k=0;k<n;k++)
	bump[k]=(float)rand()/(RAND_MAX*5.0)-0.1;
      for (j=i+1;j<e;j++) //search ncube struct for identical vertices
      {
	match[0]=1; match[1]=1;
	if (vertmem[j][0]==0)
	{
	  for (k=0;k<n;k++)
	    if (fabs(ncube[j].vert1[k]-ncube[i].vert2[k])>0.01)
	    match[0] = 0;
	  if (match[0]==1)
	  {
	    vertmem[j][0]=1;
	    for (k=0;k<n;k++)
	      ncube[j].vert1[k]+=bump[k];
	  }
	}
	if (vertmem[j][1]==0)
	{
	  for (k=0;k<n;k++)
	    if (fabs(ncube[j].vert2[k]-ncube[i].vert2[k])>0.01)
	    match[1] = 0;
	  if (match[1]==1)
	  {
	    vertmem[j][1]=1;
	    for (k=0;k<n;k++)
	      ncube[j].vert2[k]+=bump[k];
	  }
	}
      } // matches have been bumped
      for (k=0;k<n;k++)
	ncube[i].vert2[k]+=bump[k]; //vertex bumped
    }
// */
  }
// */
}

void recube_cont() // recube/contract
// rejoins unattached edge vertices
// draft took all to same point, renamed rec_cont_point()
// earlier draft took all to line, renamed rec_cont_liner()
{
  int i,i2,j;
//  int vertmem[MAX_EDGE][2];
  float diff11[MAX_DIM];
  float diff22[MAX_DIM];
  float diff12[MAX_DIM];
  float diff21[MAX_DIM];
  float diff11p, diff22p; //space metric squared
  float diff12p, diff21p; //space metric squared
  float clsst1, clsst2; //space metric, best so far
  int i2mem1, i2mem2,i2v1v,i2v2v;
// run through each vertex and have it attract the nearest vertex, in space
  for (i=0;i<e;i++)
  {
// search for this edge's vert1 and vert2 neighbors
// Do I even need to work within the edge vert1 vert2 structure when?
  for (i2=0;i2<e;i2++)
  {
    diff11p=0.0; diff22p=0.0; //space metric squared
    diff12p=0.0; diff21p=0.0; //space metric squared
    if((i2!=i))
    {
    for (j=0;j<n;j++)
    {
    diff11[j]=ncube[i].vert1[j]-ncube[i2].vert1[j];
    diff22[j]=ncube[i].vert2[j]-ncube[i2].vert2[j];
    diff12[j]=ncube[i].vert1[j]-ncube[i2].vert2[j];
    diff21[j]=ncube[i].vert2[j]-ncube[i2].vert1[j];
    diff11p+=diff11[j]*diff11[j];
    diff22p+=diff22[j]*diff22[j];
    diff12p+=diff12[j]*diff12[j];
    diff21p+=diff21[j]*diff21[j];
    }
    }
    else
    {
    for (j=0;j<n;j++)
    {
    diff11[j]=ncube[i].vert1[j]-ncube[(i+1)%e].vert1[j];
    diff22[j]=ncube[i].vert2[j]-ncube[(i+1)%e].vert2[j];
    diff12[j]=ncube[i].vert1[j]-ncube[(i+1)%e].vert2[j];
    diff21[j]=ncube[i].vert2[j]-ncube[(i+1)%e].vert1[j];
    diff11p+=diff11[j]*diff11[j];
    diff22p+=diff22[j]*diff22[j];
    diff12p+=diff12[j]*diff12[j];
    diff21p+=diff21[j]*diff21[j];
    }
    }
    if (i2==0)
    {
      if (diff11p<diff12p)
      {clsst1=diff11p;i2mem1=i2;i2v1v=1;} else
      {clsst1=diff12p;i2mem1=i2;i2v1v=2;}
      if (diff22p<diff21p)
      {clsst2=diff22p;i2mem2=i2;i2v2v=2;} else
      {clsst2=diff21p;i2mem2=i2;i2v2v=1;}
    }
    else
    {
      if(diff11p<clsst1)
      {clsst1=diff11p;i2mem1=i2;i2v1v=1;}
      if(diff12p<clsst1)
      {clsst1=diff12p;i2mem1=i2;i2v1v=2;}
      if(diff22p<clsst2)
      {clsst2=diff22p;i2mem2=i2;i2v2v=2;}
      if(diff21p<clsst2)
      {clsst2=diff21p;i2mem2=i2;i2v2v=1;}
    }
  }  // search loop done
//  vector difference / divide by 10 or diffp  add to affected vector
    for (j=0;j<n;j++)
    {
    if(i2mem1!=i) { if(i2v1v==1) {
      ncube[i].vert1[j]-=(ncube[i].vert1[j]-ncube[i2mem1].vert1[j])/(10*(1+clsst1));
      } else {
      ncube[i].vert1[j]-=(ncube[i].vert1[j]-ncube[i2mem1].vert2[j])/(10*(1+clsst1));
      } }
    if(i2mem2!=i) { if(i2v2v==1) {
      ncube[i].vert2[j]-=(ncube[i].vert2[j]-ncube[i2mem2].vert1[j])/(10*(1+clsst2));
      } else {
      ncube[i].vert2[j]-=(ncube[i].vert2[j]-ncube[i2mem2].vert2[j])/(10*(1+clsst2));
      } }
    }
  }
}

float distance(float *v1, float *v2)
{
  int j;
  float dist, distt;
  dist=0.0;
  for (j=0;j<n;j++)
  {
    distt=v1[j]-v2[j];
    dist+=distt*distt;
  }
  return dist;
}

float comp_aver_edge_length() // length in space metric squared
{
  int i;
  float dist=0.0;
  for (i=0;i<e;i++)
  {
    dist+=distance(ncube[i].vert1,ncube[i].vert2)/e;
  }
  return dist;
}

// average vertex-origin distance^2
//  float dist=0.0, origin[MAX_DIM];
//  for (i=0;i<n;i++) origin[i]=0.0;
//  for (i=0;i<e;i++)
//  {
//    dist+=distance(ncube[i].vert1,origin)/e;
//    dist+=distance(ncube[i].vert2,origin)/e;
//  }
//  return dist;
//}

void rec_cont_mas2()
// each vertex and every vertex within .01,.01... of it
//  will get equated and pushed towards the nearest other vertex
//   unless n (n-m?) or more coincident vertices, then no push
// 'B', contrmas2
{
  float dist, bdist[2*MAX_EDGE], eps;
//  int vertmchs[2*MAX_EDGE], vertmatch[2*MAX_EDGE]
  int bdisti[2*MAX_EDGE];
  int i, i2, ix2, i2x2, j, ix2p1, i2x2p1;
// init
  for (i=0;(i<(e+e));i++)
  {
    vertmchs[i]=0; vertmatch[i]=-1; // it has no matches; it isn't a match of another
    bdist[i]=MAXFLOAT; bdisti[i]=0; // easy init; safety;
  }
  eps=.0001*n;
  for (i=0;(i<(e-1));i++)
  {
    ix2=2*i; ix2p1=ix2+1;
    if ((vertmatch[ix2]==-1)){  if ((vertmatch[ix2p1]==-1)){
    dist=distance(ncube[i].vert1,ncube[i].vert2);
    if (dist<eps) {
      vertmchs[ix2]+=1;
      vertmatch[ix2p1]=ix2;
    } else {
	if (dist<bdist[ix2]) {
      bdist[ix2]=dist;
      bdisti[ix2]=ix2p1;    }
	if (dist<bdist[ix2p1]) {
      bdist[ix2p1]=dist;
      bdisti[ix2p1]=ix2;  }
    }}}
    if ((vertmatch[ix2]==-1)||(vertmatch[ix2p1]==-1)){
    for (i2=i+1;i2<e;i2++)
    {
      i2x2=2*i2; i2x2p1=i2x2+1;
      if (vertmatch[ix2]==-1) { if (vertmatch[i2x2]==-1){
      dist=distance(ncube[i].vert1,ncube[i2].vert1);
      if (dist<eps) {
	vertmchs[ix2]+=1;
	vertmatch[i2x2]=ix2;
      } else {
	if (dist<bdist[ix2])
	{
	  bdist[ix2]=dist;
	  bdisti[ix2]=i2x2;
	}
	if (dist<bdist[i2x2])
	{
	  bdist[i2x2]=dist;
	  bdisti[i2x2]=ix2;
	}
      }} if (vertmatch[i2x2p1]==-1){
      dist=distance(ncube[i].vert1,ncube[i2].vert2);
      if (dist<eps) {
	vertmchs[ix2]+=1;
	vertmatch[i2x2p1]=ix2;
      } else {
	if (dist<bdist[ix2])
	{
	  bdist[ix2]=dist;
	  bdisti[ix2]=i2x2p1;
	}
	if (dist<bdist[i2x2p1])
	{
	  bdist[i2x2p1]=dist;
	  bdisti[i2x2p1]=ix2;
	}
      }}}
      if (vertmatch[ix2p1]==-1) { if (vertmatch[i2x2]==-1){
      dist=distance(ncube[i].vert2,ncube[i2].vert1);
      if (dist<eps) {
	vertmchs[ix2p1]+=1;
	vertmatch[i2x2]=ix2p1;
      } else {
	if (dist<bdist[ix2p1])
	{
	  bdist[ix2p1]=dist;
	  bdisti[ix2p1]=i2x2;
	}
	if (dist<bdist[i2x2])
	{
	  bdist[i2x2]=dist;
	  bdisti[i2x2]=ix2p1;
	}
      }} if (vertmatch[i2x2p1]==-1){
      dist=distance(ncube[i].vert2,ncube[i2].vert2);
      if (dist<eps) {
	vertmchs[ix2p1]+=1;
	vertmatch[i2x2p1]=ix2p1;
      } else {
	if (dist<bdist[ix2p1])
	{
	  bdist[ix2p1]=dist;
	  bdisti[ix2p1]=i2x2p1;
	}
	if (dist<bdist[i2x2p1])
	{
	  bdist[i2x2p1]=dist;
	  bdisti[i2x2p1]=ix2p1;
	}
      }}}
    }}
  }
  i=e-1; ix2=2*i; ix2p1=ix2+1;
  if ((vertmatch[ix2]==-1)){  if ((vertmatch[ix2p1]==-1)){
  dist=distance(ncube[i].vert1,ncube[i].vert2);
  if (dist<eps) {
    vertmchs[ix2]+=1;
    vertmatch[ix2p1]=ix2;
  } else {
    if (dist<bdist[ix2])
    {
      bdist[ix2]=dist;
      bdisti[ix2]=ix2p1;
    }
    if (dist<bdist[ix2p1])
    {
      bdist[ix2p1]=dist;
      bdisti[ix2p1]=ix2;
    }
  }}}  // done finding coincidents and nearests
// run through all vertices
  for (i=0;i<e;i++)
  {
    ix2=2*i; ix2p1=ix2+1;
    if (vertmatch[ix2]==-1) { if (vertmchs[ix2]<(n-1))
    {
      if ((bdisti[ix2]&1)==0) {
      for (j=0;j<n;j++)
      {
	ncube[i].vert1[j]-=(ncube[i].vert1[j]-ncube[bdisti[ix2]>>1].vert1[j])/(10*(1+bdist[ix2]));
      } } else {
      for (j=0;j<n;j++)
      {
	ncube[i].vert1[j]-=(ncube[i].vert1[j]-ncube[bdisti[ix2]>>1].vert2[j])/(10*(1+bdist[ix2]));
      } }
    } } else
    {
      if ((vertmatch[ix2]&1)==0) {
      for (j=0;j<n;j++)
      {
	ncube[i].vert1[j]=ncube[vertmatch[ix2]>>1].vert1[j];
      } } else {
      for (j=0;j<n;j++)
      {
	ncube[i].vert1[j]=ncube[vertmatch[ix2]>>1].vert2[j];
      } }
    }
    if (vertmatch[ix2p1]==-1) { if (vertmchs[ix2p1]<(n-1))
    {
      if ((bdisti[ix2p1]&1)==0) {
      for (j=0;j<n;j++)
      {
	ncube[i].vert2[j]-=(ncube[i].vert2[j]-ncube[bdisti[ix2p1]>>1].vert1[j])/(10*(1+bdist[ix2p1]));
      } } else {
      for (j=0;j<n;j++)
      {
	ncube[i].vert2[j]-=(ncube[i].vert2[j]-ncube[bdisti[ix2p1]>>1].vert2[j])/(10*(1+bdist[ix2p1]));
      } }
    } } else
    {
      if ((vertmatch[ix2p1]&1)==0) {
      for (j=0;j<n;j++)
      {
	ncube[i].vert2[j]=ncube[vertmatch[ix2p1]>>1].vert1[j];
      } } else {
      for (j=0;j<n;j++)
      {
	ncube[i].vert2[j]=ncube[vertmatch[ix2p1]>>1].vert2[j];
      } }
    }
  }
}

void cont_outl()
// contract outliers to cube center
// take vector sum of all vertices / 2*e = centroid
// find average distance from vertices to center
// contract toward center vertices more than 1.2*average distance (1.2,1.1,1.07...)
{
  int i,j;
  float dist=0.0, dist2, vsum[MAX_DIM];
  for (j=0;j<n;j++)
    vsum[j]=0.0;
  for (i=0;(i<e);i++)
    for (j=0;j<n;j++)
    {
      vsum[j]+=ncube[i].vert1[j];
      vsum[j]+=ncube[i].vert2[j];
    }
  for (j=0;j<n;j++)
    vsum[j]=vsum[j]/(2.0*(float)e);
  for (i=0;(i<e);i++)
  {
    dist+=distance(vsum,ncube[i].vert1);
    dist+=distance(vsum,ncube[i].vert2);
  }
  dist=dist/(2.0*(float)e);
  for (i=0;(i<e);i++)
  {
    dist2=distance(vsum,ncube[i].vert1);
    if (dist2>(1.07*dist))
    {
      for (j=0;j<n;j++) {
      ncube[i].vert1[j]-=(ncube[i].vert1[j]-vsum[j])/(10*(1+dist2)); }
    }
    dist2=distance(vsum,ncube[i].vert2);
    if (dist2>(1.07*dist))
    {
      for (j=0;j<n;j++) {
      ncube[i].vert2[j]-=(ncube[i].vert2[j]-vsum[j])/(10*(1+dist2)); }
    }
  }
}

void cont_outld(int contoutld)
//   draft of cont_outl(), vert1 or vert2 skewed
// contract outliers to cube center
// take vector sum of all vertices / 2*e = centroid
// find average distance from vertices to center
// contract toward center vertices more than 1.2*average distance (1.2,1.1,1.07...)
{
  int i,j;
  float dist=0.0, dist2, vsum[MAX_DIM];
  for (j=0;j<n;j++)
    vsum[j]=0.0;
  for (i=0;(i<e);i++)
    for (j=0;j<n;j++)
    {
    if (contoutld==1) {
      vsum[j]+=ncube[i].vert1[j];
      vsum[j]+=ncube[i].vert1[j]; }
    else {
      vsum[j]+=ncube[i].vert2[j];
      vsum[j]+=ncube[i].vert2[j]; }
    }
  for (j=0;j<n;j++)
    vsum[j]=vsum[j]/(2.0*(float)e);
  for (i=0;(i<e);i++)
  {
    dist+=distance(vsum,ncube[i].vert1);
    dist+=distance(vsum,ncube[i].vert2);
  }
  dist=dist/(2.0*(float)e);
  for (i=0;(i<e);i++)
  {
    dist2=distance(vsum,ncube[i].vert1);
    if (dist2>(1.07*dist))
    {
      for (j=0;j<n;j++) {
      ncube[i].vert1[j]-=(ncube[i].vert1[j]-vsum[j])/(10*(1+dist2)); }
    }
    dist2=distance(vsum,ncube[i].vert2);
    if (dist2>(1.07*dist))
    {
      for (j=0;j<n;j++) {
      ncube[i].vert2[j]-=(ncube[i].vert2[j]-vsum[j])/(10*(1+dist2)); }
    }
  }
}

void rec_cont_mas2d()
// draft of rec_cont_mas2()
// each vertex and every vertex within .01,.01... of it
//  will get equated and pushed towards the nearest other vertex
//   unless n (n-m?) or more coincident vertices, then no push
// (a pretty or interesting flawed implementation of that)
// 'H', contrmas2d
{
  float dist, bdist[2*MAX_EDGE], eps;
//  int vertmchs[2*MAX_EDGE], vertmatch[2*MAX_EDGE];
  int bdisti[2*MAX_EDGE];
  int i, i2, ix2, i2x2, j, ix2p1, i2x2p1;
//
  for (i=0;(i<(e+e));i++)
  {
    vertmchs[i]=0; vertmatch[i]=-1; // it has no matches; it isn't a match of another
    bdist[i]=MAXFLOAT; bdisti[i]=0; // easy init; safety;
  }
  eps=.0001*n;
  for (i=0;(i<(e-1));i++)
  {
    ix2=2*i; ix2p1=ix2+1;
    if ((vertmatch[ix2]==-1)){  if ((vertmatch[ix2p1]==-1)){
    dist=distance(ncube[i].vert1,ncube[i].vert2);
    if (dist<eps) {
      vertmchs[ix2]+=1;
      vertmatch[ix2p1]=ix2;
    } else {
	if (dist<bdist[ix2]) {
      bdist[ix2]=dist;
      bdisti[ix2]=ix2p1;
      bdist[ix2p1]=dist;
      bdisti[ix2p1]=ix2;  }
    }}}
    if ((vertmatch[ix2]==-1)||(vertmatch[ix2p1]==-1)){
    for (i2=i+1;i2<e;i2++)
    {
      i2x2=2*i2; i2x2p1=i2x2+1;
      if (vertmatch[ix2]==-1) { if (vertmatch[i2x2]==-1){
      dist=distance(ncube[i].vert1,ncube[i2].vert1);
      if (dist<eps) {
	vertmchs[ix2]+=1;
	vertmatch[i2x2]=ix2;
      } else {
	if (dist<bdist[ix2])
	{
	  bdist[ix2]=dist;
	  bdisti[ix2]=i2x2;
	  bdist[i2x2]=dist;
	  bdisti[i2x2]=ix2;
	}
      }} if (vertmatch[i2x2p1]==-1){
      dist=distance(ncube[i].vert1,ncube[i2].vert2);
      if (dist<eps) {
	vertmchs[ix2]+=1;
	vertmatch[i2x2p1]=ix2;
      } else {
	if (dist<bdist[ix2])
	{
	  bdist[ix2]=dist;
	  bdisti[ix2]=i2x2p1;
	  bdist[i2x2p1]=dist;
	  bdisti[i2x2p1]=ix2;
	}
      }}}
      if (vertmatch[ix2p1]==-1) { if (vertmatch[i2x2]==-1){
      dist=distance(ncube[i].vert2,ncube[i2].vert1);
      if (dist<eps) {
	vertmchs[ix2p1]+=1;
	vertmatch[i2x2]=ix2p1;
      } else {
	if (dist<bdist[ix2p1])
	{
	  bdist[ix2p1]=dist;
	  bdisti[ix2p1]=i2x2;
	  bdist[i2x2]=dist;
	  bdisti[i2x2]=ix2p1;
	}
      }} if (vertmatch[i2x2p1]==-1){
      dist=distance(ncube[i].vert2,ncube[i2].vert2);
      if (dist<eps) {
	vertmchs[ix2p1]+=1;
	vertmatch[i2x2p1]=ix2p1;
      } else {
	if (dist<bdist[ix2p1])
	{
	  bdist[ix2p1]=dist;
	  bdisti[ix2p1]=i2x2p1;
	  bdist[i2x2p1]=dist;
	  bdisti[i2x2p1]=ix2p1;
	}
      }}}
    }}
  }
  i=e-1; ix2=2*i; ix2p1=ix2+1;
  if ((vertmatch[ix2]==-1)){  if ((vertmatch[ix2p1]==-1)){
  dist=distance(ncube[i].vert1,ncube[i].vert2);
  if (dist<eps) {
    vertmchs[ix2]+=1;
    vertmatch[ix2p1]=ix2;
  } else {
    if (dist<bdist[ix2])
    {
      bdist[ix2]=dist;
      bdisti[ix2]=ix2p1;
      bdist[ix2p1]=dist;
      bdisti[ix2p1]=ix2;
    }
  }}}  // done finding coincidents and nearests
// run through all vertices
  for (i=0;i<e;i++)
  {
    ix2=2*i; ix2p1=ix2+1;
    if (vertmatch[ix2]==-1) { if (vertmchs[ix2]<(n-1))
    {
      if ((bdisti[ix2]&1)==0) {
      for (j=0;j<n;j++)
      {
	ncube[i].vert1[j]-=(ncube[i].vert1[j]-ncube[bdisti[ix2]>>1].vert1[j])/(10*(1+bdist[ix2]));
      } } else {
      for (j=0;j<n;j++)
      {
	ncube[i].vert1[j]-=(ncube[i].vert1[j]-ncube[bdisti[ix2]>>1].vert2[j])/(10*(1+bdist[ix2]));
      } }
    } } else
    {
      if ((vertmatch[ix2]&1)==0) {
      for (j=0;j<n;j++)
      {
	ncube[i].vert1[j]=ncube[vertmatch[ix2]>>1].vert1[j];
      } } else {
      for (j=0;j<n;j++)
      {
	ncube[i].vert1[j]=ncube[vertmatch[ix2]>>1].vert2[j];
      } }
    }
    if (vertmatch[ix2p1]==-1) { if (vertmchs[ix2p1]<(n-1))
    {
      if ((bdisti[ix2p1]&1)==0) {
      for (j=0;j<n;j++)
      {
	ncube[i].vert2[j]-=(ncube[i].vert2[j]-ncube[bdisti[ix2p1]>>1].vert1[j])/(10*(1+bdist[ix2p1]));
      } } else {
      for (j=0;j<n;j++)
      {
	ncube[i].vert2[j]-=(ncube[i].vert2[j]-ncube[bdisti[ix2p1]>>1].vert2[j])/(10*(1+bdist[ix2p1]));
      } }
    } } else
    {
      if ((vertmatch[ix2p1]&1)==0) {
      for (j=0;j<n;j++)
      {
	ncube[i].vert2[j]=ncube[vertmatch[ix2p1]>>1].vert1[j];
      } } else {
      for (j=0;j<n;j++)
      {
	ncube[i].vert2[j]=ncube[vertmatch[ix2p1]>>1].vert2[j];
      } }
    }
  }
}

void rec_cont_mas()
// rejoins unattached edge vertices, mas
// 'mas' meaning that the rejoining force isn't spent on
//   coincident vertices
{
  int i,i2,j;
//  int vertmem[MAX_EDGE][2];
  float diff11[MAX_DIM];
  float diff22[MAX_DIM];
  float diff12[MAX_DIM];
  float diff21[MAX_DIM];
  float diff11p, diff22p; //space metric squared
  float diff12p, diff21p; //space metric squared
  float clsst1, clsst2; //space metric, best so far
  int i2mem1, i2mem2,i2v1v,i2v2v;
  float eps;
  eps=.001*n;
// run through each vertex and have it attract the nearest vertex, in space
  for (i=0;i<e;i++)
  {
  i2mem1=0;i2mem2=0;i2v1v=1;i2v2v=1;
// search for this edge's vert1 and vert2 neighbors
  for (i2=0;i2<e;i2++)
  {
    diff11p=0.0; diff22p=0.0; //space metric squared
    diff12p=0.0; diff21p=0.0; //space metric squared
    if((i2!=i))
    {
    for (j=0;j<n;j++)
    {
    diff11[j]=ncube[i].vert1[j]-ncube[i2].vert1[j];
    diff22[j]=ncube[i].vert2[j]-ncube[i2].vert2[j];
    diff12[j]=ncube[i].vert1[j]-ncube[i2].vert2[j];
    diff21[j]=ncube[i].vert2[j]-ncube[i2].vert1[j];
    diff11p+=diff11[j]*diff11[j];
    diff22p+=diff22[j]*diff22[j];
    diff12p+=diff12[j]*diff12[j];
    diff21p+=diff21[j]*diff21[j];
    }
    }
    else
    {
    for (j=0;j<n;j++)
    {
    diff11[j]=ncube[i].vert1[j]-ncube[(i+1)%e].vert1[j];
    diff22[j]=ncube[i].vert2[j]-ncube[(i+1)%e].vert2[j];
    diff12[j]=ncube[i].vert1[j]-ncube[(i+1)%e].vert2[j];
    diff21[j]=ncube[i].vert2[j]-ncube[(i+1)%e].vert1[j];
    diff11p+=diff11[j]*diff11[j];
    diff22p+=diff22[j]*diff22[j];
    diff12p+=diff12[j]*diff12[j];
    diff21p+=diff21[j]*diff21[j];
    }
    }
    if (i2==0)
    {
      if ((diff11p<diff12p)&&diff11p>eps)
      {clsst1=diff11p;i2mem1=i2;i2v1v=1;}
      if ((diff12p<diff11p)&&diff12p>eps)
      {clsst1=diff12p;i2mem1=i2;i2v1v=2;}
      if ((diff22p<diff21p)&&diff22p>eps)
      {clsst2=diff22p;i2mem2=i2;i2v2v=2;}
      if ((diff21p<diff22p)&&diff21p>eps)
      {clsst2=diff21p;i2mem2=i2;i2v2v=1;}
    }
    else
    {
      if((diff11p<clsst1)&&diff11p>eps)
      {clsst1=diff11p;i2mem1=i2;i2v1v=1;}
      if((diff12p<clsst1)&&diff12p>eps)
      {clsst1=diff12p;i2mem1=i2;i2v1v=2;}
      if((diff22p<clsst2)&&diff22p>eps)
      {clsst2=diff22p;i2mem2=i2;i2v2v=2;}
      if((diff21p<clsst2)&&diff21p>eps)
      {clsst2=diff21p;i2mem2=i2;i2v2v=1;}
    } // this isn't going to work, re 'mas',unless concident vertices get
      //  bumped together, ie search for matches as in...
      //   'won't work' as in sort of won't work
  }  // search loop done
//  vector difference / divide by 10 or diffp  add to affected vector
    for (j=0;j<n;j++)
    {
    if(i2mem1!=i) { if(i2v1v==1) {
      ncube[i].vert1[j]-=(ncube[i].vert1[j]-ncube[i2mem1].vert1[j])/(10*(1+clsst1));
      } else {
      ncube[i].vert1[j]-=(ncube[i].vert1[j]-ncube[i2mem1].vert2[j])/(10*(1+clsst1));
      } }
    if(i2mem2!=i) { if(i2v2v==1) {
      ncube[i].vert2[j]-=(ncube[i].vert2[j]-ncube[i2mem2].vert1[j])/(10*(1+clsst2));
      } else {
      ncube[i].vert2[j]-=(ncube[i].vert2[j]-ncube[i2mem2].vert2[j])/(10*(1+clsst2));
      } }
    }
  }
}

void rec_cont_point()
// recube_cont() reassembles broken cubes,
// draft took all to same point, renamed rec_cont_point()
// earlier draft took all to line, renamed rec_cont_liner()
//
// revised a bit to not make such small moves in high-d
{
  int i,i2,j;
  float diff11[MAX_DIM];
  float diff22[MAX_DIM];
  float diff12[MAX_DIM];
  float diff21[MAX_DIM];
  float diff11p=0.0, diff22p=0.0; //space metric
  float diff12p=0.0, diff21p=0.0; //space metric
  float clsst1, clsst2; //space metric, best so far
  int i2mem1, i2mem2,i2v1v,i2v2v;
// run through each vertex and have it attract the nearest vertex, in space
  for (i=0;i<e;i++)
  {
// search for this edge's vert1 and vert2 neighbors
// Do I even need to work within the edge vert1 vert2 structure when?
  for (i2=0;i2<e;i2++)
  {
    if((i2!=i))
    {
    for (j=0;j<n;j++)
    {
    diff11[j]=ncube[i].vert1[j]-ncube[i2].vert1[j];
    diff22[j]=ncube[i].vert2[j]-ncube[i2].vert2[j];
    diff12[j]=ncube[i].vert1[j]-ncube[i2].vert2[j];
    diff21[j]=ncube[i].vert2[j]-ncube[i2].vert1[j];
    diff11p+=diff11[j]*diff11[j];
    diff22p+=diff22[j]*diff22[j];
    diff12p+=diff12[j]*diff12[j];
    diff21p+=diff21[j]*diff21[j];
    }
    }
    else
    {
    for (j=0;j<n;j++)
    {
    diff11[j]=ncube[i].vert1[j]-ncube[(i+1)%e].vert1[j];
    diff22[j]=ncube[i].vert2[j]-ncube[(i+1)%e].vert2[j];
    diff12[j]=ncube[i].vert1[j]-ncube[(i+1)%e].vert2[j];
    diff21[j]=ncube[i].vert2[j]-ncube[(i+1)%e].vert1[j];
    diff11p+=diff11[j]*diff11[j];
    diff22p+=diff22[j]*diff22[j];
    diff12p+=diff12[j]*diff12[j];
    diff21p+=diff21[j]*diff21[j];
    }
    }
    if (i2==0)
    {
      if (diff11p<diff12p)
      {clsst1=diff11p;i2mem1=i2;i2v1v=1;} else
      {clsst1=diff12p;i2mem1=i2;i2v1v=2;}
      if (diff22p<diff21p)
      {clsst2=diff22p;i2mem2=i2;i2v2v=2;} else
      {clsst2=diff21p;i2mem2=i2;i2v2v=1;}
    }
    else
    {
      if(diff11p<clsst1)
      {clsst1=diff11p;i2mem1=i2;i2v1v=1;}
      if(diff12p<clsst1)
      {clsst1=diff12p;i2mem1=i2;i2v1v=2;}
      if(diff22p<clsst2)
      {clsst2=diff22p;i2mem2=i2;i2v2v=2;}
      if(diff21p<clsst2)
      {clsst2=diff21p;i2mem2=i2;i2v2v=1;}
    }
  }  // search loop done
//  vector difference / divide by 10 or diffp  add to affected vector
    if(n>4){
      clsst1=sqrt(clsst1);
      clsst2=sqrt(clsst2);}
    for (j=0;j<n;j++)
    {
    if(i2v1v==1) {
      ncube[i].vert1[j]-=(ncube[i].vert1[j]-ncube[i2mem1].vert1[j])/(5*(1+clsst1));
      } else {
      ncube[i].vert1[j]-=(ncube[i].vert1[j]-ncube[i2mem1].vert2[j])/(5*(1+clsst1));
      }
    if(i2v2v==1) {
      ncube[i].vert2[j]-=(ncube[i].vert2[j]-ncube[i2mem2].vert1[j])/(5*(1+clsst2));
      } else {
      ncube[i].vert2[j]-=(ncube[i].vert2[j]-ncube[i2mem2].vert2[j])/(5*(1+clsst2));
      }
    }
  }
}

void rec_cont_liner()
// recube_cont() reassembles broken cubes,
// draft took all to same point, renamed rec_cont_point()
// earlier draft took all to line, renamed rec_cont_liner()
//   revised a bit to not make such small moves in high-d
{
  int i,i2,j2;
//  int vertmem[MAX_EDGE][2];
  float diff1[MAX_DIM];
  float diff2[MAX_DIM];
  float diff1p=0.0, diff2p=0.0; //space metric
  float clsst1, clsst2; //space metric, best so far
  int i2mem1, i2mem2;
//  float bump[MAX_DIM];
// run through each vertex and have it attract the nearest vertex, in space
  for (i=0;i<e;i++)
  {
// search for this edge's vert1 and vert2 neighbors
// Do I even need to work within the edge vert1 vert2 structure when?
  for (i2=0;i2<e;i2++)
  {
    if((i2!=i))
    {
    for (j2=0;j2<n;j2++)
    {
    diff1[j2]=ncube[i].vert1[j2]-ncube[i2].vert1[j2];
    diff2[j2]=ncube[i].vert2[j2]-ncube[i2].vert2[j2];
    diff1p+=diff1[j2]*diff1[j2];
    diff2p+=diff2[j2]*diff2[j2];
    }
//    diff1p=sqrt(diff1p);
//    diff2p=sqrt(diff2p);
    }
    else
    {
    for (j2=0;j2<n;j2++)
    {
    diff1[j2]=ncube[i].vert1[j2]-ncube[(i+1)%e].vert1[j2];
    diff2[j2]=ncube[i].vert2[j2]-ncube[(i+1)%e].vert2[j2];
    diff1p+=diff1[j2]*diff1[j2];
    diff2p+=diff2[j2]*diff2[j2];
    }
    }
    if (i2==0)
    {
      clsst1=diff1p; i2mem1=i2;
      clsst2=diff2p; i2mem2=i2;
    }
    else
    {
      if(diff1p<clsst1)
      {
	clsst1=diff1p; i2mem1=i2;
      }
      if(diff2p<clsst2)
      {
	clsst2=diff2p; i2mem2=i2;
      }
    }
  }  // search loop done
//  vector difference / divide by 10 or diffp  add to affected vector
    if(n>4){
      clsst1=sqrt(clsst1);
      clsst2=sqrt(clsst2);}
    for (j2=0;j2<n;j2++)
    {
    ncube[i].vert1[j2]-=(ncube[i].vert1[j2]-ncube[i2mem1].vert1[j2])/(10*(1+clsst1));
    ncube[i].vert2[j2]-=(ncube[i].vert2[j2]-ncube[i2mem2].vert2[j2])/(10*(1+clsst2));
    }
  }
}

void blow_from_centr()
// expand cube from cube center, by 1/(1+distance^2) force
//
// take vector sum of all vertices / 2*e = centroid
// expand cube from cube center
{
  int i,j;
  float dist2, vsum[MAX_DIM];
  for (j=0;j<n;j++)
    vsum[j]=0.0;
  for (i=0;(i<e);i++)
    for (j=0;j<n;j++)
    {
      vsum[j]+=ncube[i].vert1[j];
      vsum[j]+=ncube[i].vert2[j];
    }
  for (j=0;j<n;j++)
    vsum[j]=vsum[j]/(2.0*(float)e);
  for (i=0;(i<e);i++)
  {
    dist2=distance(vsum,ncube[i].vert1);
    for (j=0;j<n;j++) {
      ncube[i].vert1[j]+=(ncube[i].vert1[j]-vsum[j])/(10*(1+dist2)); }
    dist2=distance(vsum,ncube[i].vert2);
    for (j=0;j<n;j++) {
      ncube[i].vert2[j]+=(ncube[i].vert2[j]-vsum[j])/(10*(1+dist2)); }
  }
}

void blow_astr()
// draft of blow_from_centr(), math mistake
//
{
  int i,j;
  float dist2, vsum[MAX_DIM];
  for (j=0;j<n;j++)
    vsum[j]=0.0;
  for (i=0;(i<e);i++)
    for (j=0;j<n;j++)
    {
      vsum[j]+=ncube[i].vert1[j];
      vsum[j]+=ncube[i].vert2[j];
    }
  for (j=0;j<n;j++)
    vsum[j]=vsum[j]/(2.0*(float)e);
  for (i=0;(i<e);i++)
  {
    dist2=distance(vsum,ncube[i].vert1);
    for (j=0;j<n;j++) {
      ncube[i].vert1[j]+=(ncube[i].vert1[j]-vsum[j])/(10*(1+dist2)); }
    dist2=distance(vsum,ncube[i].vert2);
    for (j=0;j<n;j++) {
      ncube[i].vert2[j]-=(ncube[i].vert2[j]-vsum[j])/(10*(1+dist2)); }
  }
}

void exp_astr()
// expand cube from cube center, which might not be 0,0,0...
// but use math error from  // void blow_astr()
//   draft of blow_from_centr(), math mistake
//
// take vector sum of all vertices / 2*e = centroid
// expand cube from cube center
{
  int i,j;
  float vsum[MAX_DIM];
  for (j=0;j<n;j++)
    vsum[j]=0.0;
  for (i=0;(i<e);i++)
    for (j=0;j<n;j++)
    {
      vsum[j]+=ncube[i].vert1[j];
      vsum[j]+=ncube[i].vert2[j];
    }
  for (j=0;j<n;j++)
    vsum[j]=vsum[j]/(2.0*(float)e);
  for (i=0;(i<e);i++)
  {
    for (j=0;j<n;j++) {
      ncube[i].vert1[j]+=(ncube[i].vert1[j]-vsum[j])*0.087;
      ncube[i].vert2[j]-=(ncube[i].vert2[j]-vsum[j])*0.087;
    }
  }
}

void exp_from_centr()
// expand cube from cube center, which might not be 0,0,0...
//
// take vector sum of all vertices / 2*e = centroid
// expand cube from cube center
{
  int i,j;
  float vsum[MAX_DIM];
  for (j=0;j<n;j++)
    vsum[j]=0.0;
  for (i=0;(i<e);i++)
    for (j=0;j<n;j++)
    {
      vsum[j]+=ncube[i].vert1[j];
      vsum[j]+=ncube[i].vert2[j];
    }
  for (j=0;j<n;j++)
    vsum[j]=vsum[j]/(2.0*(float)e);
  for (i=0;(i<e);i++)
  {
    for (j=0;j<n;j++) {
      ncube[i].vert1[j]+=(ncube[i].vert1[j]-vsum[j])*0.087;
      ncube[i].vert2[j]+=(ncube[i].vert2[j]-vsum[j])*0.087;
    }
  }
}

void knock_left()
{
  int i;
  for (i=0;(i<e);i++)
  {
      ncube[i].vert1[0]-=0.087;
      ncube[i].vert2[0]-=0.087;
  }
}

void knock_cent_left()
// knock verts left if within .5*sqrt(n) of centroid
//
// take vector sum of all vertices / 2*e = centroid
{
  int i,j;
  float vsum[MAX_DIM], dist2, dist;
  for (j=0;j<n;j++)
    vsum[j]=0.0;
  for (i=0;(i<e);i++)
    for (j=0;j<n;j++)
    {
      vsum[j]+=ncube[i].vert1[j];
      vsum[j]+=ncube[i].vert2[j];
    }
  for (j=0;j<n;j++)
    vsum[j]=vsum[j]/(2.0*(float)e);
  dist=.25*n;
  for (i=0;(i<e);i++)
  {
    dist2=distance(vsum,ncube[i].vert1);
    if (dist2<=dist)
    {
      ncube[i].vert1[0]-=0.087;
    }
    dist2=distance(vsum,ncube[i].vert2);
    if (dist2<=dist)
    {
      ncube[i].vert2[0]-=0.087;
    }
  }
}

void drop_pixels()
{
  register int color,x,y;
  for (x=0;x<X_RES;x++)
    for (y=0;y<Y_RES;y++)
      if ((color = get_pixel(x,y)) > 0)
      {
//	put_pixel(x,y,(color-=3)*(color>0));
	color-=3;
	if (color<0) color=0;
	put_pixel(x,y,color);
      }
}

void drop_pixels_xy()
{
  register int color,x,y;
    for (y=0;y<Y_RES;y++)
  for (x=0;x<X_RES;x++)
      if ((color = get_pixel(x,y)) > 0)
      {
//	put_pixel(x,y,(color-=3)*(color>0));
	color-=3;
	if (color<0) color=0;
	put_pixel(x,y,color);
      }
}

void drop_pixels_buff()
{
  unsigned char canbuff[(unsigned)X_RES*(unsigned)Y_RES];
  unsigned int ind_maxp1 = ((unsigned)X_RES*(unsigned)Y_RES);
  unsigned int i, val;
  // buffer the canvas
  for (i=0; i<ind_maxp1; i++) {
    canbuff[i]=canvas[i];
  }
  // drop pixels
  for (i=0; i<ind_maxp1; i++) {
    if (val=canbuff[i]) {
      if (val > 3)
	canbuff[i] = val - 3;
      else
	canbuff[i]=0;
    }
  }
  // canvas the buffer //
  for (i=0; i<ind_maxp1; i++) {
    canvas[i]=canbuff[i];
  }
}

void drop_pixels2()
//try to drop them faster, see drop_pixels()
{
  unsigned int xy;
  for (xy=0;xy<((unsigned)X_RES*(unsigned)Y_RES);xy++)
      if (canvas[xy])
      {
	if (canvas[xy]>3) canvas[xy]-=3;
	else canvas[xy]=0;
      }
}

void drop_pixels3()
//try to drop them faster, see drop_pixels2()
{
  unsigned int xy;
  for (xy=0;xy<(unsigned)X_RES*(unsigned)Y_RES;xy++)
      {
	if (canvas[xy]>3) canvas[xy]-=3;
	else canvas[xy]=0;
      }
}

void drop_pixels3v2()
//try to drop them faster, see drop_pixels2()
{
  unsigned char *p = canvas;
//      if (canvas[xy])
  for (; p<(canvas + (X_RES*Y_RES)); p++)
    if (*p)
    {
      if (*p > 3) *p-=3; else *p=0;
    }
}

void drop_pixels3v3()
//try to drop them faster, see drop_pixels2()
{
  unsigned char *p = canvas;
  unsigned char *p_max = p + (X_RES*Y_RES);
//  register unsigned char val;   /* try a register int here as well */
//  register int val;   /* try a register int here as well */
//  unsigned char val;
  int val;
  for (; p<p_max; p++) {
    if ((val=*p)) {
      if (val > 3)
	*p = val - 3;
      else
	*p=0;
    }
  }
}

void drop_pixels3v()
//try to drop them faster, see drop_pixels2()
{
  unsigned char *p = canvas;
  for (; p<(canvas + ((unsigned)X_RES*(unsigned)Y_RES)); p++) {
   if (*p > 3) *p-=3; else *p=0;
  }
}

void drop_pixels4()
//try to drop them faster, see drop_pixels3()
{
  unsigned int xy;
  for (xy=0;xy<(unsigned)X_RES*(unsigned)Y_RES;xy++)
      {
	if (canvas[xy]>3) canvas[xy]-=3;
	else if (canvas[xy]) canvas[xy]=0;
      }
}

void print_ncube()
{
  int i,j;
//  int p=23;
//  float max=0.0, min=0.0;
//  float max=-MAXFLOAT, min=MAXFLOAT;
  float max, min;
  max = ncube[0].vert1[0];   min = ncube[0].vert1[0];
  for (i=0;i<e;i++)
  {
    printf("Edge %i:",i);
    for (j=0;j<n;j++)
    {
      printf("%f.", ncube[i].vert1[j]);
      if (ncube[i].vert1[j] > max)
	max = ncube[i].vert1[j];
	else
	if (ncube[i].vert1[j] < min)
	  min = ncube[i].vert1[j];
      if (ncube[i].vert2[j] > max)
	max = ncube[i].vert2[j];
	else
	if (ncube[i].vert2[j] < min)
	  min = ncube[i].vert2[j];
    }
    printf(" - ");
    for (j=0;j<n;j++)
    {
      printf("%f.", ncube[i].vert2[j]);
    }
    printf("\n");
/*    if (i>p)
    {
      p=p+23;
      printf("Page. (Hit key)\n");
      pauser();
    }*/
  }
  printf("MaxCoordVal: %f, MinCoordVal: %f\n",max,min);
//  printf("randmax: %u\n",RAND_MAX);
// Spinmodes report, spinmode 2, ang_inc indexing report
/*    for (i=0;i<n-1;i++)
    {
      for (j=i+1;j<n;j++)
      {
	switch (2%4)
	{
	case 0: printf("%i %i %i",i,j,0);
	break;
	case 1: printf("%i %i %i",i,i,(0-j+i) % 4);
	break;
	case 2: printf("%i %i %i",i,j,( (j-i)*(i-j) ) % 6);
	break;
	case 3: printf("%i %i %i",i,j,0);
	break;
	}
      }
    } */

}

void vgapalset(int sml)
{
  int i,j;
  int gcola[3],gcolb[3];
  int vgapal[256][3]; //custom vga palette
  gcola[0]=63;gcola[1]=10;gcola[2]=0;
  gcolb[2]=63;gcolb[0]=0;gcolb[1]=2;
  sml=abs(sml);
  if ((abs(sml)%77)<1)
  {
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[i+16*j][j]=abs(63-4*i);
    vgapal[i+16*j][(j+1)%3]=0;
    vgapal[i+16*j][(j+2)%3]=0;}
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[48+i+16*j][j]=abs(63-4*i);
    vgapal[48+i+16*j][(j+1)%3]=abs(63-4*i);
    vgapal[48+i+16*j][(j+2)%3]=0;}
  for (i=0;i<16;i++) for (j=0;j<3;j++) vgapal[96+i][j]=4*i;
  for (i=0;i<32;i++) for (j=0;j<3;j++) vgapal[112+i][j]=abs(4*i-63);
  for (i=0;i<16;i++) for (j=0;j<3;j++) vgapal[144+i][j]=abs(4*i-63);
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[160+i+16*j][j]=abs(63-4*i);
    vgapal[160+i+16*j][(j+1)%3]=abs(32-2*i);
    vgapal[160+i+16*j][(j+2)%3]=abs(32-2*i);}
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[208+i+16*j][j]=abs(63-4*i);
    vgapal[208+i+16*j][(j+1)%3]=32;
    vgapal[208+i+16*j][(j+2)%3]=0;}
  } else if ((abs(sml)%77)<2)
  {
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[i+16*j][j]=abs(63*(j%2)-4*i);
    vgapal[i+16*j][(j+1)%3]=0;
    vgapal[i+16*j][(j+2)%3]=0;}
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[48+i+16*j][j]=abs(63*((j+1)%2)-4*i);
    vgapal[48+i+16*j][(j+1)%3]=abs(63*((j+1)%2)-4*i);
    vgapal[48+i+16*j][(j+2)%3]=0;}
  for (i=0;i<16;i++) for (j=0;j<3;j++) vgapal[96+i][j]=4*i;
  for (i=0;i<32;i++) for (j=0;j<3;j++) vgapal[112+i][j]=abs(4*i-63);
  for (i=0;i<16;i++) for (j=0;j<3;j++) vgapal[144+i][j]=abs(4*i-63);
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[160+i+16*j][j]=abs(63*(j%2)-4*i);
    vgapal[160+i+16*j][(j+1)%3]=abs(32*(j%2)-2*i);
    vgapal[160+i+16*j][(j+2)%3]=abs(32*(j%2)-2*i);}
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[208+i+16*j][j]=abs(63*((j+1)%2)-4*i);
    vgapal[208+i+16*j][(j+1)%3]=32;
    vgapal[208+i+16*j][(j+2)%3]=0;}
  } else if ((abs(sml)%77)<3)
  {
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[i+16*j][j]=abs(63-4*i);
    vgapal[i+16*j][(j+1)%3]=32;
    vgapal[i+16*j][(j+2)%3]=30;}
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[i+16*j][j]=abs(4*i-63);
    vgapal[i+16*j][(j+1)%3]=32;
    vgapal[i+16*j][(j+2)%3]=30;}
  for (i=0;i<16;i++) for (j=0;j<3;j++) vgapal[96+i][j]=abs(63-4*i);
  for (i=0;i<32;i++) for (j=0;j<3;j++) vgapal[112+i][j]=abs(4*i-63);
  for (i=0;i<16;i++) for (j=0;j<3;j++) vgapal[144+i][j]=abs(4*i-63);
  for (j=0;j<3;j++) vgapal[160][j]=32;
  for (i=161;i<256;i++) for (j=0;j<3;j++)
    vgapal[i][j]=(vgapal[i-1][j]+rand()*8/RAND_MAX)%64;
  } else if ((abs(sml)%77)<4)
  {
  for (j=0;j<3;j++) vgapal[1][j]=32;
  for (i=2;i<256;i++) for (j=0;j<3;j++)
    vgapal[i][j]=(vgapal[i-1][j]+rand()*8/RAND_MAX)%64;
  } else if ((abs(sml)%77)<5)
  {
    for (j=0;j<3;j++) vgapal[1][j]=32;
    for (i=2;i<256;i++)
    {
      gcola[0]=rand()*8/RAND_MAX;
      for (j=0;j<3;j++)
      {
	vgapal[i][j]=vgapal[i-1][j]+gcola[0];
      }
    }
  } else if ((abs(sml)%77)==5)
  {
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[i+16*j][j]=abs(63*(j%2)-4*i);
    vgapal[i+16*j][(j+1)%3]=0;
    vgapal[i+16*j][(j+2)%3]=0;}
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[48+i+16*j][j]=abs(63*((j+1)%2)-4*i);
    vgapal[48+i+16*j][(j+1)%3]=abs(63*((j+1)%2)-4*i);
    vgapal[48+i+16*j][(j+2)%3]=0;}
  for (i=0;i<64;i++) for (j=0;j<3;j++) vgapal[96+i][j]=i;
  gcola[0]=(rand())%3;
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[160+i+16*j][(gcola[0]+j)%3]=abs(63*((j+1)%2)-4*i);
    vgapal[160+i+16*j][(gcola[0]+j+1)%3]=abs(32*((j+1)%2)-2*i);
    vgapal[160+i+16*j][(gcola[0]+j+2)%3]=abs(32*((j+1)%2)-2*i);}
  for (j=0;j<3;j++) for (i=0;i<(16);i++){
    vgapal[208+i+16*j][j]=abs(63*((j)%2)-4*i);
    vgapal[208+i+16*j][(j+1)%3]=32;
    vgapal[208+i+16*j][(j+2)%3]=0;}
  } else if ((abs(sml)%77)==6)
  {
    for(i=1;i<256;i++)
      {     // 0 1 2 - gray
    vgapal[i][ 0       ]=(i/255.0)*gcola[0];
    vgapal[i][ 1       ]=(i/255.0)*gcola[0];
    vgapal[i][ 2       ]=(i/255.0)*gcola[0];
      }
  } else if ((abs(sml)%77)==7)
  {
    for(i=1;i<256;i++)
    { 
      for (j=0;j<3;j++)
      {
	vgapal[i][j]=defvgapal[i][j];
        vgapal[i][j]=defvgapal[i][j];
        vgapal[i][j]=defvgapal[i][j];
      }
    }
  } else if ((abs(sml)%77)<11)
  {
    for (j=0;j<3;j++) vgapal[1][j]=4;
    gcolb[0]=rand()*4.0/RAND_MAX;
    gcola[gcolb[0]%3]=rand()*90/RAND_MAX-18;
    gcola[(1+gcolb[0])%3]=rand()*rand()*gcola[gcolb[0]%3]/RAND_MAX/RAND_MAX;
    gcola[(2+gcolb[0])%3]=rand()*rand()*gcola[gcolb[0]%3]/RAND_MAX/RAND_MAX;
    for (i=2;i<256;i++)
    {
	if ((abs(gcola[0])+abs(gcola[1])+abs(gcola[2]))<1) {
      gcolb[0]=rand()*4.0/RAND_MAX;
      gcola[gcolb[0]%3]=rand()*100/RAND_MAX-18;
      gcola[(1+gcolb[0])%3]=rand()*rand()*gcola[gcolb[0]%3]/RAND_MAX/RAND_MAX;
      gcola[(2+gcolb[0])%3]=rand()*rand()*gcola[gcolb[0]%3]/RAND_MAX/RAND_MAX;    }
      for (j=0;j<3;j++)
      {
	vgapal[i][j]=(vgapal[i-1][j]+gcola[(j+sml)%3])%64;
      }
      for (j=0;j<3;j++)
	gcola[j]/=2;
//	gcola[j]>>=1;
//	gcola[j]=(gcola[j]>0)*gcola[j]/2+(gcola[j]<0)*gcola[j]/2;
    }
  } else
  for(i=1;i<256;i++)
  {
      if ((abs(sml)%77)<19) {
	for(j=0;j<3;j++)
	{
//      vgapal[i][j]=(i/255.0)*gcola[j] + (255-i)/255.0*gcolb[j];
	vgapal[i][(j+(sml%5))%3]=(i/255.0)*gcola[j] + (255-i)/255.0*gcolb[j];
	}
      } else if ((abs(sml)%77)<26)
      {     // 0 1 2 c m y
    vgapal[i][(0+(sml%5))%3]=(i/255.0)*gcola[1] + (255-i)/255.0*gcolb[2];
    vgapal[i][(1+(sml%5))%3]=(i/255.0)*gcola[0] + (255-i)/255.0*gcolb[1];
    vgapal[i][(2+(sml%5))%3]=(i/255.0)*gcola[0] + (255-i)/255.0*gcolb[2];
      } else if ((abs(sml)%77)<34)
      {     // 0 1 2 -  c r  by  mg
    vgapal[i][(0+sml)%3]=(i/255.0)*gcola[2] + (255-i)/255.0*gcolb[2];
    vgapal[i][(1+sml)%3]=(i/255.0)*gcola[0] + (255-i)/255.0*gcolb[0];
    vgapal[i][(2+sml)%3]=(i/255.0)*gcola[0] + (255-i)/255.0*gcolb[1];
      } else if ((abs(sml)%77)<39)
      {     // 0 1 2 -  c g
    vgapal[i][(0+sml)%3]=(i/255.0)*gcola[2] + (255-i)/255.0*gcolb[0];
    vgapal[i][(1+sml)%3]=(i/255.0)*gcola[0] + (255-i)/255.0*gcolb[2];
    vgapal[i][(2+sml)%3]=(i/255.0)*gcola[0] + (255-i)/255.0*gcolb[0];
      } else if ((abs(sml)%77)<45)
      {     // 0 1 2 -  c b
    vgapal[i][(0+sml)%3]=(i/255.0)*gcola[2] + (255-i)/255.0*gcolb[0];
    vgapal[i][(1+sml)%3]=(i/255.0)*gcola[0] + (255-i)/255.0*gcolb[0];
    vgapal[i][(2+sml)%3]=(i/255.0)*gcola[0] + (255-i)/255.0*gcolb[2];
      } else if ((abs(sml)%77)<51)
      {     // 0 1 2 c-gray-black  c cmy
    vgapal[i][(0+sml)%3]=(i/255.0)*gcola[1] + (255-i)/255.0*gcolb[0];
    vgapal[i][(1+sml)%3]=(i/255.0)*gcola[0] + (255-i)/255.0*gcolb[0];
    vgapal[i][(2+sml)%3]=(i/255.0)*gcola[0] + (255-i)/255.0*gcolb[0];
      } else if ((abs(sml)%77)<56)
      {     // 0 1 2 - r-gray-black  r rgb
    vgapal[i][(0+sml)%3]=(i/255.0)*gcola[2] + (255-i)/255.0*gcolb[0];
    vgapal[i][(1+sml)%3]=(i/255.0)*gcola[0] + (255-i)/255.0*gcolb[0];
    vgapal[i][(2+sml)%3]=(i/255.0)*gcola[0] + (255-i)/255.0*gcolb[0];
      } else if ((abs(sml)%77)==56)
      {     // 0 1 2 - gray
    vgapal[i][ 0       ]=(i/255.0)*gcola[0];
    vgapal[i][ 1       ]=(i/255.0)*gcola[0];
    vgapal[i][ 2       ]=(i/255.0)*gcola[0];
      } else if ((abs(sml)%77)<60)
      {     // 0 1 2 -  red-gg-dred ?.
    vgapal[i][(2+sml)%3]=(i/255.0)*gcola[0] + ((i>90)*(i<170))*gcolb[1];
    vgapal[i][(  sml)%3]=(i/255.0)*gcola[2] + (i/255.0)*((i>90)*(i<170))*gcolb[2];
    vgapal[i][(1+sml)%3]=(i/255.0)*gcola[2] + (i/255.0)*((i>90)*(i<170))*gcolb[2];
      } else if ((abs(sml)%77)<63)
      {     // 0 1 2 -  c-green-black
    vgapal[i][(0+sml)%3]=(i/255.0)*gcola[2] + ((i>90)*(i<170))*gcolb[0];
    vgapal[i][(1+sml)%3]=(i/255.0)*gcola[0] + (i/255.0)*((i>90)*(i<170))*gcolb[2];
    vgapal[i][(2+sml)%3]=(i/255.0)*gcola[0] + ((i>90)*(i<170))*gcolb[0];
      } else if ((abs(sml)%77)<66)
      {     // 0 1 2 -  c-red-black
    vgapal[i][(0+sml)%3]=(i/255.0)*gcola[2] + (i/255.0)*((i>90)*(i<170))*gcolb[2];
    vgapal[i][(1+sml)%3]=(i/255.0)*gcola[0] + ((i>90)*(i<170))*gcolb[0];
    vgapal[i][(2+sml)%3]=(i/255.0)*gcola[0] + ((i>90)*(i<170))*gcolb[0];
      } else if ((abs(sml)%77)<69)
      {     // 0 1 2 -  c-blue-black
    vgapal[i][(0+sml)%3]=(i/255.0)*gcola[2] + ((i>90)*(i<170))*gcolb[0];
    vgapal[i][(1+sml)%3]=(i/255.0)*gcola[0] + ((i>90)*(i<170))*gcolb[0];
    vgapal[i][(2+sml)%3]=(i/255.0)*gcola[0] + (i/255.0)*((i>90)*(i<170))*gcolb[2];
      } else if ((abs(sml)%77)<73)
      {     // 0 1 2 -  r-black   rgb
    vgapal[i][(0+sml)%3]=(i/255.0)*gcola[0] ;
    vgapal[i][(1+sml)%3]=(i/255.0)*gcola[2] ;
    vgapal[i][(2+sml)%3]=(i/255.0)*gcola[2] ;
      } else //if ((abs(sml)%77)<77)
      {     // 0 1 2 -  c-black
    vgapal[i][(0+sml)%3]=(i/255.0)*gcola[2] ;
    vgapal[i][(1+sml)%3]=(i/255.0)*gcola[0] ;
    vgapal[i][(2+sml)%3]=(i/255.0)*gcola[0] ;
      }
  }
/*  for(i=127;i<255;i++)
  {
//    vgapal[i][0]=63-(int)((float)i/4.05);
//    vgapal[i][1]=51-(int)((float)i/5.0);
//    vgapal[i][2]=43-(int)((float)i/6.0);
//    vgapal[i][0]=(int)((float)i/4.05);
//    vgapal[i][1]=(int)((float)i/5.0);
//    vgapal[i][2]=(int)((float)i/6.0);
//    vgapal[i][0]=(int)((float)i/5.05);
//    vgapal[i][1]=10+(int)((float)i/6.05);
//    vgapal[i][2]=30+(int)((float)i/8.0);
    vgapal[i][(1+sml)%3]=10-(int)((float)i/32.05);
    vgapal[i][(2+sml)%3]=0+(int)((float)i/8.05);
    vgapal[i][(0+sml)%3]=11+(int)((float)i/5.05);
  }
  for(i=1;i<127;i++)
  {
    vgapal[i][(1+sml)%3]=5-(int)((float)i/32.05);
    vgapal[i][(0+sml)%3]=0+(int)((float)i/14.05);
    vgapal[i][(2+sml)%3]=33-(int)((float)i/7.05);
  } */
  for(j=0;j<3;j++) {vgapal[0][j]=0;}
//  outp(0x03C6,0xff);
  outp(0x03c8, 0);
  for(i=0;i<256;i++)
  {
    for(j=0;j<3;j++)
    {
    outp(0x03c9, vgapal[i][j]);
    }
  }
}

void vertmatcher()
{
// a match is a vertex within distance^2 epsilon of another, epsilon^2=.0001*n
//  extern int vertmchs[], vertmatch[]; // the indexes to vertices are indexed by edge*2 plus 0 or 1
//int vertmchs[2*MAX_EDGE], vertmatch[2*MAX_EDGE];
  int i, i2, ix2, i2x2, j, ix2p1, i2x2p1;
  float dist, eps;
  for (i=0;(i<(e+e));i++) //init
  {
    vertmchs[i]=0; vertmatch[i]=-1; // it has no matches; it isn't a match of another (of a previous)
  }
  eps=.0001*n;
  for (i=0;(i<(e-1));i++)
  {
    ix2=2*i; ix2p1=ix2+1;
    if ((vertmatch[ix2]==-1)){  if ((vertmatch[ix2p1]==-1)){
    dist=distance(ncube[i].vert1,ncube[i].vert2);
    if (dist<eps) {
      vertmchs[ix2]+=1;
      vertmatch[ix2p1]=ix2;
    }
    }}
    if ((vertmatch[ix2]==-1)||(vertmatch[ix2p1]==-1)){
    for (i2=i+1;i2<e;i2++)
    {
      i2x2=2*i2; i2x2p1=i2x2+1;
      if (vertmatch[ix2]==-1) { if (vertmatch[i2x2]==-1){
      dist=distance(ncube[i].vert1,ncube[i2].vert1);
      if (dist<eps) {
	vertmchs[ix2]+=1;
	vertmatch[i2x2]=ix2;
      }
      } if (vertmatch[i2x2p1]==-1){
      dist=distance(ncube[i].vert1,ncube[i2].vert2);
      if (dist<eps) {
	vertmchs[ix2]+=1;
	vertmatch[i2x2p1]=ix2;
      }
      }}
      if (vertmatch[ix2p1]==-1) { if (vertmatch[i2x2]==-1){
      dist=distance(ncube[i].vert2,ncube[i2].vert1);
      if (dist<eps) {
	vertmchs[ix2p1]+=1;
	vertmatch[i2x2]=ix2p1;
      }
      } if (vertmatch[i2x2p1]==-1){
      dist=distance(ncube[i].vert2,ncube[i2].vert2);
      if (dist<eps) {
	vertmchs[ix2p1]+=1;
	vertmatch[i2x2p1]=ix2p1;
      }
      }}
    }}
  }
  i=e-1; ix2=2*i; ix2p1=ix2+1;
  if ((vertmatch[ix2]==-1)){  if ((vertmatch[ix2p1]==-1)){
  dist=distance(ncube[i].vert1,ncube[i].vert2);
  if (dist<eps) {
    vertmchs[ix2]+=1;
    vertmatch[ix2p1]=ix2;
  }
  }}  // done finding coincidents
//returns extern int vertmchs[2*MAX_EDGE], vertmatch[2*MAX_EDGE];
// vertmatch[2e+0or1], if -1 it is not a match of a previous, otherwise is a 2e+0or1 index to that match
// vertmchs[2e+0or1] is the number of matches it has if it is not a match of a previous
}

void init_n_inc(void)
{
  int old_n, old_e; // redundant with n and e over the course of the procedure, but for possible clarity
  int new_e, new_n;
  int i,j,i2,i3;  int diffverts;
  int n2, ix2, ix2p1;

  old_n=n; old_e=e;
  n2 = pow(2,old_n);
  new_e = old_e*2 + n2;
  new_n = old_n+1;

//set incrementing n destination bridge length for n_inc expansion animation
  ninc_targ_bridge_length = comp_aver_edge_length();
//  printf("targ_bridge_length^2: %f\n",ninc_targ_bridge_length);
  if (ninc_targ_bridge_length>6.0) ninc_targ_bridge_length=2.1;

//add -0.05 at coordinate new_n for all the old edges
//add edges for cube at new_n
//duplicate all the old edges with +.05 in the new dimension
  for (i=0;i<old_e;i++)
  {
    i2=i+old_e;
    for (j=0;j<old_n;j++)
    {
      ncube[i2].vert1[j]=ncube[i].vert1[j];
      ncube[i2].vert2[j]=ncube[i].vert2[j];
    }
    ncube[i].vert1[old_n]=-0.11;
    ncube[i].vert2[old_n]=-0.11; // old_n = new_n-1
    ncube[i2].vert1[old_n]=0.11;
    ncube[i2].vert2[old_n]=0.11;
  }

//add more edges for cube at new_n, bridge edges between the old cube and its duplicate
//  run a vertmatch search on the old edges
//  if there are exactly 2^n (old_n) unique verts or else
// //
// run a search of indexes into the old ncube, indexes to vertices (index by edge*2 plus 0 or 1),
//  no two indexes to a vertex equal to another
// count the multiplicity of each of those verteces;
  vertmatcher(); //returns extern int vertmchs[2*MAX_EDGE], vertmatch[2*MAX_EDGE];
// vertmatch[2e+0or1] if -1 it is not a match of a previous, otherwise is a 2e+0or1 index to that match
// vertmchs[2e+0or1] is the number of matches it has if it is not a match of a previous
// //
// if there are exactly 2^n (old_n) different verts or else in the old cube
  diffverts=0;
  for (i=0;(i<(old_e+old_e));i++)
  {
    if (vertmatch[i]==-1) diffverts++;
  }
//  printf("different verts count, old cube: %d\n",diffverts);
  if (diffverts>=n2)  // try a random distribution for the >n2 case in later revision, seniority now
  {
    i3=old_e+old_e;
    for (i=0;i<old_e;i++)
    {
      ix2=2*i; ix2p1=ix2+1;
      i2=i+old_e;
      if (vertmatch[ix2]==-1 && (i3<old_e+old_e+n2))
      {
	for (j=0;j<new_n;j++)
	{
	  ncube[i3].vert1[j] = ncube[i].vert1[j];
	  ncube[i3].vert2[j] = ncube[i2].vert1[j];
	}
        i3++;
      }
      if (vertmatch[ix2p1]==-1 && (i3<old_e+old_e+n2))
      {
	for (j=0;j<new_n;j++)
	{
	  ncube[i3].vert1[j] = ncube[i].vert2[j];
	  ncube[i3].vert2[j] = ncube[i2].vert2[j];
	}
        i3++;
      }
    }
  }
  if (diffverts<n2)
  {
//  printf("different verts count, old cube: %d\n",diffverts);
    i=0;
    i3=old_e+old_e;
    while (i3<old_e+old_e+n2)
    {
      ix2=2*i; ix2p1=ix2+1;
      i2=i+old_e;
      if (vertmatch[ix2]==-1)
      {
	for (j=0;j<new_n;j++)
	{
	  ncube[i3].vert1[j] = ncube[i].vert1[j];
	  ncube[i3].vert2[j] = ncube[i2].vert1[j];
	}
	i3++;
      }
      if ((vertmatch[ix2p1]==-1) && (i3<old_e+old_e+n2))
      {
	for (j=0;j<new_n;j++)
	{
	  ncube[i3].vert1[j] = ncube[i].vert2[j];
	  ncube[i3].vert2[j] = ncube[i2].vert2[j];
	}
	i3++;
      }
      i=(i+1)%old_e;
    }
  }

//init vector to track expansion axis for animated seperation of doubled old cube while rotating in main loop
// extern float n_inc_ref[new_n-1] = 1, rest zero;
  for (j=0;j<old_n;j++)
  {
    n_inc_ref[j]=0;
  }
  n_inc_ref[old_n]=1;

  n=new_n;
  e=new_e;
} //end init_n_inc()

int animate_n_inc(void)
// (part of the n_inc functionality; see n_inc_flag)
//run the bridge edges between the two sub cubes up to the old_cube's average edge_length, one step
//if done then return zero else return one
{
//  extern float ninc_targ_bridge_length; //incrementing n destination bridge length
//  extern float n_inc_ref[]; //incrementing n reference vector for bridge expansion animation
  int i, j, i2;  int ix2,ix2p1,i2x2,i2x2p1;
  int on2; // 2 ^ old n

  if (distance(ncube[e-1].vert1,ncube[e-1].vert2)>=ninc_targ_bridge_length) return 0;

  on2 = pow(2,n-1);
  vertmatcher();
//returns extern int vertmchs[2*MAX_EDGE], vertmatch[2*MAX_EDGE];
// vertmatch[2e+0or1], if -1 it is not a match of a previous, otherwise is a 2e+0or1 index to that match
// vertmchs[2e+0or1] is the number of matches it has if it is not a match of a previous
  for (i=e-1;i>=e-on2;i--)
  {
    for (j=0;j<n;j++)
    {
      ncube[i].vert1[j]-=n_inc_ref[j]/50.0;
      ncube[i].vert2[j]+=n_inc_ref[j]/50.0;
    }
    ix2=2*i; ix2p1=ix2+1;
    if (vertmatch[ix2]!=-1)
    {
      if ((vertmatch[ix2]&1)==0)
      {
	for (j=0;j<n;j++)
	{
	  ncube[vertmatch[ix2]>>1].vert1[j]= ncube[i].vert1[j];
	}
      } else
      {
	for (j=0;j<n;j++)
	{
	  ncube[vertmatch[ix2]>>1].vert2[j]= ncube[i].vert1[j];
	}
      }
      for (i2=vertmatch[ix2]>>1; i2<e-on2; i2++)
      {
	i2x2=i2*2; i2x2p1=i2x2+1;
	if (vertmatch[i2x2]==vertmatch[ix2])
	  for (j=0;j<n;j++)
	    ncube[i2].vert1[j]= ncube[i].vert1[j];
	if (vertmatch[i2x2p1]==vertmatch[ix2])
	  for (j=0;j<n;j++)
	    ncube[i2].vert2[j]= ncube[i].vert1[j];
      }
    }
    if (vertmatch[ix2p1]!=-1)
    {
      if ((vertmatch[ix2p1]&1)==0)
      {
	for (j=0;j<n;j++)
	{
	  ncube[vertmatch[ix2p1]>>1].vert1[j]= ncube[i].vert2[j];
	}
      } else
      {
	for (j=0;j<n;j++)
	{
	  ncube[vertmatch[ix2p1]>>1].vert2[j]= ncube[i].vert2[j];
	}
      }
      for (i2=vertmatch[ix2p1]>>1; i2<e-on2; i2++)
      {
	i2x2=i2*2; i2x2p1=i2x2+1;
	if (vertmatch[i2x2]==vertmatch[ix2p1])
	  for (j=0;j<n;j++)
	    ncube[i2].vert1[j]= ncube[i].vert2[j];
	if (vertmatch[i2x2p1]==vertmatch[ix2p1])
	  for (j=0;j<n;j++)
	    ncube[i2].vert2[j]= ncube[i].vert2[j];
      }
    }
  }
  if (distance(ncube[e-1].vert1,ncube[e-1].vert2)>=ninc_targ_bridge_length)
    return 0;
  else
    return 1;
}

void equi_edge_legth(void)
// spring-like force on edges of length other than the original edge length, 2.0,
// combined on coincident vertices
{
  float ave;
  int i, j, i2;
  int ix2,ix2p1,i2x2,i2x2p1;
  float el;

  for (i=0;i<e+e;i++)
  {
    pushvi[i]=-1;
    for (j=0;j<n;j++) pushes[i][j]=0.0;
  }
  ave =2.0;

  vertmatcher();
//returns extern int vertmchs[2*MAX_EDGE], vertmatch[2*MAX_EDGE];
// vertmatch[2e+0or1], if -1 it is not a match of a previous, otherwise is a 2e+0or1 index to that match
// vertmchs[2e+0or1] is the number of matches it has if it is not a match of a previous

// run through all vertices
// at each vertex and those vertices that are matches of that vertex, handle the edges
  for (i=0;i<e;i++)
  {
    ix2=2*i; ix2p1=ix2+1;
    if (vertmatch[ix2]==-1) // if not a match of a senior vert
    {
      pushvi[ix2]=1;
// sum the edge spring vectors for vertex and all its matches, if any, into pushes[]
      el=distance(ncube[i].vert1,ncube[i].vert2); el = sqrt(el);
// edge spring = ( e / |e| ) * ( 1/2 * (|e|-|a|) ) / 50(animation step factor)
      if (el>0.0003)
      {
	el = ( el - ave ) / (50.0 * el);
	for (j=0;j<n;j++)
	  pushes[ix2][j]=(ncube[i].vert2[j] - ncube[i].vert1[j])*el;
      }
      if (vertmchs[ix2]>0)
      {
/*    if (vertmatch[ix2p1]==ix2)  // this case means this and the previous spring vector are null, but anyway
      {
	el=distance(ncube[i].vert2,ncube[i].vert1); el = sqrt(el);
	if (el>0.0003)
	{
	  el = ( el - ave ) / (50.0 * el);
	  for (j=0;j<n;j++)
	    pushes[ix2][j]+=(ncube[i].vert1[j] - ncube[i].vert2[j])*el;
	}
      }  */
      if (i<e-1) for (i2=i+1;i2<e;i2++)
      {
	i2x2=2*i2; i2x2p1=i2x2+1;
	if (vertmatch[i2x2]==ix2)
	{
	  el=distance(ncube[i2].vert1,ncube[i2].vert2); el = sqrt(el);
	  if (el>0.0003)
	  {
	    el = ( el - ave ) / (50.0 * el);
	    for (j=0;j<n;j++)
	      pushes[ix2][j]+=(ncube[i2].vert2[j] - ncube[i2].vert1[j])*el;
	  }
	}
	if (vertmatch[i2x2p1]==ix2)
	{
	  el=distance(ncube[i2].vert2,ncube[i2].vert1); el = sqrt(el);
	  if (el>0.0003)
	  {
	    el = ( el - ave ) / (50.0 * el);
	    for (j=0;j<n;j++)
	      pushes[ix2][j]+=(ncube[i2].vert1[j] - ncube[i2].vert2[j])*el;
	  }
	}
      }
      }
    }
    if (vertmatch[ix2p1]==-1) // if not a match of a senior vert
    {
      pushvi[ix2p1]=1;
// sum the edge spring vectors for it and all its matches, if any, into pushes[]
      el=distance(ncube[i].vert2,ncube[i].vert1); el = sqrt(el);
// edge spring = ( e / |e| ) * ( 1/2 * (|e|-|a|) ) / 50(animation step factor)
      if (el>0.0003)
      {
	el = ( el - ave ) / (50.0 * el);
	for (j=0;j<n;j++)
	  pushes[ix2p1][j]=(ncube[i].vert1[j] - ncube[i].vert2[j])*el;
      }
      if (vertmchs[ix2p1]>0)
      if (i<e-1) for (i2=i+1;i2<e;i2++)
      {
	i2x2=2*i2; i2x2p1=i2x2+1;
	if (vertmatch[i2x2]==ix2p1)
	{
	  el=distance(ncube[i2].vert1,ncube[i2].vert2); el = sqrt(el);
	  if (el>0.0003)
	  {
	    el = ( el - ave ) / (50.0 * el);
	    for (j=0;j<n;j++)
	      pushes[ix2p1][j]+=(ncube[i2].vert2[j] - ncube[i2].vert1[j])*el;
	  }
	}
	if (vertmatch[i2x2p1]==ix2p1)
	{
	  el=distance(ncube[i2].vert2,ncube[i2].vert1); el = sqrt(el);
	  if (el>0.0003)
	  {
	    el = ( el - ave ) / (50.0 * el);
	    for (j=0;j<n;j++)
	      pushes[ix2p1][j]+=(ncube[i2].vert1[j] - ncube[i2].vert2[j])*el;
	  }
	}
      }
    }
  }
//  pauser();
  for (i=0;i<e;i++)
  {
    ix2=2*i; ix2p1=ix2+1;
    if (pushvi[ix2]==1)
    {
      for (j=0;j<n;j++)
	ncube[i].vert1[j]+=pushes[ix2][j];
      if (vertmchs[ix2]>0)
      {
      if (vertmatch[ix2p1]==ix2)
      {
	for (j=0;j<n;j++)
//	  ncube[i].vert2[j]+=pushes[ix2][j];
	  ncube[i].vert2[j]=ncube[i].vert1[j];
      }
      if (i<e-1) for (i2=i+1;i2<e;i2++)
      {
	i2x2=2*i2; i2x2p1=i2x2+1;
	if (vertmatch[i2x2]==ix2)
	{
	  for (j=0;j<n;j++)
//	    ncube[i2].vert1[j]+=pushes[ix2][j];
	    ncube[i2].vert1[j]=ncube[i].vert1[j];
	}
	if (vertmatch[i2x2p1]==ix2)
	{
	  for (j=0;j<n;j++)
//	    ncube[i2].vert2[j]+=pushes[ix2][j];
	    ncube[i2].vert2[j]=ncube[i].vert1[j];
	}
      }
      }
    }
    if (pushvi[ix2p1]==1)
    {
      for (j=0;j<n;j++)
	ncube[i].vert2[j]+=pushes[ix2p1][j];
      if (vertmchs[ix2p1]>0)
      if (i<e-1) for (i2=i+1;i2<e;i2++)
      {
	i2x2=2*i2; i2x2p1=i2x2+1;
	if (vertmatch[i2x2]==ix2p1)
	{
	  for (j=0;j<n;j++)
//	    ncube[i2].vert1[j]+=pushes[ix2p1][j];
	    ncube[i2].vert1[j]=ncube[i].vert2[j];
	}
	if (vertmatch[i2x2p1]==ix2p1)
	{
	  for (j=0;j<n;j++)
//	    ncube[i2].vert2[j]+=pushes[ix2p1][j];
	    ncube[i2].vert2[j]=ncube[i].vert2[j];
	}
      }
    }
  }
}

void swappc6(float *v6_2)
{
  float temp;
  int i;
  for (i=0;i<6;i++)
  {
    temp=palcol[i].b;palcol[i].b=v6_2[i];v6_2[i]=temp;
  }
  return;
}

void matmult(void)
// rot = rotm * rot;
{
  int i, j, k;
  float result[MAX_DIM][MAX_DIM];
  for (i=0;i<n;i++)
  {
    for (j=0;j<n;j++)
    {
      result[i][j]=0.0;
      for (k=0;k<n;k++)
      {
        result[i][j]+=rotm[i][k]*rot[k][j];
      }
    }
  }
  for (i=0;i<n;i++)
  {
    for (j=0;j<n;j++)
    {
      rot[i][j]=result[i][j];
    }
  }
}

void set_rot_mat_m(float a, int b1, int b2)
{
  int i;
  memset(rotm,0,mat_size);
  for (i=0;i<n;i++)
    rotm[i][i] = 1;
  rotm[b1][b1] = cos(a);
  rotm[b1][b2] = sin(a);
  rotm[b2][b1] = -sin(a);
  rotm[b2][b2] = cos(a);
}

void rotatem(float a, int b1, int b2)
{
  int i;
  set_rot_mat_m(a, b1, b2);
  matmult(); // mat. mult., rot = rotm * rot
  //rotation of cube elsewhere after building up rot with this
}

void spinsm(int spin, int *spinmodes, int *frot_change_flag, float ang_inc[])
{
  int i, j, ji, ij;
  static int old_spin=-1, old_spinmode=-1;

  if ( ( spin != old_spin ) || ( *spinmodes != old_spinmode ) || ( *frot_change_flag) )
  {
  *frot_change_flag = 0;
  old_spin = spin;
  old_spinmode = *spinmodes;
  memset(rot,0,mat_size);
  for (i=0;i<n;i++)
    rot[i][i] = 1;
  if (*spinmodes >= 3) *spinmodes = 0;
  if (*spinmodes == 0)
    switch (spin%8) { case 0: case 1: case 2: case 3:
    for (i=n-2;i>=0;i--)
    {
      for (j=n-1;j>=i+1;j--)
      {
	switch (spin%4)
	{
	case 0: rotatem(ang_inc[0], i, j); break;
	case 1: rotatem(RAD((float)(i+j+2)/10), i, j); break;
	case 2: rotatem(ang_inc[(i+j) % 4], i, j); break;
	case 3: rotatem(ang_inc[(i+j) % 6], i, j); break; //
	}
      }
    } break;
    case 4:
      if (n>3) rotatem(ang_inc[3], 1, 3); 
      if (n>4) rotatem(ang_inc[1], 2, 4);
      if (n>5) rotatem(ang_inc[1], 0, 5);
      if (n>1) rotatem(ang_inc[1], 0, 1); break;
    case 5: // spin from old ncube'6'
      if (n>3) rotatem(ang_inc[3], 0, 3); 
      if (n>4) rotatem(ang_inc[1], 3, 4);
      if (n>5) rotatem(ang_inc[1], 4, 5);
      if (n>2) rotatem(ang_inc[1], 1, 2); break;
    case 6:
      if (n>3) rotatem(ang_inc[3], 1, 3); 
      if (n>4) rotatem(ang_inc[1], 2, 4);
      if (n>5) rotatem(ang_inc[1], 2, 5);
      if (n>1) rotatem(ang_inc[0], 0, 1); break;
    case 7: break; }
  if (*spinmodes == 1)
    for (i=n-2;i>=0;i--)
    {
      for (j=n-1;j>=i+1;j--)
      {
	switch (spin%8)
	{
	case 0: rotatem(ang_inc[0]*(0-j+i), i, j); break;
	case 1: rotatem(ang_inc[abs((0-j+i) % 4)]*((0-j+i) % 4), i, j); break;
	case 2: rotatem(ang_inc[abs(( (j-i)*(i-j) ) % 6)]*(( (j-i)*(i-j) ) % 6), i, j); break;
	case 3: /*rotatem(ang_inc[sml]*(0-j+i), i, j);*/ break;
	case 4: rotatem(ang_inc[(i+j)%5]*(0-j+i), i, j); break;
	case 5: ji=(i-j)%4;rotatem(ang_inc[ji]*(ji>=0),i,j); break;
	case 6: ij=((j-i)*(i-j))%6;rotatem(ang_inc[ij]*(ij>=0), i, j); break;
	case 7: break;
	}
      }
    }
  if (*spinmodes == 2)
	switch (spin%8)
	{
	case 0: for (i=n-2;i>=0;i--) for (j=n-1;j>=i+1;j--) rotatem(ang_inc[(i+j)%5]*50, i, j ); break;
	case 1: for (i=n-2;i>=0;i--) for (j=n-1;j>=i+1;j--) if (n>2)
	    rotatem(ang_inc[(i+j)%5], (i%(n-2)), (j%(n-2)) );
	    else rotatem(ang_inc[(i+j)%5], i, j ); break;
	case 2: for (i=n-2;i>=0;i--) for (j=n-1;j>=i+1;j--) rotatem(ang_inc[(i+j)%5], (i%(n-1)), (j%(n-1)) ); break;
	case 3: // spin from old ncube'6' // may produce false rotation matrices for n<6
	    rotatem(ang_inc[3], 0, 3*(n>3));
	    rotatem(ang_inc[1], 3*(n>3), 4*(n>4));
	    rotatem(ang_inc[1], 4*(n>4), 5*(n>5));
	    rotatem(ang_inc[1], 1, 2*(n>2)); break;
	case 4: for (i=n-2;i>=0;i--) for (j=n-1;j>=i+1;j--) if((j<(n-1))&&(i<(n-1))) rotatem(ang_inc[(i+j)%5],i,j);break;
	case 5: for (i=n-2;i>=0;i--) for (j=n-1;j>=i+1;j--) if((j<(n-2))&&(i<(n-2))) rotatem(ang_inc[(i+j)%5],i,j);break;
	case 6: for (i=n-2;i>=0;i--) for (j=n-1;j>=i+1;j--) if((j<(n-1))&&(i<(n-2))) rotatem(ang_inc[(i+j)%6],i,j);break;
	case 7: break;
	}  // cases 1, 2, and 3 may cause false rotation matrices
  }
  for (i=0;i<e;i++)
    rotate_edge(i);
  //also rotate the reference vector for n_inc
  rotate_point(n_inc_ref);
}

void spins(int spin, int *spinmodes, float ang_inc[])
{
  int i, j, ji, ij;

  if (*spinmodes >= 3) *spinmodes = 0;
  if (*spinmodes == 0)
    switch (spin%8) { case 0: case 1: case 2: case 3:
    for (i=0;i<n-1;i++)
    {
      for (j=i+1;j<n;j++)
      {
	switch (spin%4)
	{
	case 0: rotate(ang_inc[0], i, j); break;
	case 1: rotate(RAD((float)(i+j+2)/10), i, j); break;
	case 2: rotate(ang_inc[(i+j) % 4], i, j); break;
	case 3: rotate(ang_inc[(i+j) % 6], i, j); break; //
	}
      }
    } break;
    case 4:
      if (n>1) rotate(ang_inc[1], 0, 1);
      if (n>5) rotate(ang_inc[1], 0, 5);
      if (n>4) rotate(ang_inc[1], 2, 4);
      if (n>3) rotate(ang_inc[3], 1, 3); break;
    case 5: // spin from old ncube'6'
      if (n>2) rotate(ang_inc[1], 1, 2);
      if (n>5) rotate(ang_inc[1], 4, 5);
      if (n>4) rotate(ang_inc[1], 3, 4);
      if (n>3) rotate(ang_inc[3], 0, 3); break;
    case 6:
      if (n>1) rotate(ang_inc[0], 0, 1);
      if (n>5) rotate(ang_inc[1], 2, 5);
      if (n>4) rotate(ang_inc[1], 2, 4);
      if (n>3) rotate(ang_inc[3], 1, 3); break;
    case 7: break; }
  if (*spinmodes == 1)
    for (i=0;i<n-1;i++)
    {
      for (j=i+1;j<n;j++)
      {
	switch (spin%8)
	{
	case 0: rotate(ang_inc[0]*(0-j+i), i, j); break;
	case 1: rotate(ang_inc[abs((0-j+i) % 4)]*((0-j+i) % 4), i, j); break;
	case 2: rotate(ang_inc[abs(( (j-i)*(i-j) ) % 6)]*(( (j-i)*(i-j) ) % 6), i, j); break;
	case 3: /*rotate(ang_inc[sml]*(0-j+i), i, j);*/ break;
	case 4: rotate(ang_inc[(i+j)%5]*(0-j+i), i, j); break;
	case 5: ji=(i-j)%4;rotate(ang_inc[ji]*(ji>=0),i,j); break;
	case 6: ij=((j-i)*(i-j))%6;rotate(ang_inc[ij]*(ij>=0), i, j); break;
	case 7: break;
	}
      }
    }
  if (*spinmodes == 2)
	switch (spin%8)
	{
	case 0: for (i=0;i<n-1;i++) for (j=i+1;j<n;j++) rotate(ang_inc[(i+j)%5]*50, i, j ); break;
	case 1: for (i=0;i<n-1;i++) for (j=i+1;j<n;j++) if (n>2)
	    rotate(ang_inc[(i+j)%5], (i%(n-2)), (j%(n-2)) );
	    else rotate(ang_inc[(i+j)%5], i, j ); break;
	case 2: for (i=0;i<n-1;i++) for (j=i+1;j<n;j++) rotate(ang_inc[(i+j)%5], (i%(n-1)), (j%(n-1)) ); break;
	case 3: // spin from old ncube'6' // may produce false rotation matrices for n<6
	    rotate(ang_inc[1], 1, 2*(n>2));
	    rotate(ang_inc[1], 4*(n>4), 5*(n>5));
	    rotate(ang_inc[1], 3*(n>3), 4*(n>4));
	    rotate(ang_inc[3], 0, 3*(n>3)); break;
	case 4: for (i=0;i<n-1;i++) for (j=i+1;j<n;j++) if((j<(n-1))&&(i<(n-1))) rotate(ang_inc[(i+j)%5],i,j);break;
	case 5: for (i=0;i<n-1;i++) for (j=i+1;j<n;j++) if((j<(n-2))&&(i<(n-2))) rotate(ang_inc[(i+j)%5],i,j);break;
	case 6: for (i=0;i<n-1;i++) for (j=i+1;j<n;j++) if((j<(n-1))&&(i<(n-2))) rotate(ang_inc[(i+j)%6],i,j);break;
	case 7: break;
	}  // cases 1, 2, and 3 may cause false rotation matrices
}

void main()
{
  int ch=2; int delnum=20;
  int col=16, i, j, wipe, quit, spin, nd=0, vertsy=0;
  int spinmodes=0, sml=83;
  int contract=0, contline=0, contpoint=0, contrmas=0, contrmas2=0, contrmas2d=0;
  int contoutl=0, contoutld=0;  int n_inc_flag=0;  int equiel=0;
  int old_draw_flag=0, up=1, faster_rot=0, faster_rot_changed=0;
  char reduce=0; // flag for toggle of reduce old pixels procedure, drop_pixels()
//  float ang[MAX_DIM];
  float ang_inc[MAX_DIM];
  float pcbalts[2][6]; //palcol[].b alts, mid-centers, high centers
  int plclbalt=0;

  canvas = SCREEN;

  printf("Pick a number 2 to %i.\n",MAX_DIM);
  while ((nd<2) || (nd>MAX_DIM))
  {
    while (!kbhit());
    nd = (getch() - 48);
  }
  init_dimension(nd);  // choose n here

  vga_mode(0x13);

  wipe=1; quit=0; spin=0; colpal=0;
  ang_inc[0] = RAD(.05*5);
  ang_inc[1] = RAD(.07*5);
  ang_inc[2] = RAD(.04*5);
  ang_inc[3] = RAD(.03*5);
  ang_inc[4] = RAD(.10*5);
  ang_inc[5] = RAD(.13*5);
  palcol[0].a= 7.5;  palcol[0].b= 23.5;
  palcol[1].a= 15.5; palcol[1].b= 34.5;
  palcol[2].a= 30.5; palcol[2].b= 34.5;
  palcol[3].a= 60.5; palcol[3].b= 64.5;
  palcol[4].a= 63.5; palcol[4].b= 190.5;
  palcol[5].a= 120.5; palcol[5].b= 127.5;
  for(i=0;i<6;i++) pcbalts[0][i]=127.5;
  for(i=0;i<6;i++) pcbalts[1][i]=255-palcol[i].b;
  outp(0x03c7, 0); //saving the default vga palette
  for(i=0;i<256;i++)
  {
    for(j=0;j<3;j++)
    {
      defvgapal[i][j]=inp(0x03c9);
    }
  } //done saving the default vga palette
  for(j=0;j<MAX_DIM;j++) n_inc_ref[j]=0.0;
  while (quit!=1)
  {
    if (faster_rot==0)
      spins(spin,&spinmodes,ang_inc);
    if (faster_rot==1)
    {
      spinsm(spin,&spinmodes,&faster_rot_changed,ang_inc);
//      faster_rot_changed=1;
    }
//    clear_screen();
    if (wipe==1) { if (vertsy) { VERT_SYNC2; } wipe_prev_cube();}
    if (vertsy) { VERT_SYNC1; }
    if (old_draw_flag)
    {
      old_draw_ncube(col);
      col+=up;
//      if (col>=31)   {    up=-1;   }
//      if (col<=16)   {    up=1;    }
      if (col>=palcol[colpal].b+palcol[colpal].a) {up=-1;}
      if (col<=palcol[colpal].b-palcol[colpal].a) {up=1; }
    }
    else
      draw_ncube(10);
    delay( delnum );
    if (kbhit())
      {
	ch=getch();
	if (ch=='w') wipe= !wipe;
	else
	if (ch=='z') old_draw_flag= !old_draw_flag;
	else
	if (ch=='r') reduce= !reduce;
	else
	if (ch=='c') clear_screen();
	else
	if (ch=='v') vertsy=!vertsy;
	else
	if (ch=='s') spin=spin+1;
	else
	if (ch=='S') delnum=(delnum+1)%28; else if(ch=='!') delnum=0;
	else
	if (ch=='M') spinmodes = spinmodes +1;
	else
	if (ch=='f') { faster_rot = (faster_rot+1)%2; faster_rot_changed = 1; }
	else
	if (ch=='A') sml += 1;
	else
	if (ch=='Q') sml -= 1;
	else
	if (ch=='W') sml = 82;
	else
	if (ch=='d') distort();
	else
	if (ch=='b') breakdist();
	else
	if (ch=='a') shrink();
	else
	if (ch=='t') grow();
	else
	if (ch=='e') exp_from_centr();
	else
	if (ch=='E') blow_from_centr();
	else
	if (ch=='Y') blow_astr();
	else
	if (ch=='y') exp_astr();
	else
	if (ch=='R') knock_left();
	else
	if (ch=='X') knock_cent_left();
	else
	if (ch=='N') {  if (n<MAX_DIM) {  n_inc_flag=1; init_n_inc(); faster_rot_changed = 1; equiel=0;
	  contract=0;contline=0;contpoint=0;contrmas=0;contrmas2=0;contrmas2d=0;contoutl=0;contoutld=0;} else pauser(); }
	else
	if (ch=='C'&&n_inc_flag==0) {contract=~contract+2;contline=0;contpoint=0;contrmas=0;contrmas2=0;contrmas2d=0;contoutl=0;contoutld=0;equiel=0;}
	else
	if (ch=='D'&&n_inc_flag==0) {contline=~contline+2;contract=0;contpoint=0;contrmas=0;contrmas2=0;contrmas2d=0;contoutl=0;contoutld=0;equiel=0;}
	else
	if (ch=='F'&&n_inc_flag==0) {contpoint=~contpoint+2;contract=0;contline=0;contrmas=0;contrmas2=0;contrmas2d=0;contoutl=0;contoutld=0;equiel=0;}
	else
	if (ch=='V'&&n_inc_flag==0) {contrmas=~contrmas+2;contract=0;contline=0;contpoint=0;contrmas2=0;contrmas2d=0;contoutl=0;contoutld=0;equiel=0;}
	else
	if (ch=='B'&&n_inc_flag==0) {contrmas2=~contrmas2+2;contract=0;contline=0;contpoint=0;contrmas=0;contrmas2d=0;contoutl=0;contoutld=0;equiel=0;}
	else
	if (ch=='H'&&n_inc_flag==0) {contrmas2d=~contrmas2d+2;contract=0;contline=0;contpoint=0;contrmas=0;contrmas2=0;contoutl=0;contoutld=0;equiel=0;}
	else
	if (ch=='G'&&n_inc_flag==0) {contoutl=~contoutl+2;contract=0;contline=0;contpoint=0;contrmas=0;contrmas2=0;contrmas2d=0;contoutld=0;equiel=0;}
	else
	if (ch=='g'&&n_inc_flag==0) {contoutld=(contoutld+1)%3;contract=0;contline=0;contpoint=0;contrmas=0;contrmas2=0;contrmas2d=0;contoutl=0;equiel=0;}
	else
	if (ch=='T'&&n_inc_flag==0) {equiel=!equiel;contract=0;contline=0;contpoint=0;contrmas=0;contrmas2=0;contrmas2d=0;contoutl=0;contoutld=0;}
	else
	if (ch=='p') colpal=(colpal+1)%6;
	else
	if (ch=='P') vgapalset(sml);
	else
	if (ch=='\'') {swappc6(pcbalts[plclbalt]);plclbalt=!plclbalt;}
	else
	if (ch=='n')
	{
	  printf("Pick a number 2 to %i.\n",MAX_DIM);
	  nd =0;
	  while ((nd<2) || (nd>MAX_DIM))
	  {
	    while (!kbhit());
	    nd = (getch() - 48);
	  }
	  init_dimension(nd);  // choose n here
	  contract=0;contline=0;contpoint=0;contrmas=0;contrmas2=0;contrmas2d=0;
	  contoutl=0;contoutld=0; n_inc_flag=0; equiel=0;
	  faster_rot_changed = 1;
	}
	else
	if (ch==0x1b) quit=1;
	else
	pauser();
      }
    if (reduce) {
      switch (abs(sml%4))
      {
	case 0: drop_pixels(); break;
	case 1: drop_pixels_buff(); break;
	case 2: drop_pixels3v3(); break;
	case 3: drop_pixels_xy(); break;
//	case 3: drop_pixels3v(); break;
//	case 4: drop_pixels2(); break;
//	case 5: drop_pixels4(); break;
      }      }
    if (n_inc_flag==1)  {  n_inc_flag = animate_n_inc();  }
    if (contract!=0) { recube_cont(); /*contract+=1;*/ }
    if (contline!=0) { rec_cont_liner(); }
    if (contpoint!=0) { rec_cont_point(); }
    if (contrmas!=0) { rec_cont_mas(); }
    if (contrmas2!=0) { rec_cont_mas2(); }
    if (contrmas2d!=0) { rec_cont_mas2d(); }
    if (contoutl!=0) { cont_outl(); }
    if (contoutld!=0) { cont_outld(contoutld); }
    if (equiel!=0) { equi_edge_legth(); }
  }
  vga_mode(0x03);
  print_ncube();
  pauser();
}