/* ========================================================================== */
/*                                                                            */
/*   ecox.c                                                                   */
/*   (c) 2004 Mauro Grassi                                                    */
/*                                                                            */
/*   Description: Coxeter group program with OpenGL...                        */
/*                including Groupoid viewer...                                */  
/*                including positive semidefinite types...                    */
/*                including hyperbolic types...                       */
/*                Version 1.0 4th November 2004                               */
/*                including presentation of N_J                               */
/* ========================================================================== */
#include <stdlib.h> 
#include <math.h>
#include <string.h>
#include <stdarg.h>
#include <stdio.h>
#include <gl/glut.h>

// Coxeter group types...
#define NONE -1
#define A 1
#define B 2
#define D 3
#define E 4
#define F 5
#define H 6 
#define I 7
#define X 8  // infinite types here...
#define C 9
#define G 10
#define HYPBASE 11 // base for hyperbolic groups...
#define CUSTOM 83
#define CUSTOMRANK 10
// Main program parameters...
#define MAXRANK 20         // defines the maximum number of nodes...
#define MAXSELECT 128
#define MAXNODES 1000      // defines the maximum number of nodes in groupoid...
#define MAXRELATIONLEN 500 // maximum length of relation...
#define MAXRELATIONS 10000   // maximum number of relations...
#define MAXDICT MAXRANK*MAXNODES // maximum number of labels...
#define MAXTRANSLATIONS 10000// maximum number of translations to avoid infinite loops...

// data for hyperbolic groups, split into ranks...
int hyperbolic=0;  // index to current one...
int hyprank[72]={ 4,4,4,4,4,4,4,4,4,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,8,8,8,8,9,9,9,9,10,10,10 };
float hyppoints[72][10][3]={
// 1 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 2 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 3 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 4 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 5 - a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 6 - a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 7 - a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 8 - a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 9 - a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 10 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 11 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 12 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 13 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 14 - trapezoid
0.0,0.0,0.0,1.0,1.0,0.0,2.0,1.0,0.0,3.0,0.0,0.0,1.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 15- a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 16- a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 17- a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 18- a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 19- a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 20- a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 21- a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 22- a square
0.0,0.0,0.0,0.0,1.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 23 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 24 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 25 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 26 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 27 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 28 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 29 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 30 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 31 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 32 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 33 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 34 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 35 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 36 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 37 - triangle with midpoint
0.0,0.0,0.0,1.0,1.0,0.0,2.0,0.0,0.0,1.0,0.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 38 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 39 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 40 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 41 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 42 - diamond with line...
0.0,0.0,0.0,1.0,0.0,0.0,2.0,1.0,0.0,3.0,0.0,0.0,2.0,-1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 43 - diamond with line...
0.0,0.0,0.0,1.0,0.0,0.0,2.0,1.0,0.0,3.0,0.0,0.0,2.0,-1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 44 - four pointed cross with line...
0.0,0.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,1.0,-1.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,
// 45 - four pointed cross with line...
0.0,0.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,1.0,-1.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,
// 46 - trapezoid
0.0,0.0,0.0,1.0,1.0,0.0,2.0,1.0,0.0,3.0,0.0,0.0,1.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 47 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 48 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 49 - a line
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0,
// 50 - a fork 
0.0,1.0,0.0,0.0,-1.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 51 - E6 type
2.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 52 - E6 type
2.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 53 - four pointed cross with line...
0.0,0.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,1.0,-1.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,
// 54 - four pointed cross with line...
0.0,0.0,0.0,1.0,1.0,0.0,1.0,0.0,0.0,1.0,-1.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,
// 55 - five pointed star...
0.0,0.0,0.0,1.0,1.0,0.0,2.0,0.0,0.0,1.5,-1.0,0.0,0.5,-1.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
// 56 - big diamond with line
0.0,0.0,0.0,1.0,1.0,0.0,2.0,1.0,0.0,1.0,-1.0,0.0,2.0,-1.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,
// 57 - big diamond with line
0.0,0.0,0.0,1.0,1.0,0.0,2.0,1.0,0.0,1.0,-1.0,0.0,2.0,-1.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,
// 58 - chimney with line
0.0,1.0,0.0,1.0,1.0,0.0,0.0,-1.0,0.0,1.0,-1.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,
// 59 - E6 type
2.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 60 - double E6 type
1.0,1.0,0.0,2.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,
// 61 - big diamond with line
0.0,0.0,0.0,1.0,1.0,0.0,2.0,1.0,0.0,1.0,-1.0,0.0,2.0,-1.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,
// 62 - E6 type
2.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 63 - double double E6 type...
1.0,1.0,0.0,3.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,
// 64 - double fork
0.0,2.0,0.0,0.0,-2.0,0.0,1.0,1.0,0.0,1.0,-1.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,
// 65 - double chimney
0.0,1.0,0.0,1.0,1.0,0.0,2.0,1.0,0.0,0.0,-1.0,0.0,1.0,-1.0,0.0,2.0,-1.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,
// 66 - skew E6 type
3.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 67 - E6 type
2.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 68 - triple double E6 type
1.0,1.0,0.0,4.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,
// 69 - double diamond with line
0.0,0.0,0.0,1.0,1.0,0.0,2.0,1.0,0.0,3.0,1.0,0.0,1.0,-1.0,0.0,2.0,-1.0,0.0,3.0,-1.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,
// 70 - E6 type
2.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 71 - E6 type
2.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,
// 72 - quad double E6 type
1.0,1.0,0.0,5.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0};

int customgram[45]={ 0,2,2,2,2,2,2,2,2,0,2,2,2,2,2,2,2,0,2,2,2,2,2,2,0,2,2,2,2,2,0,2,2,2,2,0,2,2,2,0,2,2,0,2,0};
float custompoints[10][3]={
0.0,0.0,0.0,1.0,0.0,0.0,2.0,0.0,0.0,3.0,0.0,0.0,4.0,0.0,0.0,5.0,0.0,0.0,6.0,0.0,0.0,7.0,0.0,0.0,8.0,0.0,0.0,9.0,0.0,0.0};

int hypgram[72][45]={
// 1- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      4,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,5,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
// 2- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,2,2,2,2,2,2,2,5,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
// 3- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      5,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,5,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
// 4- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,5,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
// 5- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,4,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
// 6- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,4,2,2,2,2,2,2,4,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
// 7- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,4,2,2,2,2,2,2,5,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
// 8- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,5,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
// 9- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,5,2,2,2,2,2,2,5,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//10- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      4,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,5,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//11- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,5,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//12- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      5,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,5,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//13- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,5,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//14- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,2,3,2,2,2,2,2,4,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//15- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      4,2,4,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//16- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      4,2,4,2,2,2,2,2,2,3,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//17- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      4,2,4,2,2,2,2,2,2,4,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//18- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,6,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//19- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,4,2,2,2,2,2,2,6,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//20- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,5,2,2,2,2,2,2,6,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//21- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,6,2,2,2,2,2,2,6,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//22- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,3,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//23- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//24- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//25- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,5,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//26- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,6,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//27- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      4,2,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//28- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      4,2,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//29- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      4,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,6,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//30- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      5,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,6,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//31- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,6,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//32- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,2,2,2,2,2,2,2,6,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//33- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      6,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,6,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//34- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,6,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//35- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,4,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//36- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,4,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//37- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,3,3,2,2,2,2,2,2,3,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//38- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,3,2,2,2,2,2,2,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//39- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,4,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//40- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,4,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//41- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,4,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//42- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,2,2,2,2,2,2,2,3,2,3,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//43- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      4,2,2,2,2,2,2,2,2,3,2,3,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//44- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,4,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//45- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,3,3,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//46- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      4,2,2,3,2,2,2,2,2,3,2,2,2,2,2,2,2,4,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//47- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//48- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,4,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//49- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//50- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,4,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//51- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//52- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,3,2,2,2,2,2,2,4,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//53- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,3,2,2,2,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//54- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,3,2,2,2,2,2,2,2,2,2,2,2,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//55- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//56- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,3,2,2,2,2,4,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//57- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,3,2,2,2,2,2,2,4,2,2,2,2,2,2,2,2,2,3,2,2,2,2,4,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//58- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,3,2,2,2,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
//59- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,4,2,2,2,2,2,2,2,2,2,
//60- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,3,2,2,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,2,2,2,2,2,2,
//61- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,3,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,2,2,2,2,2,2,
//62- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,4,2,2,2,2,2,
//63- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,3,2,2,2,2,2,2,2,2,2,3,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,3,2,2,2,2,2,
//64- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,3,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,3,2,2,2,2,2,
//65- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,3,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,3,2,2,2,2,2,
//66- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,2,3,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,3,2,2,3,2,2,
//67- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,3,2,2,4,2,2,
//68- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,3,2,2,2,2,2,2,2,2,2,2,3,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,3,2,2,3,2,2,
//69- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      3,2,2,3,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,3,2,2,3,2,2,2,2,3,2,2,2,3,2,2,3,2,2,
//70- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,3,2,2,3,2,3,
//71- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,3,2,2,3,2,4,
//72- 1,2,3,4,5,6,7,8,9,2,3,4,5,6,7,8,9,3,4,5,6,7,8,9,4,5,6,7,8,9,5,6,7,8,9,6,7,8,9,7,8,9,8,9,9
      2,2,3,2,2,2,2,2,2,2,2,2,2,2,3,2,2,3,2,2,2,2,2,2,3,2,2,2,2,2,3,2,2,2,2,3,2,2,2,3,2,2,3,2,3 };

typedef struct {
  int rank;    // total rank
  int type;    // type for each component
  int parameter; // used only for type I
  int gram[MAXRANK][MAXRANK];
  float points[MAXRANK][3];
  int vertex[MAXRANK];
  int order[MAXRANK];
  int vconvert[MAXRANK];
  int length; 
  } TYPE;

typedef struct {
 int loop;  // 0=not loop 1=loop
 int label; // a number. x and -x are inverses...
} LABELTYPE;
  
typedef struct {
 int len;           // length of relation...
 LABELTYPE rel[MAXRELATIONLEN];   // actual label relation...
 int tp;            // type of the relation... 
                    // tp=1 means R2c
                    // tp=2 means R2b
                    // tp=3 means R2a
 int eliminate;     // 0=don't 1=do eliminate bit...
} RELATION;         // this is a relation...

typedef struct {
  int subs[MAXRANK];// subset
  int num;          // its number if -1 then empty record...
  int adj[MAXRANK]; // a list of adjacent nodes...
  int expansion;    // 0=not expanded 1=expanded...
  float x,y,z;      // coordinates of the node to display...
  int self;               // number of loops for this node... 
  int len[MAXRANK];       // lengths of the v[a,I]'s
  LABELTYPE lab[MAXRANK]; // labels on the edges...
  int span[MAXRANK];  // whether in spanning tree or not...
} NODE;

RELATION dict[MAXDICT]; // this is the dictionary
                        // that translates labels
                        // whether loop or edge
                        // into "group elements"
RELATION relation[MAXRELATIONS]; // relations in labels...
RELATION argument;    // argument to relation procedures...
char message[255];
int posted=0; 
int relationnum=0;    // current number of relations...
int error;      // =1 if an error occurs for eg. exceed limits...
NODE groupoid[MAXNODES]; // the groupoid object...
NODE relgroupoid[MAXNODES]; // used for relations...
int nodeindex=1; // index to next free entry in relgroupoid...
int freenode=1;     // index to next free entry in groupoid...
int edgelabel=1; // numbers non-loop edges...
int looplabel=1; // numbers loop edges...
int spanning=0;  // 0=do not show spanning tree...
                 // 1=show spanning tree in colour...
int showlabels=0;  // 0=do not show labels...
                   // 1=show labels on edges...
TYPE big; // information for big irreducible component...
float distance=(float)MAXRANK;
float distanceg=(float)MAXRANK/2;
GLuint selectBuf[MAXSELECT];
GLint windW = 450, windH = 350;
GLint windWg = 800, windHg=400;
GLint vp[4]; 
GLint vpg[4];
int win,wing;           // window identifiers...
int subset[MAXRANK];    // selected reducible vertices...
int subset0[MAXRANK];   // provisional hold for subset...
int subset1[MAXRANK];   // provisional hold for subset...
int subset2[MAXRANK];   // provisional hold for subset...
int component[MAXRANK]; // a connected component...
TYPE reducible[MAXRANK];// a reducible type up to MAXRANK irreducible components...
int pivot=0;            // holds vertex number of pivot point...
GLfloat colour[6][3]={ 1.0,0.0,0.0,0.0,1.0,0.0,0.3,0.5,1.0,0.1,0.9,0.9,0.4,0.4,0.2,0.8,0.3,0.7};
float roll=0.0;   
float yaw=0.0;
float pitch=0.0;
float rollg=0.0;
float yawg=0.0;
float pitchg=0.0;
float originx=0.0;
float originy=0.0;
float originxg=0.0;
float originyg=0.0;
int mouseX0,mouseY0,motion; 
int mouseX0g,mouseY0g,motiong;
char oldkey='a';
int currentrank=6;
int subsetlength=0; 
int totallength=0; // length of v[a,I] worked out by getnextsubset()
int vailength=0;

void
erasemessage(void){
message[0]='\0';
posted=0;
}

void
printmessage(void){
char *p;

p=&message[0];
while((*p)!='\0'){
 printf("%c",(*p));
 p++;
 }
}

void
postmessage(char *p){
int i;
i=0;
if (posted==0){
while((*p)!='\0'){
 message[i]=(*p);
 i++;
 p++;
 }
message[i]='\0';
}
posted=1;
}

void printlabel(LABELTYPE l){

if(l.loop==1)printf("L"); else printf("E");
printf("(%d)",l.label);
}

void printrelationsshort(void){
// print the totals of relations...
int i,j;
int a[4];

 i=0;
 a[1]=0;
 a[2]=0;
 a[3]=0;
 while(i<relationnum){
 a[relation[i].tp]++; 
 i++; 
 }
 printf("R2a: %d R2b: %d R2c: %d \n",a[3],a[2],a[1]);
 printf("Total relations found: %d\n",relationnum);
} 

void printrelations(void){
// print the list of relations...
int i,j;
int a[4];

 printf("\nRelations...\n");
 printf("\n");
 i=0;
 a[0]=0;
 a[1]=0;
 a[2]=0;
 a[3]=0;
 while(i<relationnum){
 
 printf("Rel no.: %d Rel len: %d Rel type: R2%c\n",i+1,relation[i].len,(char)('d'-relation[i].tp));
 a[relation[i].tp]++; 
 if(relation[i].eliminate==1)a[0]++;
 for(j=0;j<relation[i].len;j++)printlabel(relation[i].rel[j]);
 printf("\n");
 i++; 
 }
 printf("\n");
 printf("R2a: %d R2b: %d R2c: %d\n",a[3],a[2],a[1]);
 printf("Total relations:   %d\n",relationnum);
 
} 

void
printargument(void){
// print argument relation...
int i,j;
 printf("Argument relation is:\n");
 printf("Rel len: %d Rel type: R2%c\n",argument.len,(char)('d'-argument.tp));
 for(j=0;j<argument.len;j++)printlabel(argument.rel[j]);
 printf("\n");
}

void initializerelations(void){
relationnum=0;
}

void 
invertargument(void){
// argument becomes argument^{-1}
int i,k,lo,la,n;

n=argument.len;
k=(int)(n/2);
// first swap order of terms...
for(i=0;i<k;i++){
 lo=argument.rel[i].loop;
 la=argument.rel[i].label;
 argument.rel[i].loop=argument.rel[n-1-i].loop;
 argument.rel[i].label=argument.rel[n-1-i].label;
 argument.rel[n-1-i].loop=lo;
 argument.rel[n-1-i].label=la;
 }
// then invert each term...
for(i=0;i<n;i++){
 if(argument.rel[i].loop==0){
  argument.rel[i].label=-argument.rel[i].label;
  }
 }
 // done!!
}

int iscyclicallyequal(int p){
// return 1 if argument=relation[p]
// or a cyclic conjugate...
// else return 0
int equal=0,s,j,k,l,n;

if ((argument.len==relation[p].len)&&(argument.tp==relation[p].tp)){   

 // is it cyclically equal or not??
 n=argument.len;
 for(k=0;k<n;k++){
  s=0;
  for(j=0;j<n;j++){
   l=(j+k)%n;
   if((argument.rel[l].loop==relation[p].rel[j].loop)&&(argument.rel[l].label==relation[p].rel[j].label))s++;
   }
  if(s==n)equal=1;  
  }
 }
return equal;
}

int isequal(int p){
// return 1 if argument=relation[p]
// or a cyclic conjugate or an inverse of...
// else return 0
int equal=0,s,j,e1,e2;

e1=iscyclicallyequal(p);
invertargument();
e2=iscyclicallyequal(p);
invertargument();
if ((e1==1)||(e2==1))equal=1;
return equal;
}

int searchrelationfrom(int y){
// searches for the relation argument in relation[y+]
// returns the index or -1 if not found...
int i,j,s,found=-1;

i=y;
while(i<relationnum){
 if (isequal(i)==1){ found=i; break; }
 i++;
 }
 return found;
}

int searchrelation(void){
// searches for the relation argument in relation[..]
// returns the index or -1 if not found...
 return searchrelationfrom(0);
}

void
replacerelation(int p){
int i;
     
     relation[p].len=argument.len;
     relation[p].tp=argument.tp;
     for(i=0;i<argument.len;i++){
     relation[p].rel[i].loop=argument.rel[i].loop;
     relation[p].rel[i].label=argument.rel[i].label;
     }
     
}

void addrelation(void){
// adds argument to list of relations...
// only if it is not in the list already...
int j,found;
 found=searchrelation();
 if ((found==-1)&&(argument.len>0)){
 // add it only if it is not in the list already...
 relation[relationnum].len=argument.len;
 relation[relationnum].tp=argument.tp;
 relation[relationnum].eliminate=argument.eliminate;
 for(j=0;j<argument.len;j++){
 relation[relationnum].rel[j].loop=argument.rel[j].loop;
 relation[relationnum].rel[j].label=argument.rel[j].label;
 }  
 relationnum++;
 if(relationnum>=MAXRELATIONS){ 
    error=1; 
    relationnum=MAXRELATIONS-1;
    postmessage("Error. Too many relations.\n");
    }
}
// do nothing if relation is already there...
}

void
deleterelation(int p){
// deletes relation numbered p
// ie. deletes relation[p]...
int i,j;

if(p<0)p=0;
i=p+1;
while(i<relationnum){
   relation[i-1].len=relation[i].len;
   relation[i-1].tp=relation[i].tp;
   relation[i-1].eliminate=relation[i].eliminate;
   for(j=0;j<relation[i].len;j++){
    relation[i-1].rel[j].loop=relation[i].rel[j].loop;
    relation[i-1].rel[j].label=relation[i].rel[j].label;   
   }  
   i++;  
   }
if(p<relationnum)relationnum--;
}

void delrelation(void){
// deletes the relation argument...
// only if it is in the list already...
int found;
  found=searchrelation();
  if(found!=-1){
  // we must delete this entry then...
  // entry found must be deleted...
  deleterelation(found);
  }
// do nothing if not found...
}

void freereduceargument(void){
// freely reduce argument...
int change,i,j,n,u;

change=1;
while(change==1){
  change=0;
  if(argument.len>1){

  i=0;
  n=argument.len;
  while(i<n){
  
 if(((argument.rel[i%n].loop==1)&&(argument.rel[(i+1)%n].loop==1)&&(argument.rel[i%n].label==argument.rel[(i+1)%n].label))||((argument.rel[i%n].loop==0)&&(argument.rel[(i+1)%n].loop==0)&&(argument.rel[i%n].label==-argument.rel[(i+1)%n].label))){
     // it is either two loops together or two inverse edges together so delete these...
     
     if((i+1)!=n){
     j=i+2;
     while(j<argument.len){
      argument.rel[(j-2)].loop=argument.rel[j].loop;
      argument.rel[(j-2)].label=argument.rel[j].label;
      j++;
      } 
     } else {
     j=1;
     while(j<argument.len){
      argument.rel[(j-1)].loop=argument.rel[j].loop;
      argument.rel[(j-1)].label=argument.rel[j].label;
      j++;
      }
     }
     
     argument.len-=2;       
     change=1;
     break;
     } else i++;
    }
   }
  }
}

void
freereduce(int p){
// freely reduce relation numbered p ie. relation[p]
// into argument...
int i;

// first copy the relation to argument...
for(i=0;i<relation[p].len;i++){
  argument.rel[i].loop=relation[p].rel[i].loop;
  argument.rel[i].label=relation[p].rel[i].label;
 }
argument.len=relation[p].len;
argument.tp=relation[p].tp;
// the relation has been copied, now to reduce it...
freereduceargument();
// we are now finished. the freely reduced relation
// corresponding to relation[p] is now in argument...
}

void
printdict(void){
int i,j;
// prints out the dictionary...

printf("\nDictionary...\n");
printf("\n");
for(i=0;i<MAXDICT;i++){

 if (dict[i].len>0){  
  if (i<=looplabel)printf("L(%d): ",i); else printf("E(%d): ",i-looplabel);
  for(j=0;j<dict[i].len;j++)printlabel(dict[i].rel[j]);
  printf("\n");
  }
 }
printf("\n");
}

void
initializedict(void){
int i;
// initializes the dictionary...
// initially, the dictionary just says that
// label=label
// we map the label (loop,label)
// to    dict[label]           if loop=1
// or to dict[label+looplabel] if loop=0

// this to be empty to begin with
// not really necessary if we are not printing
// the dictionary but necessary to print
// the dictionary...
for(i=0;i<MAXDICT;i++)dict[i].len=-1;

for(i=1;i<looplabel;i++){
 dict[i].len=1;
 dict[i].rel[0].loop=1;
 dict[i].rel[0].label=i;
 }
for(i=1;i<edgelabel;i++){
 dict[i+looplabel].len=1;
 dict[i+looplabel].rel[0].loop=0;
 dict[i+looplabel].rel[0].label=i;
 }
 // that's all to be done for now...
}

void
convertargument(void){
// converts the relation argument
// through the dictionary to the 
// relation argument...
int i,j,k,n,l,kloop,klabel,inverse,change,same,limit;
RELATION temp;

change=1;
limit=0;

while((limit<MAXTRANSLATIONS)&&(change==1)&&(error==0)){

limit++;
// first copy it to temp...
temp.len=argument.len;
temp.tp=argument.tp;
for(i=0;i<temp.len;i++){
 temp.rel[i].loop=argument.rel[i].loop;
 temp.rel[i].label=argument.rel[i].label;
}

argument.len=0;
n=0;
for(i=0;i<temp.len;i++){
  inverse=0;
  kloop=temp.rel[i].loop;
  klabel=temp.rel[i].label;
  if((kloop==0)&&(klabel<0)){ inverse=1; klabel=-klabel; }
  k=klabel+(1-kloop)*looplabel;
  // k is the entry in the dictionary now...
  j=0;
  while((error==0)&&(j<dict[k].len)){
   
   if(inverse==0){
   
   if((n+j)<MAXRELATIONLEN){
   argument.rel[n+j].loop=dict[k].rel[j].loop;
   argument.rel[n+j].label=dict[k].rel[j].label;
   } else { 
            error=1; 
            postmessage("Error. Relation too long to convert...\n"); 
          }
    } else {
   // it is an inverse...
   l=n+dict[k].len-1-j;
   if(l<MAXRELATIONLEN){
   argument.rel[l].loop=dict[k].rel[j].loop;
   if(dict[k].rel[j].loop==0){ 
   argument.rel[l].label=-dict[k].rel[j].label;
   } else {
   argument.rel[l].label=dict[k].rel[j].label;
     }     
    } else { error=1; 
             postmessage("Error. Relation too long to convert...\n");
            }
   }
   j++;
   }
   if(error==0)n+=dict[k].len;
  }
argument.len=n;

change=1;
// compare temp with argument see if changed or not...
if(temp.len==argument.len){
  same=0;
  for(i=0;i<temp.len;i++){
    if((temp.rel[i].loop==argument.rel[i].loop)&&(temp.rel[i].label==argument.rel[i].label))same++;
    }
  if(same==temp.len)change=0;
  }
 }
 if (limit>=MAXTRANSLATIONS){
  error=1;
  postmessage("Error. Maximum number of translations exceeded.\n");
  }
 if(error==0){  freereduceargument(); } else 
 { argument.len=0; }
  // done! now argument has been translated fully and reduced...
}

char converttochar(int t)
{
char r;

r='X';
switch(t){
    case A:
        r='A';
        break;
    case B:
        r='B';
        break;
    case D:
        r='D';
        break;
    case E:
        r='E';
        break;
    case F:
        r='F';
        break;
    case H:
        r='H';
        break;
    case I:
        r='I';
        break;
    case X:
        r='+';
        break;
   }
return r;
}

TYPE recognisecomponent(void){
// returns the TYPE of component subset "component"...
TYPE r;
TYPE *t;
int i,j,k,l,m,n,s,u,found,last;
int three[MAXRANK];
int four[MAXRANK];
int five[MAXRANK];
int other[MAXRANK];
int vertexconvert[MAXRANK];
int onethree=0,onefour=0,onefive=0,twothree=0,onefouronethree=0,threethree=0,onefiveonethree=0;
int vertexthreethree=-1; // number of three three vertex...
int vertexonethree[3]={-1,-1,-1};
int provtype=-1; // provisional type...
int numbered[MAXRANK]; // a tally of the ones we've numbered...
int horder[4]={1,2,3,4}; // for ordering type H

  t=&big;
  for(i=0;i<MAXRANK;i++){r.vertex[i]=0;}
  k=0;
  
  for(i=0;i<t->rank;i++){
   if (component[i]==1){
     r.vertex[i]=1;
     r.points[k][0]=t->points[i][0];
     r.points[k][1]=t->points[i][1];
     r.points[k][2]=t->points[i][2];   
   l=0;
   for(j=0;j<t->rank;j++){
    if (component[j]==1){
    r.gram[k][l]=t->gram[i][j];
    r.gram[l][k]=t->gram[i][j];
    l++;
    } 
   }
   k++;
    } 
  }

  // k is now the rank of the connected component...
  r.rank=k;
  
  // the rank, gram matrix and points fields are all done...
  // now we must do valency counts to determine the type...
  // we now work out the three,four and five valencies...
  for(i=0;i<r.rank;i++){
    three[i]=0;
    four[i]=0;
    five[i]=0;
    other[i]=0;
  }
  for(i=0;i<r.rank;i++){
    for(j=0;j<r.rank;j++)
    {
     if (r.gram[i][j]==3)three[i]++;
     else
     if (r.gram[i][j]==4)four[i]++;
     else
     if (r.gram[i][j]==5)five[i]++;
     else 
     if ((r.gram[i][j]!=2)&&(r.gram[i][j]!=1))other[i]++;
    }
  }
  m=0;
  for (i=0;i<r.rank;i++){
  
  if ((other[i]==0)&&(three[i]==1)&&(four[i]==0)&&(five[i]==0)){vertexonethree[m]=i; m++; onethree++;}
  if ((other[i]==0)&&(three[i]==2)&&(four[i]==0)&&(five[i]==0))twothree++;
  if ((other[i]==0)&&(three[i]==0)&&(four[i]==1)&&(five[i]==0))onefour++;
  if ((other[i]==0)&&(three[i]==1)&&(four[i]==1)&&(five[i]==0))onefouronethree++;
  if ((other[i]==0)&&(three[i]==3)&&(four[i]==0)&&(five[i]==0)){ vertexthreethree=i;threethree++;}
  if ((other[i]==0)&&(three[i]==1)&&(four[i]==0)&&(five[i]==1))onefiveonethree++;
  if ((other[i]==0)&&(three[i]==0)&&(four[i]==0)&&(five[i]==1))onefive++;
  }
  if (vertexthreethree!=-1){
    n=0;
    for (m=0;m<3;m++){
     if (r.gram[vertexthreethree][vertexonethree[m]]==3)n++; }
    if (n>1)provtype=D;
    if (n==1)provtype=E;  
  }
  
  if (r.rank<2){ r.type=A; }
   else
  if (r.rank==2){
    if (r.gram[0][1]==3){
      r.type=A;
      
    } else
    if (r.gram[0][1]==4){
      r.type=B;} else
    if (r.gram[0][1]>4){
      r.parameter=r.gram[0][1];
      r.type=I;} else { r.type=X; }
    } 
   else    
  if ((r.rank==4)&&(onefive==1)&&(onefiveonethree==1)&&(threethree==0)&&(onefouronethree==0)&&(onefour==0)&&(twothree==1)&&(onethree==1)){ r.type=H; }
   else
  if ((r.rank==3)&&(onefive==1)&&(onefiveonethree==1)&&(threethree==0)&&(onefouronethree==0)&&(onefour==0)&&(twothree==0)&&(onethree==1)){ r.type=H; }
   else
  if ((r.rank==4)&&(onefive==0)&&(onefiveonethree==0)&&(threethree==0)&&(onefouronethree==2)&&(onefour==0)&&(twothree==0)&&(onethree==2)){ r.type=F; }
   else
  if ((provtype==D)&&(onefive==0)&&(onefiveonethree==0)&&(threethree==1)&&(onefouronethree==0)&&(onefour==0)&&(twothree==r.rank-4)&&(onethree==3)){ r.type=D; }
   else
  if ((provtype==E)&&(r.rank==8)&&(onefive==0)&&(onefiveonethree==0)&&(threethree==1)&&(onefouronethree==0)&&(onefour==0)&&(twothree==4)&&(onethree==3)){ r.type=E; }
   else
  if ((provtype==E)&&(r.rank==7)&&(onefive==0)&&(onefiveonethree==0)&&(threethree==1)&&(onefouronethree==0)&&(onefour==0)&&(twothree==3)&&(onethree==3)){ r.type=E; }
   else
  if ((provtype==E)&&(r.rank==6)&&(onefive==0)&&(onefiveonethree==0)&&(threethree==1)&&(onefouronethree==0)&&(onefour==0)&&(twothree==2)&&(onethree==3)){ r.type=E; }
   else
  if ((onefive==0)&&(onefiveonethree==0)&&(threethree==0)&&(onefouronethree==1)&&(onefour==1)&&(twothree==r.rank-3)&&(onethree==1)){ r.type=B; }
   else 
  if ((onethree==2)&&(twothree==r.rank-2)){ r.type=A; } 
   else { r.type=X; }

  // now we finally work out the ordering consistent
  // with Franzen's thesis...
  for(i=0;i<MAXRANK;i++)r.order[i]=0;
  
  n=0;
  for(i=0;i<MAXRANK;i++){
   if (r.vertex[i]==1){
     vertexconvert[n]=i;
     r.vconvert[n]=i;
     n++;
   }
  }
  for(i=0;i<MAXRANK;i++)numbered[i]=r.vertex[i];
    
  if (r.type==A){
  
    if (r.rank<2){
      r.order[vertexconvert[0]]=1;    
      r.length=1;
    } else 
    {
    // find a vertex with a one three valency and number
    // it one then the adjacent two and so on...
    for(i=0;i<r.rank;i++){
      if ((three[i]==1)&&(four[i]==0)&&(five[i]==0))break;    
    }
    // i is now the vertex with a one three valency
    // so number it one...
    r.order[vertexconvert[i]]=1;
    numbered[vertexconvert[i]]=0; // we've numbered this one...
    last=i;
    for (n=1;n<r.rank;n++){
     for (i=0;i<r.rank;i++){
     // now number the next ones...
      if ((numbered[vertexconvert[i]]==1)&&(r.gram[last][i]>2)){
        numbered[vertexconvert[i]]=0; r.order[vertexconvert[i]]=n+1;
        last=i;
        break;           
      } 
    }
    }
    r.length=(r.rank*(r.rank+1))/2;
    }
    // done!!
   }
  else
  if (r.type==B){
   
    // find a vertex with a one four valency and number
    // it one then the adjacent two and so on...
    for(i=0;i<r.rank;i++){
      if ((three[i]==0)&&(four[i]==1)&&(five[i]==0))break;    
    }
    // i is now the vertex with a one four valency
    // so number it one...
    r.order[vertexconvert[i]]=1;
    numbered[vertexconvert[i]]=0; // we've numbered this one...
    last=i;
    for (n=1;n<r.rank;n++){
     for (i=0;i<r.rank;i++){
     // now number the next ones...
      if ((numbered[vertexconvert[i]]==1)&&(r.gram[last][i]>2)){
        numbered[vertexconvert[i]]=0; r.order[vertexconvert[i]]=n+1;
        last=i;
        break;           
      } 
    }
    }
    r.length=(r.rank)*(r.rank);
   // done!!
  }
  else
  if (r.type==D){
    // vertexthreethree is the vertex with three three's
    // so we number it 3.
    numbered[vertexconvert[vertexthreethree]]=0;
    r.order[vertexconvert[vertexthreethree]]=3;
    // now we number the fork...
    n=0;
    for (i=0;i<onethree;i++){
    if ((n<2)&&(r.gram[vertexthreethree][vertexonethree[i]]>2)){
    numbered[vertexconvert[vertexonethree[i]]]=0;
    r.order[vertexconvert[vertexonethree[i]]]=n+1;
    n++;
    }
    }
    // now we must number the remaining vertices...
    last=vertexthreethree;
    for (n=3;n<r.rank;n++){
     for (i=0;i<r.rank;i++){
     // now number the next ones...
      if ((numbered[vertexconvert[i]]==1)&&(r.gram[last][i]>2)){
        numbered[vertexconvert[i]]=0; r.order[vertexconvert[i]]=n+1;
        last=i;
        break;           
      } 
    }
    }
    // done!!
    r.length=(r.rank)*(r.rank-1);
  }
  else
  if (r.type==E){
     // vertexthreethree is the vertex with three three's coming out of it...
     // this vertex must be numbered 3
     numbered[vertexconvert[vertexthreethree]]=0;
     r.order[vertexconvert[vertexthreethree]]=3;
     
     
     s=0;
     for(i=0;i<3;i++){
     if(r.gram[vertexthreethree][vertexonethree[i]]==3)
       s=i; }
     
     numbered[vertexconvert[vertexonethree[s]]]=0;
     r.order[vertexconvert[vertexonethree[s]]]=4;
     
     // the one adjacent to three and adjacent to a 
     // onethree we must number two...
     found=0;
     for(i=0;i<r.rank;i++){
     // we test i to see if it is adjacent to the 
     // vertex numbered three and adjacent to a 
     // onethree in the diagram...
      if (r.gram[i][vertexthreethree]==3){
       u=-1;
       for(s=0;s<3;s++){
         if (r.gram[vertexonethree[s]][i]==3)u=s;
          }
       if(u!=-1){ 
        numbered[vertexconvert[i]]=0;
        r.order[vertexconvert[i]]=2;
        found=1;
        break;
          }
         }
       }
     if (found==1){ 
     // i is the vertex numbered 2...
     // so 2,3,4 are numbered...
     // now number 1 which is adjacent to 2...
     for(u=0;u<r.rank;u++){
       if((r.gram[i][u]==3)&&(numbered[vertexconvert[u]]==1))break; 
     } 
     numbered[vertexconvert[u]]=0;
     r.order[vertexconvert[u]]=1; 
     for(u=0;u<r.rank;u++){
       if((r.gram[u][vertexthreethree]==3)&&(numbered[vertexconvert[u]]==1))break;
    }
    numbered[vertexconvert[u]]=0;
    r.order[vertexconvert[u]]=5;
    
    last=u;
    s=6;
    for(j=0;j<r.rank;j++){
    k=-1;
    for(i=0;i<r.rank;i++){
      if((r.gram[i][last]==3)&&(numbered[vertexconvert[i]]==1)){ k=0; break;}
     }
    if(k==0){
    
    numbered[vertexconvert[i]]=0;
    r.order[vertexconvert[i]]=s;
    last=i; s++;}
        }
     
   // done!!
    if (r.rank==6)r.length=36;
    else
    if (r.rank==7)r.length=63;
    else r.length=120;
    } else { 
    
    r.type=X;
   for(i=0;i<r.rank;i++){
     numbered[vertexconvert[i]]=0;
     r.order[vertexconvert[i]]=i+1;   
   }
   r.length=-1; //-1 indicates there is no longest element in infinite type... 
     }
  }
  else
  if (r.type==F){
     for(i=0;i<r.rank;i++){
       numbered[vertexconvert[i]]=0;
       r.order[vertexconvert[i]]=i+1;
     }  
     // done!!
     r.length=24;
  }
  else
  if (r.type==H){
    // we use the array of orderings horder[4]
    for (i=0;i<r.rank;i++){
     numbered[vertexconvert[i]]=0;
     r.order[vertexconvert[i]]=horder[i];
    }
    // done!!
    if(r.rank==3)r.length=15; else r.length=60;
  }
  else
  if (r.type==I){
    for(i=0;i<r.rank;i++){
      numbered[vertexconvert[i]]=0;
      r.order[vertexconvert[i]]=i+1;
    }    
    r.length=r.parameter; 
  }
  else
  if (r.type==X){
   for(i=0;i<r.rank;i++){
     numbered[vertexconvert[i]]=0;
     r.order[vertexconvert[i]]=i+1;   
   }
   r.length=-1; //-1 indicates there is no longest element in infinite type...
  }
  
  return r;
}

void getsubsettype(void){
// fills in the reducible type structure with the
// irreducible components of subset...
int i,j,k,l;
TYPE* t;
int com[MAXRANK];
int sub[MAXRANK];
int s,z;

   t=&big;
   for(i=0;i<MAXRANK;i++){
     reducible[i].rank=0;
   }
  // we have initialized the components to have all 0 rank...
  // we now work out the irreducible components of subset...
  // get a copy of subset...
  for(i=0;i<t->rank;i++)sub[i]=subset[i];
  s=0;
  for(i=0;i<t->rank;i++){
  // i is the vertex number we are working the component for...
  if (sub[i]==1){ 
  // work out the component that i belongs to...
  
  for(j=0;j<MAXRANK;j++)com[j]=0; 
  com[i]=1;
  
  for(z=0;z<t->rank;z++){
  for(k=0;k<t->rank;k++){
   if (com[k]==1){
    for(l=0;l<t->rank;l++){
     if ((sub[l]==1)&&(t->gram[k][l]!=1)&&(t->gram[k][l]!=2)){
       com[l]=1;
       sub[l]=0;
     }
    }
   } 
  }
  }
  // now com is the component that i belongs to...
  for(j=0;j<MAXRANK;j++)component[j]=com[j];
  // copy it to component...
  reducible[s]=recognisecomponent();
  s++;
   }
  } 
  subsetlength=totallength;
  totallength=0;
  for(j=0;j<MAXRANK;j++){
    if((reducible[j].rank>0)&&(reducible[j].type==X)){
     totallength=-1;
    }
  } 
  if(totallength==0){
  
  for(j=0;j<MAXRANK;j++){
   if(reducible[j].rank>0)totallength+=reducible[j].length;  
  }
  
  }  
  
}

void printsubsettype(void){
// print out the type of the subset...
int i,j;
  getsubsettype();
  printf("Type selected: ");
  
  for(i=0;i<MAXRANK;i++){
    if (reducible[i].rank!=0){
      printf("+%c%d",converttochar(reducible[i].type),reducible[i].rank);   
      if (reducible[i].type==I)printf("(%d)",reducible[i].parameter);
     //printf("\n");
     //for(j=0;j<MAXRANK;j++){
     // printf("%d ",reducible[i].order[j]);
     //}      
    }
  }
  printf(" Longest element: %d\n",totallength);
}

int correcttype(TYPE *t){
// correct the type according to the rules:
// if Type A: any rank
// if Type B: rank>=2, B1=A1
// if Type D: rank>=4, D1=A1, D2=A2, D3=A3
// if Type E: rank 6,7,8, any other is converted to type A
// if Type F: rank 4 only, any other converted to type A
// if Type H: rank 3 or 4 only, any other converted to type A
// if Type I: rank 2 only, any other converted to type A
if (t->rank<1)t->rank=1;
if (t->rank>MAXRANK)t->rank=MAXRANK;
 
switch(t->type){
    case A:
       break;
    case B:
       if (t->rank==1)t->type=A;
       break;
    case D:
       if (t->rank<4)t->type=A;
       break;
    case E:
       if ((t->rank<6)||(t->rank>8))t->type=A;
       break;
    case F:
       if (t->rank!=4)t->type=A;
       break;
    case H:
       if ((t->rank<3)||(t->rank>4))t->type=A;
       break;
    case I:
       if (t->rank!=2)t->type=A;
       if (t->parameter<=4)t->parameter=5;
       break;
    case X:
       if (t->rank==1)t->type=A;
       
       if (t->parameter==A){
       
        if(t->rank==2)t->parameter=0;
       
       }
       else
       if ((t->rank==2)&&(t->parameter>2)){t->type=I; break;}
       else
       if (t->parameter==B){
         if (t->rank<4)t->rank=4;
       }
       else
       if (t->parameter==C){
            
       }
       else
       if (t->parameter==D){
         if (t->rank<5)t->rank=5;
       
       
       }
       else
       if (t->parameter==E){
       
         if((t->rank<7)||(t->rank>9))t->rank=7;
       
       }
       else
       if (t->parameter==F){
       
         t->rank=5;
       
       }
       else
       if (t->parameter==G){
       
         t->rank=3;
       
       } else
       if ((t->parameter>=HYPBASE)&&(t->parameter<CUSTOM)){
        t->rank=hyprank[t->parameter-HYPBASE];
       }
       else
       if (t->parameter==CUSTOM){
        t->rank=CUSTOMRANK;
       }
       else
       if (t->parameter==0){
        t->rank=2;
       }
       break;
  }
return 0;
}

int inittype(TYPE *t)
{
// initializes a TYPE structure for irreducible component...
// given the rank and type it fills in the gram matrix
// and the point coordinates...
int i,j,k,n;
float a;
  correcttype(t);
  correcttype(t);
  // yes do it twice, this is not an error...
  for(i=0;i<t->rank;i++){t->vertex[i]=1; t->order[i]=i;}
  printf("Irreducible Component of Type ");
  printf("%c%d",converttochar(t->type),t->rank);
  if (t->type==I)printf("(%d)\n",t->parameter); else printf("\n");

  
  for(i=0;i<t->rank;i++){
    for(j=0;j<t->rank;j++)
      if(i==j)t->gram[i][j]=1; else t->gram[i][j]=2; } 

  switch(t->type){
     case A:
     // type A
     // 0 is joined to 1 is joined to 2 ... with 3's
        for(i=0;i<t->rank-1;i++){
          t->gram[i][i+1]=3;
          t->gram[i+1][i]=3;
        }
        a=((float)MAXRANK)/((float)t->rank);
        for(i=0;i<t->rank;i++){
          t->points[i][0]=(float)a*i;
          t->points[i][1]=0.0;
          t->points[i][2]=0.0;
        }
        break;
      case B:
      // type B
      // like type A but with a 4 linking 0 to 1...
        for(i=0;i<t->rank-1;i++){
          t->gram[i][i+1]=3;
          t->gram[i+1][i]=3;
        }
        t->gram[0][1]=4;
        t->gram[1][0]=4;
        a=((float)MAXRANK)/((float)t->rank);
        for(i=0;i<t->rank;i++){
          t->points[i][0]=(float)a*i;
          t->points[i][1]=0.0;
          t->points[i][2]=0.0;
        }
        break;
      case D:
      // type D
      // assuming rank is at least 4...
        for (i=1;i<t->rank-1;i++){
          t->gram[i][i+1]=3;
          t->gram[i+1][i]=3;
        }
        t->gram[0][2]=3;
        t->gram[2][0]=3;
        a=((float)MAXRANK)/((float)(t->rank-1));
        for(i=2;i<t->rank;i++){
          t->points[i][0]=(float)a*(i-1);
          t->points[i][1]=0.0;
          t->points[i][2]=0.0;
        }
        t->points[0][0]=0.0;
        t->points[0][1]=2.0;
        t->points[0][2]=0.0;
        t->points[1][0]=0.0;
        t->points[1][1]=-2.0;
        t->points[1][2]=0.0;       
        break;
      case E:
      // type E
      // rank must be 6,7 or 8.
        for(i=1;i<8;i++){
          t->gram[i][i+1]=3;
          t->gram[i+1][i]=3;
        }
        t->gram[0][3]=3;
        t->gram[3][0]=3;
        a=((float)MAXRANK)/((float)t->rank);
        for(i=1;i<t->rank;i++){
          t->points[i][0]=(float)a*i;
          t->points[i][1]=0.0;
          t->points[i][2]=0.0;
        }
        t->points[0][0]=t->points[3][0];
        t->points[0][1]=2.0;
        t->points[0][2]=0.0;
        break;
      case F:
      // type F
      // rank must be 4
      t->gram[0][1]=3;
      t->gram[1][0]=3;
      t->gram[1][2]=4;
      t->gram[2][1]=4;
      t->gram[2][3]=3;
      t->gram[3][2]=3;
      a=((float)MAXRANK)/((float)t->rank);
      for(i=0;i<t->rank;i++){
        t->points[i][0]=(float)a*i;
        t->points[i][1]=0.0;
        t->points[i][2]=0.0;
      }
        break;
      case H:
      // type H 
      // rank must be 3 or 4...
      t->gram[0][1]=5;
      t->gram[1][0]=5;
      t->gram[1][2]=3;
      t->gram[2][1]=3;
      t->gram[2][3]=3;
      t->gram[3][2]=3;
      a=((float)MAXRANK)/((float)t->rank);
      for(i=0;i<t->rank;i++){
        t->points[i][0]=(float)a*i;
        t->points[i][1]=0.0;
        t->points[i][2]=0.0;
      }
      break;
      case I:
      // type I
      // rank must be 2...
      t->gram[0][1]=t->parameter;
      t->gram[1][0]=t->parameter;
      a=((float)MAXRANK)/((float)t->rank);
      for(i=0;i<t->rank;i++){
        t->points[i][0]=(float)a*i;
        t->points[i][1]=0.0;
        t->points[i][2]=0.0;
      }
      break;  
      case X:
      // infinite type ...
      if(t->parameter==0){
       t->gram[0][1]=0;
       t->gram[1][0]=0;
       t->rank=2;
       a=((float)MAXRANK)/((float)t->rank);
       t->points[0][0]=0.0;
       t->points[0][1]=0.0;
       t->points[0][2]=0.0;
       t->points[1][0]=a;
       t->points[1][1]=0.0;
       t->points[1][2]=0.0;
      }
      else      
      if (t->parameter==A){
      
      for(i=0;i<t->rank-1;i++){
          t->gram[i][i+1]=3;
          t->gram[i+1][i]=3;
        }
        t->gram[t->rank-1][0]=3;
        t->gram[0][t->rank-1]=3;
        
        a=((float)MAXRANK)/((float)t->rank);
        for(i=1;i<t->rank;i++){
          t->points[i][0]=(float)a*i;
          t->points[i][1]=0.0;
          t->points[i][2]=0.0;
        }
        t->points[0][0]=a*((float)t->rank/2);
        t->points[0][1]=2.0;
        t->points[0][2]=0.0;
      }
      else
      if (t->parameter==B){
      
      if (t->rank>3){
      for (i=1;i<t->rank-1;i++){
          t->gram[i][i+1]=3;
          t->gram[i+1][i]=3;
        }
        t->gram[0][2]=3;
        t->gram[2][0]=3;
        t->gram[t->rank-2][t->rank-1]=4;
        t->gram[t->rank-1][t->rank-2]=4;
        a=((float)MAXRANK)/((float)(t->rank-1));
        for(i=2;i<t->rank;i++){
          t->points[i][0]=(float)a*(i-1);
          t->points[i][1]=0.0;
          t->points[i][2]=0.0;
        }
        t->points[0][0]=0.0;
        t->points[0][1]=2.0;
        t->points[0][2]=0.0;
        t->points[1][0]=0.0;
        t->points[1][1]=-2.0;
        t->points[1][2]=0.0;
        } else {
        t->rank=3;
        t->gram[0][1]=4; 
        t->gram[1][0]=4;
        t->gram[1][2]=4;
        t->gram[2][1]=4;
        a=((float)MAXRANK)/((float)t->rank);
        for(i=0;i<t->rank;i++){
         t->points[i][0]=(float)a*i;
         t->points[i][1]=0.0;
         t->points[i][2]=0.0;
         }
        }
      } 
      else
      if (t->parameter==C){
      
       for(i=0;i<t->rank-1;i++){
          t->gram[i][i+1]=3;
          t->gram[i+1][i]=3;
        }
        t->gram[0][1]=4;
        t->gram[1][0]=4;
        t->gram[t->rank-2][t->rank-1]=4;
        t->gram[t->rank-1][t->rank-2]=4;
        a=((float)MAXRANK)/((float)t->rank);
        for(i=0;i<t->rank;i++){
          t->points[i][0]=(float)a*i;
          t->points[i][1]=0.0;
          t->points[i][2]=0.0;
        }
      }
      else
      if (t->parameter==D){
      
       for (i=1;i<t->rank-2;i++){
          t->gram[i][i+1]=3;
          t->gram[i+1][i]=3;
        }
        t->gram[0][2]=3;
        t->gram[2][0]=3;
        t->gram[t->rank-3][t->rank-1]=3;
        t->gram[t->rank-1][t->rank-3]=3;
        a=((float)MAXRANK)/((float)(t->rank-1));
        for(i=2;i<t->rank-2;i++){
          t->points[i][0]=(float)a*(i-1);
          t->points[i][1]=0.0;
          t->points[i][2]=0.0;
        }
        t->points[0][0]=0.0;
        t->points[0][1]=2.0;
        t->points[0][2]=0.0;
        t->points[1][0]=0.0;
        t->points[1][1]=-2.0;
        t->points[1][2]=0.0;  
        t->points[t->rank-1][0]=(float)a*(t->rank-3);
        t->points[t->rank-1][1]=2.0;
        t->points[t->rank-1][2]=0.0;
        t->points[t->rank-2][0]=(float)a*(t->rank-3);
        t->points[t->rank-2][1]=-2.0;
        t->points[t->rank-2][2]=0.0;  
      }
      else
      if (t->parameter==E){
      
      a=((float)MAXRANK)/((float)t->rank);
      for(i=1;i<t->rank;i++){
       t->points[i][0]=(float)a*i;
       t->points[i][1]=0.0;
       t->points[i][2]=0.0;
      }
      
      if (t->rank==7){
      t->points[0][0]=t->points[4][0];
      t->points[0][1]=t->points[4][1]+4.0;
      t->points[0][2]=0.0;
      t->points[1][0]=t->points[4][0];
      t->points[1][1]=t->points[4][1]+2.0;
      t->points[1][2]=0.0;
      for(i=2;i<t->rank-1;i++){
        t->gram[i][i+1]=3;
        t->gram[i+1][i]=3;
      }
      t->gram[0][1]=3;
      t->gram[1][0]=3;
      t->gram[1][4]=3;
      t->gram[4][1]=3;
     
      }
      else
      if (t->rank==8){
      t->points[0][0]=t->points[4][0];
      t->points[0][1]=t->points[4][1]+2.0;
      t->points[0][2]=0.0;
      for(i=1;i<t->rank-1;i++){
        t->gram[i][i+1]=3;
        t->gram[i+1][i]=3;
      }
      t->gram[0][4]=3;
      t->gram[4][0]=3;
        }
      else {
      t->points[0][0]=t->points[3][0];
      t->points[0][1]=t->points[3][1]+2.0;
      t->points[0][2]=0.0;
      for(i=1;i<t->rank-1;i++){
        t->gram[i][i+1]=3;
        t->gram[i+1][i]=3;
      }
      t->gram[0][3]=3;
      t->gram[3][0]=3;
       }
      }
      else
      if (t->parameter==F){
        for(i=0;i<t->rank-1;i++){
          t->gram[i][i+1]=3;
          t->gram[i+1][i]=3;
        }
        t->gram[2][3]=4;
        t->gram[3][2]=4;
        a=((float)MAXRANK)/((float)t->rank);
        for(i=0;i<t->rank;i++){
          t->points[i][0]=(float)a*i;
          t->points[i][1]=0.0;
          t->points[i][2]=0.0;
        }
      } 
      else
      if (t->parameter==G){
      
      for(i=0;i<t->rank-1;i++){
          t->gram[i][i+1]=3;
          t->gram[i+1][i]=3;
        }
        t->gram[1][2]=6;
        t->gram[2][1]=6;
        a=((float)MAXRANK)/((float)t->rank);
        for(i=0;i<t->rank;i++){
          t->points[i][0]=(float)a*i;
          t->points[i][1]=0.0;
          t->points[i][2]=0.0;
        }
      }
      else
      if ((t->parameter>=HYPBASE)&&(t->parameter<CUSTOM)){
      
       a=((float)MAXRANK)/((float)t->rank);
       for(i=0;i<t->rank;i++){
        t->points[i][0]=(float)a*hyppoints[t->parameter-HYPBASE][i][0];
        t->points[i][1]=(float)2*hyppoints[t->parameter-HYPBASE][i][1];
        t->points[i][2]=hyppoints[t->parameter-HYPBASE][i][2];
       }
       n=0;
       // now the gram matrix...
       for(i=0;i<10;i++){
        for(j=i+1;j<10;j++){
          t->gram[i][j]=hypgram[t->parameter-HYPBASE][n];
          t->gram[j][i]=t->gram[i][j];
          n++;                  
        }
       }    
      }
      else
      if(t->parameter==CUSTOM){
      
      a=((float)MAXRANK)/((float)t->rank);
      for(i=0;i<t->rank;i++){
      t->points[i][0]=(float)a*custompoints[i][0];
      t->points[i][1]=(float)2*custompoints[i][1];
      t->points[i][2]=custompoints[i][2];
      }
      n=0;
      for(i=0;i<10;i++){
       for(j=i+1;j<10;j++){
        t->gram[i][j]=customgram[n];
        t->gram[j][i]=t->gram[i][j];
        n++; 
        }
       }
      }
      break;
  }
  return 0;
}

void
stroke_output(GLfloat x, GLfloat y, GLfloat z,char *format,...)
{
  va_list args;
  char buffer[200], *p;

  va_start(args, format);
  vsprintf(buffer, format, args);
  va_end(args);
  glPushMatrix();
  glTranslatef(x, y, z);
  glScalef(0.005, 0.005, 0.005);
  for (p = buffer; *p; p++)
    glutStrokeCharacter(GLUT_STROKE_ROMAN, *p);
  glPopMatrix();
}

void
stroke_outputsmall(GLfloat x, GLfloat y, GLfloat z,char *format,...)
{
  va_list args;
  char buffer[200], *p;

  va_start(args, format);
  vsprintf(buffer, format, args);
  va_end(args);
  glPushMatrix();
  glTranslatef(x, y, z);
  glScalef(0.0035, 0.0035, 0.0035);
  for (p = buffer; *p; p++)
    glutStrokeCharacter(GLUT_STROKE_ROMAN, *p);
  glPopMatrix();
}

void myinit(void)
{    
    glClearColor (0.0, 0.0, 0.0, 0.0);
    glEnable(GL_DEPTH_TEST);
    glDepthFunc(GL_LESS);
}

void
Render(GLenum mode)
{
  GLint i,j,s,n,letter;
  float x,y,z;
  char buffer[MAXRANK*10];
  char auxbuf[20];
  char *mybuffer;
  TYPE* t;
  
  t=&big;
    
    for(i=0;i<t->rank;i++){
      if (mode == GL_SELECT) {
      glLoadName(i);
      }
    if (subset[i]==1)
     glColor3f(colour[0][0],colour[0][1],colour[0][2]); 
    else  
     glColor3f(colour[1][0],colour[1][1],colour[1][2]);
     
     
    glPushMatrix();
    glTranslatef(t->points[i][0],t->points[i][1],t->points[i][2]);
    glNewList(i, GL_COMPILE);
    glutSolidCube(0.2);
    glEndList();
    glCallList(i);
    glPopMatrix();   
    
    glColor3f(colour[3][0],colour[3][1],colour[3][2]);
    x=t->points[i][0];
    y=t->points[i][1];
    z=t->points[i][2];
    if (subset[i]==1){
     // find the reducible component of which this vertex
     // is a part...
     s=0;
     for(n=0;n<MAXRANK;n++){
       if (reducible[n].rank>0){
         if (reducible[n].vertex[i]==1){
          s=n;
          break;
          }              
       }
     }
    // s is the reducible component of which this
    // vertex is a part...
    itoa(reducible[s].order[i],&mybuffer[0],10);
    stroke_output(x,y-1.0,z, mybuffer); 
      }
    }
     glColor3f(colour[1][0],colour[1][1],colour[1][2]);
    if (mode==GL_SELECT){
      glLoadName(t->rank);
    }
    for(i=0;i<t->rank;i++){
      for(j=i+1;j<t->rank;j++){
        if ((t->gram[i][j]<1)||(t->gram[i][j]>2)){
          if ((subset[i]==1)&&(subset[j]==1))
            glColor3f(colour[0][0],colour[0][1],colour[0][2]);
          else
            glColor3f(colour[1][0],colour[1][1],colour[1][2]);
          glBegin(GL_LINES);
          glVertex3f(t->points[i][0],t->points[i][1],t->points[i][2]);
          glVertex3f(t->points[j][0],t->points[j][1],t->points[j][2]);
          glEnd();        
        }
        if (t->gram[i][j]>3){
         glColor3f(colour[2][0],colour[2][1],colour[2][2]);
         x=t->points[i][0]+t->points[j][0];
         y=t->points[i][1]+t->points[j][1];
         z=t->points[i][2]+t->points[j][2];        
         itoa(t->gram[i][j],&mybuffer[0],10);
         stroke_output(x/2,y/2+0.1,z/2, mybuffer); 
        } else
       if (t->gram[i][j]<1){
         glColor3f(colour[2][0],colour[2][1],colour[2][2]);
         x=t->points[i][0]+t->points[j][0];
         y=t->points[i][1]+t->points[j][1];
         z=t->points[i][2]+t->points[j][2];        
         stroke_output(x/2,y/2+0.1,z/2, "INF"); 
        }
      }
    }
    letter=0;
    buffer[letter]='\0';
    for(i=0;i<MAXRANK;i++){
    if (reducible[i].rank>0){
      buffer[letter]=converttochar(reducible[i].type);
      letter++;
      itoa(reducible[i].rank,&auxbuf[0],10);
      n=0;
      while (auxbuf[n]!='\0'){ buffer[letter]=auxbuf[n];
       letter++; n++; 
      }
      if (reducible[i].type==I){
        buffer[letter]='(';
        letter++;
        itoa(reducible[i].parameter,&auxbuf[0],10);
        n=0;
        while (auxbuf[n]!='\0'){ buffer[letter]=auxbuf[n]; 
         letter++; n++;
       }      
       buffer[letter]=')';
       letter++;
      }
      buffer[letter]='+';
      letter++;
    }
   }
   buffer[--letter]='\0';
   glColor3f(colour[2][0],colour[2][1],colour[2][2]);
   stroke_output(t->points[0][0],t->points[0][1]-10.0,0.0,&buffer[0]); 
   glutSwapBuffers();  
}

void getnextsubset(void){
int i,j,k,l,s,source,target;

// reducible has the type of the union of pivot with
// the selected subset...
    // first clear the subset result
    // subset0 holds the source vertices...
    // reducible holds all the information needed...
    for(i=0;i<MAXRANK;i++)subset[i]=0;
    
    for(i=0;i<MAXRANK;i++){
    
    if(subset0[i]==1){
     // now find the reducible component of which this
     // vertex i is a part...
     s=0;
     for(j=0;j<MAXRANK;j++){
      if(reducible[j].rank>0){
       // look in this non empty component...
        if(reducible[j].vertex[i]==1){ s=j; break;}
       }
      }
      // now s is the reducible component of which
      // this vertex i is a part...
      source=reducible[s].order[i];
      
      // source can never be 0, if 0 then error
      if (source==0){
      
      //printf("Error!\n"); 
      //for(j=0;j<MAXRANK;j++){
      //printf("%d ",reducible[s].order[j]);
      //}
      //printf("         %d\n",i);
      postmessage("Error. Source cannot be 0.\n");
      error=1;
      }
      
      if (reducible[s].type==A){
        // type A swaps
        target=reducible[s].rank+1-source;
      }
      else
      if (reducible[s].type==B){
        target=source;
        // type B is central      
      }
      else
      if (reducible[s].type==D){
        // type D depends on whether even or odd...
        k=(int)(reducible[s].rank/2);
        if ((2*k)==reducible[s].rank){
          // case when rank is even...
          target=source;
          // when even it is central...        
        }
        else {
         // case when rank is odd...
         if (source<3){ target=3-source; }
          else { target=source; }
        }
      }
      else
      if (reducible[s].type==E){
       // depending on the rank...
       // central when rank is 7 or 8
       if (reducible[s].rank==6){
         if ((source==3)||(source==4)){ target=source; }
         else { target=7-source; }
       }
       else { target=source; }
      }
      else
      if (reducible[s].type==F){
       // F4 is central...
       target=source;
      }
      else
      if (reducible[s].type==H){
        // both H's are central...
        target=source;
      }
      else
      if (reducible[s].type==I){
        // depends whether parameter is even or odd...
        k=(int)(reducible[s].parameter/2);
        if ((2*k)==reducible[s].parameter){
         // central when even...
         target=source;
        }
        else { 
         target=3-source;
         }
      }
      else
      if (reducible[s].type==X){ target=source; }
      
      // now target is the vertex we must put in subset...
      for(k=0;k<MAXRANK;k++){
       if(reducible[s].order[k]==target)subset[k]=1;
      }
      // done!! next vertex...
     }
    }
}

void initializegroupoid(void){
// initializes the groupoid object...
// it sets the first node to be the current
// subset...
int i;
  for (i=1;i<MAXNODES;i++){
   groupoid[i].num=-1;
  }
  groupoid[0].num=1;
  groupoid[0].expansion=0;
  groupoid[0].x=0.0;
  groupoid[0].y=0.0;
  groupoid[0].z=0.0;
  groupoid[0].self=0;
  for(i=0;i<MAXRANK;i++){
   groupoid[0].subs[i]=0;
   groupoid[0].adj[i]=0;
   groupoid[0].len[i]=0;
  }
  for(i=0;i<MAXRANK;i++){
   groupoid[0].subs[i]=subset[i];
  }
  freenode=1;
}
void printgroupoid(void)
{
// prints out the groupoid information...
int i,j;
printf("Groupoid Structure...\n");
for(i=0;i<freenode;i++){
  printf("  Node number: %d \n",groupoid[i].num);
  printf("   Subset: { ");
  for(j=0;j<big.rank;j++)if(groupoid[i].subs[j]==1)printf(" %d ",j);
  printf("}\n");
  printf("   Self: %d \n",groupoid[i].self);
  printf("   Adjacent: { ");
  for(j=0;j<MAXRANK;j++)if(groupoid[i].adj[j]!=0)printf(" %d ",groupoid[i].adj[j]);
  printf("}\n   Lengths:  { ");
  for(j=0;j<MAXRANK;j++)if(groupoid[i].adj[j]!=0)printf(" %d ",groupoid[i].len[j]);
  printf("}\n   Labels :  { ");
  for(j=0;j<MAXRANK;j++)if(groupoid[i].adj[j]!=0){printlabel(groupoid[i].lab[j]); printf(" ");}
  printf("}\n   Span   :  { ");
  for(j=0;j<MAXRANK;j++)if(groupoid[i].adj[j]!=0)printf(" %d ",groupoid[i].span[j]);
  printf("}\n\n");
 }
}
 
void printrelgroupoid(void)
{
// prints out the groupoid information...
int i,j;
printf("Relgroupoid Structure...\n");
for(i=0;i<nodeindex;i++){
  printf("  Node number: %d\n",relgroupoid[i].num);
  printf("   Subset: { ");
  for(j=0;j<big.rank;j++)if(relgroupoid[i].subs[j]==1)printf(" %d ",j);
  printf("}\n");
  printf("   Self: %d\n",relgroupoid[i].self);
  printf("   Adjacent: { ");
  for(j=0;j<MAXRANK;j++)if(relgroupoid[i].adj[j]!=0)printf(" %d ",relgroupoid[i].adj[j]);
  printf("}\n   Lengths:  { ");
  for(j=0;j<MAXRANK;j++)if(relgroupoid[i].adj[j]!=0)printf(" %d ",relgroupoid[i].len[j]);
  printf("}\n   Labels :  { ");
  for(j=0;j<MAXRANK;j++)if(relgroupoid[i].adj[j]!=0){printlabel(relgroupoid[i].lab[j]); printf(" ");}
  printf("}\n\n"); 
 }
 
}
 
void fillingroupoid(void){
// fills in the current groupoid structure...
int i,j,k,l,s,found,vai;
int infinite;
int complement[MAXRANK];
float line,linelength;

// first we initialize the groupoid and save
// the subset information and initialize the
// edgelabel and looplabel variables to be 1..
  error=0;
  erasemessage();
  edgelabel=1;
  looplabel=1;
  for(i=0;i<MAXRANK;i++)subset1[i]=subset[i];
  
  initializegroupoid();
  srand(1.0);
  
  for(i=0;i<MAXNODES;i++){
  
  if ((groupoid[i].num>0)&&(groupoid[i].expansion==0)){
   
  // we work out the complement of this subset...
  for(j=0;j<big.rank;j++)complement[j]=1-groupoid[i].subs[j];
  for(j=big.rank;j<MAXRANK;j++)complement[j]=0;
  linelength=0.0;
  for(j=0;j<big.rank;j++)if(complement[j]==1)linelength+=1.0;
  line=0.0; 
  for(k=0;k<big.rank;k++){
    
    if (complement[k]==1){
     pivot=k; 
     for(j=0;j<MAXRANK;j++){ 
             subset[j]=groupoid[i].subs[j];
             subset0[j]=subset[j]; }       
     subset[pivot]=1;  // add the pivot to subset...
     infinite=0;      // infinite flag tell us whether the union is finite or not
     getsubsettype();
    // now reducible has the types of the irreducible
    // components of the union of selected subset and
    // the pivot vertex...
     for (j=0;j<MAXRANK;j++){
     
     if (reducible[j].rank>0){
      if (reducible[j].type==X)infinite=1;
      }     
     }
          
     getnextsubset();
     getsubsettype();
     if(infinite!=1)vai=subsetlength-totallength;
     else vai=-1;
    // work out the appropriate transformation...
    // now subset is the next subset...
    // so we must add it to the list...
    // only if it is not in the list already!!
    //printf("Searching for subset... {");
    //for(j=0;j<big.rank;j++)if(subset[j]==1)printf(" %d ",j);
    //printf("} in: \n");
    //printgroupoid(); 
    found=-1;
    for(j=0;j<freenode;j++){
      s=0;
      for(l=0;l<MAXRANK;l++){ 
       if (groupoid[j].subs[l]==subset[l])s++;
      }
      if (s==MAXRANK)found=j;
    }
    if (found==-1){ j=freenode; } else j=found;
    
    // j is now the number of the node, or 
    // j=freenode if it was not found in the list...
    if ((!infinite)&&(j==freenode)&&(freenode<MAXNODES)){
    // we must add the node to the list...
      groupoid[j].num=freenode+1;
      for(l=0;l<MAXRANK;l++){
      groupoid[j].subs[l]=subset[l]; 
      groupoid[j].adj[l]=0;
      groupoid[j].len[l]=0;
      } 
      groupoid[j].expansion=0;
      groupoid[j].self=0;
      groupoid[j].y=groupoid[i].y+3.0*line;
      groupoid[j].x=groupoid[i].x+2.0+((float)rand()/RAND_MAX)*2.0;
      groupoid[j].z=line*2.0-(groupoid[i].x);
      line+=1.0;
      freenode++;  
    }    
    
    if (freenode>=MAXNODES){ error=1; 
     postmessage("Error. Too high number of nodes.\n");
    }
    
    if (!infinite){
    groupoid[i].adj[pivot]=groupoid[j].num;
    groupoid[i].len[pivot]=vai;
    }
    }
   }     
   groupoid[i].expansion=1; // has now been expanded...
   }
  }

for(i=0;i<freenode;i++){
  s=0;
    for(j=0;j<MAXRANK;j++){
      if(groupoid[i].adj[j]==groupoid[i].num)s++;
    }
    groupoid[i].self=s;
  }

// the groupoid has been built.
// all except the labels and spanning tree 
// which we must now process...
// so clear the expansion bit again...
for(i=0;i<freenode;i++){
 groupoid[i].expansion=0;
}

for(i=0;i<freenode;i++){
 if(groupoid[i].expansion==0){
  for(j=0;j<MAXRANK;j++){
   
   if (groupoid[i].adj[j]!=0){
    // there are two cases...
    // either it is a loop or an edge that is not a loop...
    if (groupoid[i].adj[j]==groupoid[i].num){
    // it is a loop...
    // so number it with a loop...
    groupoid[i].lab[j].loop=1;
    groupoid[i].lab[j].label=looplabel;
    looplabel++;    
    }
    else
    {
    // it is not a loop...
    // see if it is an inverse of something already
    // labelled or not...
    k=groupoid[i].adj[j];
    if(groupoid[k-1].expansion==1){
    // it is an inverse
    // we are at groupoid i so num is i+1
    // it is an inverse to num k-1
    s=-1;
    for(l=0;l<MAXRANK;l++)if(groupoid[k-1].adj[l]==i+1)s=l;
    // s cannot be -1 if everything's ok...
    if(s==-1){
    //printf("CRITICAL ERROR. s cannot be -1 here...\n"); 
    // this can only happen if freenode has exceeded
    // the capacity...
    postmessage("Error. Too high number of nodes in groupoid.\n");
    error=1; 
    }
    if(s>=0){
     groupoid[i].lab[j].loop=groupoid[k-1].lab[s].loop;
     groupoid[i].lab[j].label=-groupoid[k-1].lab[s].label;
     }
    }
    else
    {
     // it is a new label
     groupoid[i].lab[j].loop=0;
     groupoid[i].lab[j].label=edgelabel;
     edgelabel++;
     }    
    } 
   }
   // done!
  }
  groupoid[i].expansion=1;
 }
}
// now the edges have been numbered...
// there is an identity that must always hold,
// namely, looplabel+2*edgelabel=valency*numberofnodes...
// (assuming of course constant valency at each node...)

// we now work out a spanning tree (it is not unique)
// so clear expansion which here represents that
// the node is in the tree...
// among the positive edge labels...
// we include those that are between two things such
// that not both are "expanded"...

for(i=0;i<freenode;i++){
 groupoid[i].expansion=0;
}

for(i=0;i<freenode;i++){

 for(j=0;j<MAXRANK;j++){
 
 if(groupoid[i].adj[j]!=0){
 
  if((groupoid[i].lab[j].loop==0)&&(groupoid[i].lab[j].label>0)){
    
    // it is a positive edge...
    k=groupoid[i].adj[j]-1;
    if((groupoid[k].expansion==1)&&(groupoid[i].expansion==1)){
    // it is a positive edge between two things in tree
    // so don't include it or its inverse...
    groupoid[i].span[j]=0;
    for(l=0;l<MAXRANK;l++)if(groupoid[k].adj[l]==groupoid[i].num)break;
    groupoid[k].span[l]=0;    
    
    } else {    
    // include it!! and its inverse...
    groupoid[i].span[j]=1;
    for(l=0;l<MAXRANK;l++)if(groupoid[k].adj[l]==groupoid[i].num)break;
    groupoid[k].span[l]=1;    
    groupoid[i].expansion=1;
    groupoid[k].expansion=1;
    // done!!
    }  
  } else {
    // it is a loop or an inverse so not in spanning tree if a loop...
    if(groupoid[i].lab[j].loop==1)groupoid[i].span[j]=0; 
    
   }
  }
 }
}

// finally we restore the subset information...  
  for(i=0;i<MAXRANK;i++)subset[i]=subset1[i];
  getsubsettype();
}

int
searchrelgroupoid(int p){
int found=-1,i;
 for(i=0;i<nodeindex;i++){
  if(relgroupoid[i].num==p)found=i;
 }
 return found;
}

void
workoutrelation(int n,int a,int b){
//------------------------------------------
// this function works out the relation subgroupoid
// for node=n in groupoid[.] with vertices a and b...
// it puts the information in relgroupoid[]
// so that relgroupoid[] is the relation
// subgroupoid after this function returns...
// -----------------------------------------
int m,aa,bb,i,j,k,casenum,t,u,v,w,tt;
int temp[MAXNODES]; // used in some of the cases...
int tempcount=0;    // tempcount for indexing temp[.]
int errorhere=0;
int join[MAXRANK];
int joinaabb[MAXRANK];
int lengthvabJ=0; // length of v[{a,b},J]

// save subset information...
for(i=0;i<MAXRANK;i++){
  subset2[i]=subset[i];
  subset[i]=groupoid[n].subs[i];
}
getsubsettype();
subset[a]=1;
subset[b]=1;
getsubsettype();   // work out J\join {a,b}...
lengthvabJ=totallength-subsetlength;
// do only if the length is non-negative...
// if this length is negative or 0 it indicates
// that there is no relation to consider...
// initialize the structure to be empty to begin with...
for(i=0;i<MAXNODES;i++)relgroupoid[i].num=-1;
nodeindex=1;

j=0;
for(i=0;i<MAXRANK;i++){
 if((reducible[i].rank>0)&&(reducible[i].type==X))j=1;
}
if(j==1)lengthvabJ=-1; // if infinite j=1.

if (lengthvabJ>=0){

for(i=0;i<MAXRANK;i++)join[i]=groupoid[n].subs[i];
if ((join[a]==1)||(join[b]==1))errorhere=1;

// errorhere=1 when {a,b} not disjoint from J

if (errorhere==0){

join[a]=1;
join[b]=1;
// now start with subset n...
relgroupoid[0].num=n+1;
relgroupoid[0].expansion=0;
// the first node is initialized now...
// we now must work out two edges coming out of this
// node corresponding to vertices a and b...
for(i=0;i<MAXNODES;i++){
 if ((relgroupoid[i].num>0)&&(relgroupoid[i].expansion==0)){
 // we must expand node i with aa and bb...
 // node i corresponds to node .num-1 in groupoid
 // we must follow the edges aa and bb...
 m=relgroupoid[i].num-1;
 for(j=0;j<MAXRANK;j++){
  relgroupoid[i].adj[j]=0;
  relgroupoid[i].len[j]=0;
  relgroupoid[i].subs[j]=groupoid[m].subs[j];
  if(relgroupoid[i].subs[j]==1){
    joinaabb[j]=0;
  } else joinaabb[j]=join[j];
 }
 // we now work out what aa and bb are...
 // {aa,bb}=join \setminus subs
  k=0;
  for(j=0;j<MAXRANK;j++){
   if((k==1)&&(joinaabb[j]==1)){ bb=j; k++; }
   // order of these lines is important! do not change.
   if((k==0)&&(joinaabb[j]==1)){ aa=j; k++; }
  }
 // aa and bb are now set...
 relgroupoid[i].self=0;
 relgroupoid[i].x=groupoid[m].x;
 relgroupoid[i].y=groupoid[m].x;
 relgroupoid[i].z=groupoid[m].z;
 relgroupoid[i].adj[aa]=groupoid[m].adj[aa];
 relgroupoid[i].adj[bb]=groupoid[m].adj[bb];
 relgroupoid[i].lab[aa]=groupoid[m].lab[aa];
 relgroupoid[i].lab[bb]=groupoid[m].lab[bb];
 relgroupoid[i].len[aa]=groupoid[m].len[aa];
 relgroupoid[i].len[bb]=groupoid[m].len[bb];
 // now we must add to relgroupoid the nodes
 // corresponding to adj[aa] and adj[bb] only 
 // if they are not in the relgroupoid already...
 j=searchrelgroupoid(groupoid[m].adj[aa]); 
 // this returns -1 if not in list
 if (j==-1){
  relgroupoid[nodeindex].num=groupoid[m].adj[aa];
  relgroupoid[nodeindex].expansion=0;
  nodeindex++;
  if(nodeindex>(MAXNODES-1)){ 
    nodeindex=MAXNODES-1; 
    error=1;
    postmessage("Error. Exceeded nodeindex!\n");
  }
  // has been added...
 }
 j=searchrelgroupoid(groupoid[m].adj[bb]);
 if (j==-1){
   relgroupoid[nodeindex].num=groupoid[m].adj[bb];
   relgroupoid[nodeindex].expansion=0;
   nodeindex++;
   if(nodeindex>(MAXNODES-1)){
     nodeindex=MAXNODES-1;
     error=1;
     postmessage("Error. Exceeded nodeindex!\n");
   }
   // has been added...
   } 
  relgroupoid[i].expansion=1;
  }
 }
// finally work out the self fields...
for(i=0;i<nodeindex;i++){
 if(relgroupoid[i].num>0){
 k=0;
 for(j=0;j<MAXRANK;j++)if(relgroupoid[i].adj[j]==relgroupoid[i].num)k++; 
 relgroupoid[i].self=k; 
 }
 }
 // we are here to consider three separate cases
 // for the relgroupoid structure...
 // Case 1: One node with two loops.
 // Case 2: More than one node with loops 1,0,0...0,0,1.
 // Case 3: At least three nodes with loops all 0.
 //
 // we work it out as follows...
 // first test for one node if yes then case 1
 // second test for all loops 0 and at least three nodes then case 3
 // else case 2...
  casenum=2; // this is the default case...
  if(nodeindex==1)casenum=1;
  if(nodeindex>2){
    m=0;
    for(i=0;i<nodeindex;i++){
      if(relgroupoid[i].self==0)m++;
    } 
    if(m==nodeindex)casenum=3;
  }
  // now casenum=1,2 or 3 corresponding to the above
  // three cases...
  // we now must work out what the relation is
  // and add it to the list of label relations...
  
  // we now work out what the relation is...
  // first initialize it to be empty
  argument.len=0;
  argument.tp=casenum;
  argument.eliminate=0;
  m=0;    // m is a counter of length of argument built up...
   
  switch(casenum){
       case 1:
         // one node with two loops...
         v=-1;
         w=-1;
         for(i=0;i<MAXRANK;i++)if(relgroupoid[0].adj[i]!=0){
           v=w;
           w=i;
         }     
         if ((v<0)||(w<0)){
          postmessage("Error. Case 1. v or w not found.\n");
          error=1;
         } else {     
         // v and w are the indices of the two
         // loops...
         // t is the total length done so far...
         t=0;
         u=v;
         // everything ok so far...
          while(t<(2*lengthvabJ)){
            // a braid relation...
            if(m<MAXRELATIONLEN){
              argument.rel[m].loop=relgroupoid[0].lab[u].loop;
              argument.rel[m].label=relgroupoid[0].lab[u].label;
              m++;
              t+=relgroupoid[0].len[u];
              u=(v+w)-u; // alternate...
            }
            else { 
              error=1; 
              postmessage("Error. Case 1. Too high length of relation.\n"); 
              break; }
            }
            if(t==(2*lengthvabJ)){
            argument.len=m;
            } else { 
             
              error=1;
              postmessage("Error. Case 1. Lengths do not add up!\n");
              
            }
           }
         break;
       case 2:
       // here we are in case 2
       // where there are loops 1,0,0,0...0,1
       // we first work out where the two loops are...
       // we let these be v and w...
         v=-1;
         w=-1;
         for(i=0;i<nodeindex;i++)if(relgroupoid[i].self==1){
           v=w;
           w=i;
         }     
         // now v and w are the two loops...
         // t is the total length so far...
         u=v; // u is the current node we are at...
         relgroupoid[v].expansion=0; // this one visited...         
         temp[0]=v;
         tempcount=1; 
         for(i=0;i<MAXRANK;i++)if(relgroupoid[v].adj[i]==relgroupoid[v].num){
          tt=relgroupoid[v].len[i]; 
          argument.rel[m].loop=relgroupoid[v].lab[i].loop;
          argument.rel[m].label=relgroupoid[v].lab[i].label;
         } 
         t=0; // t counts the straight path length       
         m=1;
         while(u!=w){
          // we are at node u
          // and we wish to visit an adjacent node
          // without visiting one already visited...
          
          for(i=0;i<MAXRANK;i++){ 
           if(relgroupoid[u].adj[i]!=0){
            
            j=searchrelgroupoid(relgroupoid[u].adj[i]);
            // j is the index ...
            if(relgroupoid[j].expansion!=0){
             temp[tempcount]=j;
             relgroupoid[j].expansion=0;
             t+=relgroupoid[u].len[i];
             argument.rel[m].loop=relgroupoid[u].lab[i].loop;
             argument.rel[m].label=relgroupoid[u].lab[i].label;
             m++;
             if(m>MAXRELATIONLEN-1){
              m=MAXRELATIONLEN-1;
              error=1;
              postmessage("Error. Case 2. Too high length of relation.\n");
             }
             tempcount++;
             }           
            }
           }
          u=temp[tempcount-1];
          }
          temp[tempcount]=w;
          tempcount++;
          
          u=m-1;
          for(i=0;i<MAXRANK;i++)if(relgroupoid[w].adj[i]==relgroupoid[w].num){
          argument.rel[m].loop=relgroupoid[w].lab[i].loop;
          argument.rel[m].label=relgroupoid[w].lab[i].label;
          tt+=relgroupoid[w].len[i];
          m++;
          if(m>MAXRELATIONLEN-1){
           m=MAXRELATIONLEN-1;
           error=1;
           postmessage("Error. Case 2. Too high length of relation.\n");
           }         
          }
          for(i=u;i>0;i--){
           argument.rel[m].loop=argument.rel[i].loop;
           argument.rel[m].label=-argument.rel[i].label;
           m++;
           if(m>MAXRELATIONLEN-1){
            m=MAXRELATIONLEN-1;
            error=1;
            postmessage("Error. Case 2. Too high length of relation.\n");
           }
          }
    // almost done!!
    // now tt is length of both loops
    // and t is length of straight path once...
    // so the quantity tt+2*t must divide 2*lengthvabJ
        u=m;
        v=tt+2*t;
        w=(int)((2*lengthvabJ)/v);
        
        // if w is more than one then the two loops
        // will appear more than once so it is
        // actually an R2c relation not R2b!!
        // so we fix this by...
        if(w>1)argument.tp=1;
        
        if((w*v)==(2*lengthvabJ)){
           
           // now we must repeat argument w times...
           // argument from 0 to u
           i=0;
           while(i<(w-1)){
           
             for(j=0;j<u;j++){
              argument.rel[m].loop=argument.rel[j].loop;
              argument.rel[m].label=argument.rel[j].label;
              m++;
              if (m>MAXRELATIONLEN-1){
               m=MAXRELATIONLEN-1;
               error=1;
               postmessage("Error. Case 2. Too high length of relation\n");
              }
             
             }
             i++;
           }
          argument.len=m;
        }  
        else {
          // an error has occurred...
          error=1;
          postmessage("Case 2. Error. Lengths do not add up!\n"); 
          }       
         break;
       case 3:
       // here we have a circuit with no loops...
         u=0; // u is the current node we are at...
         relgroupoid[0].expansion=0; // this one visited...         
         temp[0]=u;
         tempcount=1; 
         t=0; // t counts the straight path length       
         m=0;
          //v is last one we were at...
          for(i=0;i<MAXRANK;i++){
                    if(relgroupoid[0].adj[i]!=0)v=searchrelgroupoid(relgroupoid[0].adj[i]);
          }
          
          
         while(t<(2*lengthvabJ)){
          // we are at node u
          // and we wish to visit an adjacent node
          // without visiting one already visited...
          
          for(i=0;i<MAXRANK;i++){ 
           if(relgroupoid[u].adj[i]!=0){
            
            j=searchrelgroupoid(relgroupoid[u].adj[i]);
            // j is the index ...
            if(j!=v){
             temp[tempcount]=j;
             relgroupoid[j].expansion=0;
             t+=relgroupoid[u].len[i];
             argument.rel[m].loop=relgroupoid[u].lab[i].loop;
             argument.rel[m].label=relgroupoid[u].lab[i].label;
             m++;
             if(m>MAXRELATIONLEN-1){
              m=MAXRELATIONLEN-1;
              error=1;
              postmessage("Error. Case 3. Too high length of relation.\n");
             }
             tempcount++;
             }           
            }
           }
          v=u;
          u=temp[tempcount-1];
          }
          if (t!=(2*lengthvabJ)){
           error=1;
           postmessage("Error. Case 3. Lengths do not add up!\n");
           
          } else {
          // lengths add up
           argument.len=m;
           if(u!=0){
            error=1;
            postmessage("Error. Case 3. We do not end up at the beginning.\n");
           }
          }      
         break;
  }
  // we are almost done!!
  // now add the relation...
  // the relation is in argument so add it to the list...
  addrelation(); 
  // done!!
 } else { 
          // an error has occured {a,b} not disjoint from J
          nodeindex=0;
          postmessage("Error. In workoutrelation, {a,b} not disjoint from J\n"); 
          error=1; 
          }

} else { 
         // empty relgroupoid structure since we've
         // hit a non-spherical J\join {a,b} or
         // {a,b} is not disjoint from J...
         nodeindex=0; 
         postmessage("Some relations may be missing.\n");
       }
// restore subset information...
for(i=0;i<MAXRANK;i++){
  subset[i]=subset2[i];
}
getsubsettype();
}
 
void
workoutallrelations(void){
// work out all of the relations...
// first erase all previous relations...
int comp[MAXRANK];    // this is to work out "complement"
int i,j,k,a,b,n;

// first work out the groupoid...
fillingroupoid();
// then erase all relations...
initializerelations();

for(i=0;i<freenode;i++){
 //do for each node in the groupoid[.]
 n=0;
 for(j=0;j<MAXRANK;j++){
  if(j<big.rank){
    if(groupoid[i].subs[j]==0){
    comp[n]=j;
    n++;
    }
  } 
  else comp[j]=0;
 }
 // comp now is the "numbered" complement of length n
 // and n is the cardinality of complement...
 if(n>1){
 // do only if complement cardinality is greater than 1...
 for(a=0;a<(n-1);a++){
  for(b=a+1;b<n;b++){
    if(error==0)workoutrelation(i,comp[a],comp[b]);
    }
   }
 }
}
}

void
cullrelations(void){
// delete all tagged relations.
// ie. those with eliminate bit equal to 1...
int i,j;

i=0;
while(i<relationnum){

 if(relation[i].eliminate==1){
  deleterelation(i);
  i--;
  }  
 i++;
 }
// done!!
}

void freereduceall(void){
int i,j,k,l,change;
// change measures changes to dictionary...

change=1;

while(change==1){
i=0;
change=0;
while(i<relationnum){
// freely reduce it, delete it, then include it.
  argument.len=relation[i].len;
  argument.tp=relation[i].tp;
  argument.eliminate=0;
  for(j=0;j<argument.len;j++){
  argument.rel[j].loop=relation[i].rel[j].loop;
  argument.rel[j].label=relation[i].rel[j].label;
   }  
  convertargument();
  l=searchrelationfrom(i+1);
  // if not found then...
  if(l==-1){
  if(argument.len>0){

  if(argument.len>1){
  replacerelation(i);
  i++;
   } else {
   
   k=argument.rel[0].label;
   if(k<0)k=-k;
   if(argument.rel[0].loop==0)k+=looplabel;
   dict[k].len=0;   
   deleterelation(i);
   change=1;
    }
   } else { deleterelation(i); }
  } else {
  // it was found so we must delete it...
  deleterelation(i);  
   }
  
  }
 }
}

void
putinspanningtree(void){
// put a 1 label in dictionary for
// spanning tree edges...
int i,j,k;

i=0;
while(i<freenode){

 for(j=0;j<MAXRANK;j++){
  if(groupoid[i].adj[j]!=0){
  
   if(groupoid[i].span[j]==1){
   
    k=groupoid[i].lab[j].label;
    if(k>0){
    
     k+=looplabel;
     dict[k].len=0;  
     }
   }
  }
 }
 i++;
 }
}

void 
simplifyR2brelations(void){
// simplify all the R2b relations...
// put the relation in the dictionary...
// rewrite the relations 
// cull the ones that are trivial...
// rewrite the dictionary to get rid of the loop...
int i,j,k,n,l,u,v;

i=0;
while(i<relationnum){

if((relation[i].len>0)&&(relation[i].eliminate==0)&&(relation[i].tp==2)){
  
  // it is an R2b relation...
  // but we can't assume that it starts with
  // a loop...
  // so find the loop first...
  j=0;
  v=-1;
  u=-1;
  while(j<relation[i].len){
   if(relation[i].rel[j].loop==1){ v=u; u=j; }
   j++;
  }
  // u,v are the loops now...
  // u,v cannot be -1 here...
  if((u==-1)||(v==-1)){
  error=1;
  postmessage("Error. No two loops in R2b relation.\n");
  
  } 
    
  if(relation[i].rel[u].label==relation[i].rel[v].label){
  
  // same loops in R2b relation so keep the relation
  
  
  } else 
  {
  k=relation[i].rel[u].label;
  
  for(j=u+1;j<relation[i].len;j++){
  dict[k].rel[j-u-1].loop=relation[i].rel[j].loop;
  dict[k].rel[j-u-1].label=relation[i].rel[j].label;
  }
  dict[k].len=relation[i].len-1;
  for(j=0;j<u;j++){
  dict[k].rel[relation[i].len-u-1+j].loop=relation[i].rel[j].loop;
  dict[k].rel[relation[i].len-u-1+j].label=relation[i].rel[j].label;  
  }
  // we've now put it in the dictionary...
  // so to rewrite all of the relations now...
  j=0;
  while(j<relationnum){
  
   if(relation[j].len<=0)relation[j].eliminate=1;
   
   if(relation[j].eliminate==0){
  
   n=relation[j].len;
   argument.len=n;
   for(l=0;l<n;l++){
    argument.rel[l].loop=relation[j].rel[l].loop;
    argument.rel[l].label=relation[j].rel[l].label;
   }
   // it is copied now to translate it...
   convertargument();
   if(argument.len<=0)relation[j].eliminate=1;  
   else { 
   // re write the relation...
   relation[j].len=argument.len;
   for(l=0;l<argument.len;l++){
   relation[j].rel[l].loop=argument.rel[l].loop;
   relation[j].rel[l].label=argument.rel[l].label;
   }  
   
   }
   // it is now culled...
   }
  j++;
  }
  // all the relevant relations have been culled...
  // so now to trivialise this dict entry...
  dict[k].len=0;
  }
 }
i++;
}
// done!!
}

void 
simplifyR2arelations(void){
// simplify all the R2a relations...
// put the relation in the dictionary...
// rewrite the relations 
// cull the ones that are trivial...
// rewrite the dictionary to get rid of the loop...
int i,j,k,n,l,u,v,inverse;

i=0;
while(i<relationnum){

if((relation[i].len>1)&&(relation[i].eliminate==0)&&(relation[i].tp==3)){
  
  // it is an R2a relation...
  // but we need to see if it is circuit=1
  // or circuit^2=1.
  // only reduce if it is circuit=1.
  u=0;
  for(j=1;j<relation[i].len;j++){
   if((relation[i].rel[0].loop==relation[i].rel[j].loop)&&(relation[i].rel[0].label==relation[i].rel[j].label))u++;
  }
  
  // do reduction only if u==0;
      
  if(u!=0){
  
  // it is a circuit^2=1 R2a relation so keep it...
  
  
  } else 
  {
  
  inverse=0;
  k=relation[i].rel[0].label;
  if(k<0){ k=-k; inverse=1; }
    
  // copy to argument...
  argument.len=relation[i].len-1;
  for(j=1;j<relation[i].len;j++){
   argument.rel[j-1].loop=relation[i].rel[j].loop;
   argument.rel[j-1].label=relation[i].rel[j].label;
  }
  
  if(inverse==0)invertargument();
  
  // now argument is to be put in dictionary...
  k+=looplabel; 
  
  for(j=0;j<argument.len;j++){
  dict[k].rel[j].loop=argument.rel[j].loop;
  dict[k].rel[j].label=argument.rel[j].label;
  }
  dict[k].len=argument.len;
  
  // we've now put it in the dictionary...
  // so to rewrite all of the relations now...
  j=0;
  while(j<relationnum){
  
   if(relation[j].len<=0)relation[j].eliminate=1;
   
   if(relation[j].eliminate==0){
  
   n=relation[j].len;
   argument.len=n;
   for(l=0;l<n;l++){
    argument.rel[l].loop=relation[j].rel[l].loop;
    argument.rel[l].label=relation[j].rel[l].label;
   }
   // it is copied now to translate it...
   convertargument();
   if(argument.len<=0)relation[j].eliminate=1;  
   else { 
   // re write the relation...
   relation[j].len=argument.len;
   for(l=0;l<argument.len;l++){
   relation[j].rel[l].loop=argument.rel[l].loop;
   relation[j].rel[l].label=argument.rel[l].label;
   }  
   
   }
   // it is now culled...
   }
  j++;
  } 
  // all the relevant relations have been culled...
  // so now to trivialise this dict entry...
  dict[k].len=0;
  }
 }
i++;
}
// done!!
}

void
showgroupoid(void){
// shows the groupoid...
 fillingroupoid();
 printgroupoid();
 printf("Status: "); 
 if(posted==1)printmessage(); else printf("Ok.\n");
}

void
showpresentation(void){
// works out and shows the presentation...
  printf("\nPresentation of N_J...\n");
  fillingroupoid();
  printf("Initial relations...\n");
  workoutallrelations(); 
  printrelationsshort();
  printf("Status: "); 
  if(posted==1)printmessage(); else printf("Ok.\n");
  initializedict();
  putinspanningtree();
  freereduceall();
  printf("After simplifying R2b relations...\n");
  simplifyR2brelations();
  cullrelations(); 
  printrelationsshort();
  printf("Status: ");
  if(posted==1)printmessage(); else printf("Ok.\n");
  //after spanning tree...
  printf("After simplifying R2a relations\n");
  simplifyR2arelations();
  cullrelations();
  printrelationsshort();
  printf("Status: ");
  if(posted==1)printmessage(); else printf("Ok.\n");
  //free reduce all
  printf("Reducing...\n");
  freereduceall();
  printf("Final presentation...\n");
  printrelations();
  printdict();
  printf("Status: ");
  if(posted==1)printmessage(); else printf("Ok.\n");
}

void
Renderg(GLenum mode)
{
  GLint i,j,s,n,letter;
  float x,y,z;
  char buffer[MAXRANK*10];
  char auxbuf[20];
  char *mybuffer;
  
   for(i=0;i<freenode;i++){
      if (mode == GL_SELECT) {
      glLoadName(i);
      }
    
    glColor3f(colour[0][0],colour[0][1],colour[0][2]); 
    glPushMatrix();
    glTranslatef(groupoid[i].x,groupoid[i].y,groupoid[i].z);
    glNewList(i, GL_COMPILE);
    glutSolidCube(0.2);
    glEndList();
    glCallList(i);
    glPopMatrix();   
    
    glColor3f(colour[3][0],colour[3][1],colour[3][2]);
    x=groupoid[i].x;
    y=groupoid[i].y;
    z=groupoid[i].z;
    itoa(groupoid[i].self,&mybuffer[0],10);
    stroke_output(x+0.5,y-0.5,z, mybuffer); 
    
    if(showlabels==1){
    
    glColor3f(colour[5][0],colour[5][1],colour[5][2]);
    itoa(groupoid[i].num,&mybuffer[0],10);
    stroke_outputsmall(x,y-0.5,z,mybuffer);
     }
    
    
    }
    
    if (mode==GL_SELECT){ 
      glLoadName(freenode);
    }
    for(i=0;i<freenode;i++){
      for(j=0;j<MAXRANK;j++){
        if (groupoid[i].adj[j]!=0){
          s=groupoid[i].adj[j];
          if (i+1<s){
          
          if((spanning==1)&&(groupoid[i].span[j]==1)){
            glColor3f(colour[0][0],colour[0][1],colour[0][2]);
          } else { 
          glColor3f(colour[1][0],colour[1][1],colour[1][2]);
          }
          
          x=(groupoid[i].x+groupoid[s-1].x)/2;
          y=(groupoid[i].y+groupoid[s-1].y)/2;
          z=(groupoid[i].z+groupoid[s-1].z)/2;
                    
          glBegin(GL_LINES);
          glVertex3f(groupoid[i].x,groupoid[i].y,groupoid[i].z);
          glVertex3f(groupoid[s-1].x,groupoid[s-1].y,groupoid[s-1].z);
          glEnd();  
         
         if(showlabels==1){ 
          glColor3f(colour[4][0],colour[4][1],colour[4][2]);
          itoa(groupoid[i].lab[j].label,&mybuffer[0],10);
          stroke_outputsmall(x,y+0.2,z,mybuffer); }
          
          }      
        }
      }
    }
   glColor3f(colour[2][0],colour[2][1],colour[2][2]);
   x=groupoid[0].x;
   y=groupoid[0].y;
   z=groupoid[0].z;
   if (error==0){
   itoa(freenode,&mybuffer[0],10);
   stroke_output(x-6.0,y,z,mybuffer);
   itoa(looplabel-1,&mybuffer[0],10);
   stroke_output(x-4.0,y,z,mybuffer);
   itoa(edgelabel-1,&mybuffer[0],10);
   stroke_output(x-2.0,y,z,mybuffer);
   } else { 
   stroke_output(x-5.0,y,z,"Error");
   }
   glutSwapBuffers();  
}

void display(void)
{
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glTranslatef((-((float)(MAXRANK)/2)-originx), originy, -distance);
    glRotatef(roll,0.0,1.0,0.0);
    glRotatef(yaw,0.0,0.0,1.0);
    glRotatef(pitch,1.0,0.0,0.0);
    glPixelStorei(GL_UNPACK_ALIGNMENT, 1); 
    Render(GL_RENDER);
}

void myReshape(int w, int h)
{
    glViewport(0, 0, w, h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluPerspective(60.0, 1.0*(GLfloat)w/(GLfloat)h, 0.01, (float)MAXRANK*2);
    glGetIntegerv(GL_VIEWPORT, vp);
    windH=h; windW=w;
    display();
}

void displayg(void){
float ox,oy;

    fillingroupoid();
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    ox=-distanceg/2.0;
    oy=-distanceg/2.0;
    glTranslatef(-originxg+ox, oy+originyg, -distanceg);
    glRotatef(rollg,0.0,1.0,0.0);
    glRotatef(yawg,0.0,0.0,1.0);
    glRotatef(pitchg,1.0,0.0,0.0);
    glPixelStorei(GL_UNPACK_ALIGNMENT, 1); 
    Renderg(GL_RENDER);
}

void idle(void){
int w;
w=glutGetWindow();
glutSetWindow(win);
display();
glutSetWindow(wing);
displayg();
glutSetWindow(w);
}

GLint
DoSelect(GLint x, GLint y)
{
  GLint hits;

  glSelectBuffer(MAXSELECT, selectBuf);
  glRenderMode(GL_SELECT);
  glInitNames();
  glPushName(-1);

  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPickMatrix(x, windH - y, 5, 5, vp);
  gluPerspective(60.0, 1.0*(GLfloat)windW/(GLfloat)windH, 0.01, (float)MAXRANK*2);   

  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity(); 
  glTranslatef((-((float)(MAXRANK)/2)-originx), originy, -distance);
  glRotatef(roll,0.0,1.0,0.0);
  glRotatef(yaw,0.0,0.0,1.0);
  glRotatef(pitch,1.0,0.0,0.0);
  glPixelStorei(GL_UNPACK_ALIGNMENT, 1); 
   //myReshape(windW,windH);
  Render(GL_SELECT);
  
  hits = glRenderMode(GL_RENDER);
  if (hits <= 0) {
    return -1;
  }
  return selectBuf[(hits - 1) * 4 + 3];
}

GLint
DoSelectg(GLint x, GLint y)
{
  GLint hits;
  float ox,oy;

  glSelectBuffer(MAXSELECT, selectBuf);
  glRenderMode(GL_SELECT);
  glInitNames();
  glPushName(-1);

  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPickMatrix(x, windHg - y, 5, 5, vpg);
  gluPerspective(60.0, 1.0*(GLfloat)windWg/(GLfloat)windHg, 0.01, (float)MAXRANK*2);   

  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity(); 
  ox=-distanceg/2.0;
  oy=-distanceg/2.0;
  glTranslatef(-originxg+ox,oy+originyg, -distanceg);
  glRotatef(rollg,0.0,1.0,0.0);
  glRotatef(yawg,0.0,0.0,1.0);
  glRotatef(pitchg,1.0,0.0,0.0);
  glPixelStorei(GL_UNPACK_ALIGNMENT, 1); 
   //myReshape(windW,windH);
  Renderg(GL_SELECT);
  
  hits = glRenderMode(GL_RENDER);
  if (hits <= 0) {
    return -1;
  }
  return selectBuf[(hits - 1) * 4 + 3];
}

void
Motion(int x,int y){
 if(motion==1) {
 originx=(float)(mouseX0-x)/distance;
 originy=(float)(mouseY0-y)/distance;
 } else
 if (motion==2){
 
 if((mouseX0-x)*(mouseX0-x)<900)yaw+=(float)0.05*(mouseY0-y)/distance; 
  else
 if((mouseY0-y)*(mouseY0-y)<900)roll-=(float)0.05*(mouseX0-x)/distance; 
  else pitch+=(float)0.00005*(mouseX0-x)*sqrt((mouseX0-x)*(mouseX0-x)+(mouseY0-y)*(mouseY0-y))/distance;
 
 
 }
 display();
}

void
Motiong(int x,int y){
 if(motiong==1) {
 originxg=(float)(mouseX0g-x)/distanceg;
 originyg=(float)(mouseY0g-y)/distanceg;
 } else
 if (motiong==2){
 
 if((mouseX0g-x)*(mouseX0g-x)<900)yawg+=(float)0.05*(mouseY0g-y)/distanceg; 
  else
 if((mouseY0g-y)*(mouseY0g-y)<900)rollg-=(float)0.05*(mouseX0g-x)/distanceg; 
  else pitchg+=(float)0.00005*(mouseX0g-x)*sqrt((mouseX0g-x)*(mouseX0g-x)+(mouseY0g-y)*(mouseY0g-y))/distanceg;
 }
 displayg();
}

void
Mouse(int button, int state, int mouseX, int mouseY)
{
   GLint i,j,hit;
   TYPE* t;
      
   t=&big;
   if (state == GLUT_DOWN) {
    hit = DoSelect((GLint) mouseX, (GLint) mouseY);
    if (hit != -1) {
      if (button == GLUT_RIGHT_BUTTON) {
        hit=selectBuf[3];
        if (hit==t->rank){ mouseX0=mouseX; mouseY0=mouseY; }
        if (hit<t->rank){
             subset[hit]=1-subset[hit]; 
             printsubsettype();}
      } 
      if (button==GLUT_LEFT_BUTTON){
        hit=selectBuf[3];
        pivot=t->rank;
        if (hit<t->rank){
            pivot=hit; motion=0; 
        for(i=0;i<MAXRANK;i++)subset0[i]=subset[i];
        subset[pivot]=1;  // add the pivot to subset...
        getsubsettype();
        // now reducible has the types of the irreducible
        // components of the union of selected subset and
        // the pivot vertex...
        if(totallength>=0)vailength=totallength-subsetlength;
        else vailength=-1;
        printf("Length of v[a,I]: %d\n",vailength);
        getnextsubset();
        // work out the appropriate transformation...
        getsubsettype();
            
            }
        else 
        if (hit==t->rank){ 
          mouseX0=mouseX;
          mouseY0=mouseY;
          motion=1;
        } 
        
        
      }   
  } else { mouseX0=mouseX; mouseY0=mouseY; motion=2; }
  } else motion=0;
  
  myReshape(windW,windH);  
} 

void myReshapeg(int w, int h)
{
    glViewport(0, 0, w, h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluPerspective(60.0, 1.0*(GLfloat)w/(GLfloat)h, 0.01, (float)MAXRANK*2);
    glGetIntegerv(GL_VIEWPORT, vpg);
    windHg=h; windWg=w;
    displayg();
}

void
Mouseg(int button, int state, int mouseX, int mouseY)
{
   GLint i,j,hit;
   TYPE* t;
   
   t=&big;
   if (state == GLUT_DOWN) {
    hit = DoSelectg((GLint) mouseX, (GLint) mouseY);
    if (hit != -1) {
       
      if (button==GLUT_LEFT_BUTTON){
        hit=selectBuf[3];
        if ((hit>=0)&&(hit<freenode)){
        motiong=0; 
        // display the information for this node...
  i=hit;
  printf("  Node number: %d \n",groupoid[i].num);
  printf("   Subset: { ");
  for(j=0;j<big.rank;j++)if(groupoid[i].subs[j]==1)printf(" %d ",j);
  printf("}\n");
  printf("   Self: %d \n",groupoid[i].self);
  printf("   Adjacent: { ");
  for(j=0;j<MAXRANK;j++)if(groupoid[i].adj[j]!=0)printf(" %d ",groupoid[i].adj[j]);
  printf("}\n   Lengths:  { ");
  for(j=0;j<MAXRANK;j++)if(groupoid[i].adj[j]!=0)printf(" %d ",groupoid[i].len[j]);
  printf("}\n   Labels :  { ");
  for(j=0;j<MAXRANK;j++)if(groupoid[i].adj[j]!=0){printlabel(groupoid[i].lab[j]); printf(" ");}
  printf("}\n   Span   :  { ");
  for(j=0;j<MAXRANK;j++)if(groupoid[i].adj[j]!=0)printf(" %d ",groupoid[i].span[j]);
  printf("}\n\n");
            }
        else 
        if (hit>=freenode){ 
          mouseX0g=mouseX;
          mouseY0g=mouseY;
          motiong=1;
        } 
        
        
      }   
  } else { mouseX0g=mouseX; mouseY0g=mouseY; motiong=2; }
  } else motiong=0;
  
  myReshapeg(windWg,windHg);  
} 

void keyboard(unsigned char key,int x,int y){
int i,j,s;

if (key==27)exit(0);
else
if (key=='a'){

  for(i=0;i<MAXRANK;i++)subset[i]=0;
  if (oldkey!=key){
    big.rank=6; 
    big.type=A; 
    big.parameter=3;
     } 
    else { big.rank++; 
     if (big.rank>MAXRANK)big.rank=1;
    }
    inittype(&big);  
    getsubsettype();
}
else
if (key=='b'){
for(i=0;i<MAXRANK;i++)subset[i]=0;
    if(oldkey!=key){
    big.rank=6; 
    big.type=B; 
    big.parameter=3;
    } else { big.rank++; 
     if (big.rank>MAXRANK)big.rank=2; }
    inittype(&big);  
    getsubsettype();
}
else
if (key=='d'){

for(i=0;i<MAXRANK;i++)subset[i]=0;
   if (oldkey!=key){
    big.rank=6; 
    big.type=D; 
    big.parameter=3;
    } else { big.rank++; 
     if (big.rank>MAXRANK)big.rank=4; }
    inittype(&big); 
    getsubsettype(); 
}
else
if (key=='e'){
for(i=0;i<MAXRANK;i++)subset[i]=0;
  if (oldkey!=key){
    big.rank=6; 
    big.type=E; 
    big.parameter=3;
    } else { big.rank++; 
     if (big.rank>8)big.rank=6;  }
    inittype(&big); 
    getsubsettype(); 
}
else
if (key=='f'){
for(i=0;i<MAXRANK;i++)subset[i]=0;
    big.rank=4; 
    big.type=F; 
    big.parameter=3;
    inittype(&big); 
    getsubsettype(); 
}
else
if (key=='h'){
for(i=0;i<MAXRANK;i++)subset[i]=0;
  if(oldkey!=key){
    big.rank=3; 
    big.type=H; 
    big.parameter=3;
    } else { 
     big.rank=7-big.rank;
    }
    inittype(&big); 
    getsubsettype(); 
}
else
if (key=='i'){
    for(i=0;i<MAXRANK;i++)subset[i]=0;
    big.rank=2; 
    big.type=I; 
    if (oldkey!=key){ big.parameter=5; }
    else { big.parameter++; 
     if (big.parameter>MAXRANK)big.parameter=5; }
    inittype(&big); 
    getsubsettype(); 
} 
else
if (key=='+'){ 
   for(i=0;i<MAXRANK;i++)subset[i]=0;
   big.rank++; if(big.rank>MAXRANK)big.rank=1; 
   inittype(&big);
   getsubsettype();
   }
else
if (key=='-'){
   for(i=0;i<MAXRANK;i++)subset[i]=0;
   big.rank--; if(big.rank<1)big.rank=MAXRANK;
   inittype(&big);
   getsubsettype();
} 
else
if (key=='p'){
   showpresentation();  
}
else
if (key=='g'){
   showgroupoid();
}
else
if (key=='r'){
// raw printout of relations without reducing them...
  printf("\nRaw presentation of N_J...\n");
  fillingroupoid();
  printf("Initial relations...\n");
  workoutallrelations(); 
  printrelations();
  printf("Status: "); 
  if(posted==1)printmessage(); else printf("Ok.\n");
}
else
if (key=='A'){
for(i=0;i<MAXRANK;i++)subset[i]=0;
  if (oldkey!=key){
    big.rank=3; 
    big.type=X; 
    big.parameter=A;
     } 
    else { big.rank++; 
     if (big.rank>MAXRANK){ big.rank=2; key=' ';}
    }
    inittype(&big);  
    getsubsettype();
}
else
if (key=='B'){
for(i=0;i<MAXRANK;i++)subset[i]=0;
  if (oldkey!=key){
    big.rank=6; 
    big.type=X; 
    big.parameter=B;
     } 
    else { big.rank++; 
     if (big.rank>MAXRANK)big.rank=3;
    }
    inittype(&big);  
    getsubsettype();
}
else
if (key=='C'){
for(i=0;i<MAXRANK;i++)subset[i]=0;
  if (oldkey!=key){
    big.rank=6; 
    big.type=X; 
    big.parameter=C;
     } 
    else { big.rank++; 
     if (big.rank>MAXRANK)big.rank=3;
    }
    inittype(&big);  
    getsubsettype();
}
else
if (key=='D'){
for(i=0;i<MAXRANK;i++)subset[i]=0;
  if (oldkey!=key){
    big.rank=6; 
    big.type=X; 
    big.parameter=D;
     } 
    else { big.rank++; 
     if (big.rank>MAXRANK)big.rank=3;
    }
    inittype(&big);  
    getsubsettype();
}
else
if (key=='E'){
for(i=0;i<MAXRANK;i++)subset[i]=0;
  if (oldkey!=key){
    big.rank=7; 
    big.type=X; 
    big.parameter=E;
     } 
    else { big.rank++; 
     if (big.rank>9)big.rank=7;
    }
    inittype(&big);  
    getsubsettype();
}
else
if (key=='F'){
for(i=0;i<MAXRANK;i++)subset[i]=0;
    big.rank=5; 
    big.type=X; 
    big.parameter=F;
    inittype(&big);  
    getsubsettype();
}
else
if (key=='G'){
for(i=0;i<MAXRANK;i++)subset[i]=0;
   big.rank=3;
   big.type=X;
   big.parameter=G;
   inittype(&big);
   getsubsettype();
}
else
if (key=='['){
 for(i=0;i<MAXRANK;i++)subset[i]=0;
 hyperbolic++;
 if(hyperbolic>71)hyperbolic=0;
 big.rank=hyprank[hyperbolic];
 big.type=X;
 big.parameter=HYPBASE+hyperbolic;
 inittype(&big);
 getsubsettype();
}
else
if(key==']'){
 for(i=0;i<MAXRANK;i++)subset[i]=0;
 hyperbolic--;
 if(hyperbolic<0)hyperbolic=71;
 big.rank=hyprank[hyperbolic];
 big.type=X;
 big.parameter=HYPBASE+hyperbolic;
 inittype(&big);
 getsubsettype();
}
else
if (key=='c'){
 for(i=0;i<MAXRANK;i++)subset[i]=0;
 big.rank=CUSTOMRANK;
 big.type=X;
 big.parameter=CUSTOM;
 inittype(&big);
 getsubsettype();
}
else
if (key=='s'){
 // toggle display of spanning tree...
 spanning=1-spanning;
}
else
if(key=='l'){
// toggle display of labels...
 showlabels=1-showlabels;
}
else
if(key=='x'){
 distance=MAXRANK;
 distanceg=MAXRANK/2;
 originx=0.0;
 originy=0.0;
 originxg=0.0;
 originyg=0.0;
 yaw=0.0;
 pitch=0.0;
 roll=0.0;
 yawg=0.0;
 pitchg=0.0;
 rollg=0.0;
}
oldkey=key;
display();
} 

int 
main(int argc, char** argv)
{
int i;
float si,memory;
    si=(float)sizeof(int);
    memory=si*(MAXRELATIONS*(1+MAXRELATIONLEN*2)+MAXDICT*(1+2*MAXRELATIONLEN)+MAXRANK*(MAXRANK+4)+MAXNODES*16);
    printf("Coxeter group viewer, Nov. 2004\n");
    printf("Author: Mauro Grassi.\n\n");
    printf("Current limits are: \n");
    printf("Maximum rank: %d\n",MAXRANK);
    printf("Maximum nodes in groupoid: %d\n",MAXNODES);
    printf("Maximum number of labels: %d\n",MAXDICT);
    printf("Maximum length of a relation: %d\n",MAXRELATIONLEN);
    printf("Maximum number of relations: %d\n",MAXRELATIONS);
    printf("Estimated total data memory used: %.1f MBytes\n",memory/1000000.0);
    printf("\n");    
    
    
    for(i=0;i<MAXRANK;i++)subset[i]=0;
    big.rank=6; 
    big.type=A; 
    big.parameter=3;
    inittype(&big);  
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
    glutInitWindowPosition(1000-windWg,728-windHg);
    glutInitWindowSize(windWg,windHg);
    wing=glutCreateWindow("Groupoid Viewer, Mauro Grassi Nov. 2004");
    myinit();
    glutDisplayFunc(displayg);
    glutReshapeFunc(myReshapeg);
    glutMouseFunc(Mouseg);
    glutMotionFunc(Motiong);
    glutInitWindowPosition(0,0);
    glutInitWindowSize(windW,windH);
    win=glutCreateWindow("Coxeter v[a,I] Viewer, Mauro Grassi Nov. 2004");
    myinit();
    glutReshapeFunc(myReshape);
    glutDisplayFunc(display);
    glutMouseFunc(Mouse);
    glutMotionFunc(Motion);
    glutKeyboardFunc(keyboard); 
    glutIdleFunc(idle);
    glutSetWindow(win);
    glutMainLoop(); 
    return 0;  /* ANSI C requires main to return int. */
}
 
