/* This C++ program counts and displays the 1 1 2 5 15 54 212 908 4011 18260 incongruental non-labbeled undirected trees on square lattices with n=1,2,3,.., nodes and n-1 edges. http://oeis.org/A056841 Richard J. Mathar, 2006-04-09 */ #include <limits.h> #include <iostream> #include <vector> #include <algorithm> #include <unistd.h> using namespace std ; /** A node is one point on the square lattice. * A node contains simply the two integer coordinates of a point on the * square lattice. */ class node { public: / /* xy[0] is the x, and xy[1] the y-component of the coordinates */ short xy[2] ; node() { xy[0]=xy[1] = SHRT_MIN ; } /** Constructor with a pair of coordinates. * @param x the x coordinate of the point * @param y the y coordinate of the point */ node(short x, short y) { xy[0]=x ; xy[1] = y ; } /** Constructor with a given original node translated into a direction. * @param n the node to start from * @param dir a parameter from 0 to 4 to indicate in which * direction of the node \c n this node that is created here * is to be placed. */ node(const node & n, int dir) { /** For values of the \c dir argument from 0 to 3 * we create the new node one step into the East (+x), * North (+y), West (-x) or South (-y) from the original one. */ switch(dir) { case 0: xy[0] = n.xy[0]+1; xy[1] = n.xy[1] ; break ; case 1: xy[0] = n.xy[0]; xy[1] = n.xy[1]+1 ; break ; case 2: xy[0] = n.xy[0]-1; xy[1] = n.xy[1] ; break ; case 3: xy[0] = n.xy[0]; xy[1] = n.xy[1]-1 ; break ; } } /** Create a new node derived from this one by one of the eight * symmetry operations of the square. * @param mode a value from 0 to 7 indicative of which symmetry operation * of the square is used to place the imaged node. */ node congru(const int mode) const { /** parameters \c mode from 0 to 7 generate a new node that is * either the original one (mode=0), rotated around the (0,0) center * of the coordinates by multiples of 90 degrees (mode=1 to 4), or * created in a two-step process of first reflecting the node across the * main (x=y) diagonal of the coordinates and then rotating it to * a place away a multiple of 90 degrees. */ switch(mode) { case 0 : return *this ; break ; case 1 : // 90 degrees rotated return node(-xy[1],xy[0]) ; break ; case 2 : // 180 degrees rotated return node(-xy[0],-xy[1]) ; break ; case 3 : // 270 degrees rotated return node(xy[1],-xy[0]) ; break ; case 4 : // mirror main diagonal return node(xy[1],xy[0]) ; break ; case 5 : // mirror main diagonal then 90 degrees rotated return node(-xy[0],xy[1]) ; break ; case 6 : // mirror main diagonal then 180 degrees rotated return node(-xy[1],-xy[0]) ; break ; case 7 : // mirror main diagonal then 270 degrees rotated return node(xy[0],-xy[1]) ; break ; } } protected: private: } ; /** Node ordering operator. * This introduces a ordering of a list of nodes by ordering with respect * to the y coordinate (smaller y is interpreted as smaller node), and with * respect to the x coordinate if the y coordinates are equal. * @param n1 the node left from the minor sign. * @param n2 the node right from the minor sign. */ bool operator< (const node &n1, const node &n2) { if (n1.xy[1] < n2.xy[1]) return true ; else if (n1.xy[1] > n2.xy[1]) return false ; else if (n1.xy[0] < n2.xy[0]) return true ; else return false ; } /** Node ordering operator. * Two nodes are equal if they have the same x and the same y coordinate. * @param n1 the node left from the equal sign. * @param n2 the node right from the equal sign. */ bool operator== (const node &n1, const node &n2) { return (n1.xy[0] == n2.xy[0] && n1.xy[1] == n2.xy[1]) ; } /** The edge class is a pair of nodes connected by an edge. * The connectivity is implicit, and only the two coordinates of the * pair of nodes are stored. */ class edge { public: /* The pair of nodes, assumed ordered such that nodes[0] < nodes[1]. * To form an edge, nodes[0]!=nodes[1] in all cases. */ node nodes[2] ; /** The constructor with a pair of nodes. * @param n1 one of the two nodes. * @param n2 the other one of the two nodes. */ edge(const node n1, const node n2) { /* To ensure correct ordering, we call an interface function * of this class which swaps the two arguments if necessary */ init(n1,n2) ; } /** A constructor with a node and an indicator of a direction * in which the other node is found. * @param n one node to be used for the edge. * @param dir a direction indicator from 0 to 3. This ressembles a * wind rose indicator for the four potential directions on the * square lattice. */ edge(const node n, int dir) { #ifdef OLD /* This is in essence the slow version of the construction * since it involves an additional comparison with the init() function. */ switch(dir) { case 0: init(n, node(n.xy[0]+1,n.xy[1])) ; break ; case 1: init(n, node(n.xy[0],n.xy[1]+1)) ; break ; case 2: init(n, node(n.xy[0]-1,n.xy[1])) ; break ; case 3: init(n, node(n.xy[0],n.xy[1]-1)) ; break ; } #else /* For dir values from 0 to 3 we create the pairing node in * the East, North, West and South direction (+x,+y,-x,-y) from the node \c n that * was provided explicitly. */ switch(dir) { case 0: nodes[0] = n ; nodes[1] = node(n.xy[0]+1,n.xy[1]) ; break ; case 1: nodes[0] = n ; nodes[1] = node(n.xy[0],n.xy[1]+1) ; break ; case 2: nodes[1] = n ; nodes[0] = node(n.xy[0]-1,n.xy[1]) ; break ; case 3: nodes[1] = n ; nodes[0] = node(n.xy[0],n.xy[1]-1) ; break ; } #endif } /** This constructs a congruent version of the existing node. * @param mode a value form 0 to 7 used in nodes.congru() representing * any of the eight ways (including the indentity/copy) to generate an * image of the node by rotation and reflection. * @return the new edge rotated and/or reflected with respect to the * origin. */ edge congru(const int mode) const { return edge( nodes[0].congru(mode), nodes[1].congru(mode)) ; } /** Translate the edge. * On return, the edge is moved by \c dx and dy into the two directions on the square grid. * @param dx the integer value for motion into the x direction * @param dy the integer value for motion into the y direction */ void shift(short dx, short dy) { nodes[0].xy[0] += dx ; nodes[1].xy[0] += dx ; nodes[0].xy[1] += dy ; nodes[1].xy[1] += dy ; } protected: private: /** Copy throug the nodes ensuring that the smaller node is placed first. * @param n1 the first of the nodes to be used for the new edge * @param n2 the second of the nodes to be used for the new edge */ void init(const node & n1, const node & n2) { /* The node comparison operator is involved to ensure a consistent * ordering of the two nodes in the short vector of nodes within the edge. */ if ( n1 < n2) { nodes[0]=n1; nodes[1]=n2; } else { nodes[0]=n2; nodes[1]=n1; } } } ; /** Edge ordering operator. * The edge with the smaller first node is considered the smaller of the * two edges. If the first node is shared between the two edges, the * comparison is done on the second node. * @param e1 the first of the two edges to be compared. * @param e2 the second of the two edges to be compared. * @return true if the first edge is considered the smaller one. */ bool operator< (const edge &e1, const edge &e2) { if ( e1.nodes[0] < e2.nodes[0]) return true ; else if ( e1.nodes[0] == e2.nodes[0]) { if ( e1.nodes[1] < e2.nodes[1]) return true ; else return false ; } else return false ; } /** Edge ordering operator. * Edges are the same if the nodes are the same. * @param e1 the first of the two edges to be compared. * @param e2 the second of the two edges to be compared. * @return true if the edges are the same. */ bool operator== (const edge &e1, const edge &e2) { /* Since edges are sorted in their nodes, it's sufficient to compare the first and the * second nodes separately, not all four pairs of them... This is an attempt to speed * up the program. */ return ( e1.nodes[0] == e2.nodes[0] && e1.nodes[1] == e2.nodes[1]) ; } /** A tree is a list/vector of edges, all of which are different, and which * share nodes such that one can transverse the graph on a trail that * may reach each of the nodes.. */ class tree { public: /* The vector of all edges is the only full representation of the tree. */ vector<edge> edges ; /* The \c bbox contains a bounding box, such that bbox[0] is the smallest x, * bbox[1] the smallest y, bbox[2] the largest x and bbox[3] the largest y value of * all nodes of the \c edges. */ short bbox[4] ; /* The default constructor creates an empty tree. * It is useful to generated vectors of trees and has no other significance. */ tree() { } /** The function congru() created a congruental image of the current tree * @param mode a number from 0 to 7 in the sense of edge::congru() which * specifies which of the 8 symmetry operations (mode=0 is a copy) of the square * is used to generate the new tree. * @return the rotated/reflected tree */ tree congru(const int mode) const { tree resul ; /* The algorithm is simply to create the image of each edge separately, * and add it to the new tree with addedge(). */ for(int e=0; e < int(edges.size()); e++) { const edge cedge = edges[e].congru(mode) ; resul.addedge(cedge) ; } /* The typical rotation would have created an image in any quadrant of * the coordinate system. For quicker comparison of congruent pairs (which are * to be eliminated according to the algorithm), we translate each tree such * it is placed in the upper-right quadrant with the smallest x and the smallest * y coordinate of all nodes being both 0. */ resul.normalize() ; return resul ; } /** Create one more edge in the tree. This is done without * any further check such as ensuring that the edge is connected to the * already existing truncs or such as ensuring it is not yet part of the tree. * @param e the edge to be added */ void addedge(const edge & e) { edges.push_back(e) ; } /** Investigate whether the two trees are congruent. * The result is true if all the edges of this instance and the other instance * can be placed on each other by translation, rotation or mirror operation. * @param oth the other tree that is to be compared with the current instance. * It is assumed that its \c bbox parameter is up-to-date. * @return true if this tree and the other one are congruent. */ bool iscongru(const tree & oth) { /* To avoid some overhead during the comparison, we compared some * rough parameters first: if the numbers of edges in the two trees differ, * the trees differ. */ if ( edges.size() != oth.edges.size() ) return false ; /* The congruent operations of rotations and reflections can swap * the bounding box x and y coordinates, but not change them otherwise. * So we conclude that the two trees differ if this test of matching fails. */ if ( bbox[2] != oth.bbox[2] && bbox[2] != oth.bbox[3]) return false ; if ( bbox[3] != oth.bbox[2] && bbox[3] != oth.bbox[3]) return false ; /* The final more constly operation of deciding on congruence is to create * all 8 congruent images of \c oth and compare each of them with the * current instance. To speed up the comparison, we sort the edge list * of the current instance and each edge list of the 8 images created of the other * instance. This ought result in a comparison time which is only linear in the * number of edges. */ sort(edges.begin(),edges.end() ) ; /* This is the loop over all possible 8 congruental versions of the other tree. */ for(int c=0; c < 8 ;c++) { /* The temporary mage copy of the other tree by calling congru(). */ tree imag = oth.congru(c) ; /* Sorting the edge list of the image copy and figuring out whether the * image lists are then the same is done with two calls of the STL library. * If we found that any of the 8 versions equals the current instance, * we do not need to investigate the other versions and may conclude * that both trees are congruental. */ sort(imag.edges.begin(),imag.edges.end() ) ; if ( equal(edges.begin(),edges.end(),imag.edges.begin()) ) return true ; } /* If we have fallen through the loop over the 8 images, none of these * was congruental and the result is a 'no'. */ return false ; } /** A function to generate all trees with one more edge than the existing instance. * The current instance is the parent from which all possible ways of adding * one more edge to any existing node are tested. The children are reduced to * a set of mutually incongruental trees and returned in a vector. * @return a vector of the incongruental trees generated. */ vector<tree> childs() { /* \c resul collects all new trees with one more edge than the current tree */ vector<tree> resul ; /* For some efficiency, we generate a list of all nodes of the existing tree, * which are typically shared by between 1 and 4 edges. The idea is to spawn * new edges starting from any of these into any of the four directions * admitted by the square lattice. */ const vector<node> thisnod = allnodes() ; //cout << thisnod.size() << " nodes " << edges.size() << " edges\n" ; /* In a loop over all existing nodes we can generate all child trees * that have the same skeleton as the current one, but one more edge. * This works since a tree is a connected graph, so edges cannot be isolated. */ for(int n=0; n < int(thisnod.size()) ; n ++) { /* For each of the four directions of the node we test whether this * can be the directino of a new edge for a potentially new child tree. */ for(int dir=0; dir <4 ; dir++) { #ifdef OLD /* This is a more costly (slower) version of the algorithm below */ const edge tste( thisnod[n],dir) ; if ( tstedge(tste) ) { tree candi(*this) ; candi.addedge(tste) ; candi.normalize() ; bool congru=false ; for(int pres=0; pres < int(resul.size()); pres++) { if ( resul[pres].iscongru(candi) ) { congru= true ; break ; } } if ( !congru) resul.push_back(candi) ; } #else /* For this particular node and direction, we test whether * the node and the new node are a candidate for extension. We generate * a new node into the particular direction, and admit this as a new * edge if this new site is not yet already occupied. */ const node tste( thisnod[n],dir) ; /* Test on occupation of the point in the square lattice uses havnode(). */ if ( !havnode(tste) ) { /* If this node is not yet occupied (allowed since a tree is cycle-free), * we create the new child by first copying this tree, and adding * a new edge that spreads between the current node and the new one. */ tree candi(*this) ; candi.addedge(edge(thisnod[n],tste)) ; /* Adding the edge may have extended the new tree into the * quadrant to the left or below, so we shift the tree to the right * or up with normalize() to place it into the first quadrant * of the coordinate system. */ candi.normalize() ; /* This new child tree may already exist in the list of trees * generated before. We test wheter it is congruent to any fo the * child trees already accumulated so far in the result vector. */ bool congru=false ; for(int pres=0; pres < int(resul.size()); pres++) { if ( resul[pres].iscongru(candi) ) { congru= true ; break ; } } /* If no version of the new child tree is found in the * current list of children, we put it into the result list. */ if ( !congru) resul.push_back(candi) ; } #endif } } /* The result list of children that are mutualy not congruental is returned */ return resul ; } protected: private: /* The bounding box \c bbox is updated to reflect the actual * size of the tree. */ void updateBB() { /* The algorithm scans each edge and both nodes in the edge * for the maximum and minimum coordinate in x and y, and updates * the four bounding box parameters where the edge sprawls outside * the bounding box seen so fare. */ if ( edges.size() ) { bbox[2] = bbox[0] = edges[0].nodes[0].xy[0] ; bbox[3] = bbox[1] = edges[0].nodes[0].xy[1] ; } for(int e=0; e < int(edges.size()); e++) { bbox[0] = min(bbox[0],edges[e].nodes[0].xy[0]) ; bbox[0] = min(bbox[0],edges[e].nodes[1].xy[0]) ; bbox[1] = min(bbox[1],edges[e].nodes[0].xy[1]) ; // bbox[1] = min(bbox[1],edges[e].nodes[1].xy[1]) ; nodes are sorted: sufficient to look at nodes[0] for the minimum y-component bbox[2] = max(bbox[2],edges[e].nodes[0].xy[0]) ; bbox[2] = max(bbox[2],edges[e].nodes[1].xy[0]) ; //bbox[3] = max(bbox[3],edges[e].nodes[0].xy[1]) ; nodes are sorted: sufficient to look at nodes[1] for the maximum y bbox[3] = max(bbox[3],edges[e].nodes[1].xy[1]) ; } } /** Normalization means rigid translation of the tree coordinates * (of all nodes in all edges) such that the tree is tightly fitting into * the first quadrant and the minimum x and y coordinates over all nodes are zero and zero. */ void normalize() { updateBB() ; if ( bbox[0] || bbox[1] ) { const short dx = -bbox[0] ; const short dy = -bbox[1] ; for(int e=0; e < int(edges.size()); e++) edges[e].shift(dx,dy) ; bbox[0] =bbox[1] = 0 ; bbox[2] += dx ; bbox[3] += dy ; } } /* Test whether \c n is already part of the tree. * @param n the node to be found or not found in the current instance. * @return true if the node is part of any edge in the tree. */ bool havnode(const node & n) const { /* In a simple linear search strategy, all edges and both nodes of each * of them is compared with the edge \c n. */ for(int e=0; e < int(edges.size()); e++) if ( edges[e].nodes[0] == n || edges[e].nodes[1] == n) /* If we found a node in the tree that matches the candidate, * we report this as early as possible for efficiency. */ return true ; /* At this point of the program, we have fallen through the double loop over * edges and nodes and have not found any match. So the result is a 'no'. */ return false ; } #ifdef OLD /** Test wheter 'e' is an acceptable/compatible choice for a new edge */ bool tstedge(const edge & e) const { const bool haven1 = havnode( e.nodes[0]) ; const bool haven2 = havnode( e.nodes[1]) ; /* return true if have node 0 but not node 1 or vice versa */ return ( haven1 != haven2) ; } #endif /** Test an edge as a candidate for a new edge in the current tree. * @param n one node of current instance of a tree. * @param dir a lookup direction (0 to 3) for extension along a new edge. * @return true (equalent to admitted) indicating that the new site * for the partner node of \c n is not yet occupied by another node * in the current tree. */ bool tstedge(const node & n, int dir) const { return !havnode(node(n,dir)) ; } /** Generate a list of all nodes in the current tree. * @return the vector of all nodes (duplicates removed) of the current tree. */ vector<node> allnodes() const { /** This is the result vector, empty at creation */ vector<node> resul ; /* We linearly search through all edges that host two nodes each */ for(int e=0; e < int(edges.size()) ;e++) { /* The algorithm::find() operation of the STL is used to compare * the node of the current each with the set of nodes already seen * by visiting the previous edges. */ vector<node>::iterator fi=resul.begin() ; vector<node>::iterator la=resul.end() ; if ( find(fi,la,edges[e].nodes[0]) == resul.end() ) resul.push_back( edges[e].nodes[0]) ; fi=resul.begin() ; la=resul.end() ; if ( find(fi,la,edges[e].nodes[1]) == resul.end() ) resul.push_back( edges[e].nodes[1]) ; } return resul ; } } ; /** Tree comparison operator. * @param t1 the first tree, supposed to be normalized in the first quadrant. * @param t2 the second tree, supposed to be normalized in the first quadrant * @return true if the two trees overlap completely (share the same list of edges). */ bool operator== (tree &t1, tree &t2) { /* Trees of differnt edge counts are different */ if ( t1.edges.size() != t2.edges.size() ) return false ; /* Trees of different bounding boxes are also different */ else if ( t1.bbox[2] != t2.bbox[2] && t1.bbox[2] != t2.bbox[3] || t1.bbox[3] != t2.bbox[2] && t1.bbox[3] != t2.bbox[3] ) return false ; else { /* Trees of same edge count and bounding box are the same if the lists * of edges, after sorting according to the same criteria of comparison operators, * are the same. */ sort(t1.edges.begin(),t1.edges.end() ) ; sort(t2.edges.begin(),t2.edges.end() ) ; return equal(t1.edges.begin(),t1.edges.end(),t2.edges.begin()) ; } } /** A human readable ASCII output version of the tree. * A matrix of dots (empty places), oh's (nodes) and dashes (up or sideways) * is plotted in ASCII art fashion. * @param t the tree to be printed. * @param os the output stream to be used. */ ostream & operator<<(ostream & os, const tree &t) { const int xdim = t.bbox[2]+1 ; // bbox only indicates the hightes occupied member const int ydim = t.bbox[3]+1 ; char **arr = new char *[2*ydim-1] ; for(int r=0 ; r < 2*ydim-1; r++) { arr[r] = new char [2*xdim-1] ; for(int c=0; c < 2*xdim-1; c++) arr[r][c]='.' ; } for(int e=0; e < int(t.edges.size()) ; e++) { int x0= t.edges[e].nodes[0].xy[0] ; int y0= t.edges[e].nodes[0].xy[1] ; arr[2*y0][2*x0]='o' ; int x1= t.edges[e].nodes[1].xy[0] ; int y1= t.edges[e].nodes[1].xy[1] ; arr[2*y1][2*x1]='o' ; if ( x0 == x1) arr[y0+y1][2*x0]='|' ; else arr[2*y0][x0+x1]='-' ; } for(int r=0; r < 2*ydim-1 ; r++) { for(int c=0; c < 2*xdim-1; c++) os << arr[r][c] ; os << endl ; } os << endl ; for(int r=0 ; r < 2*ydim-1; r++) delete [] arr[r] ; delete [] arr ; return os ; } /** Create a new generation of child trees by generation of all in-congruental * childs of the trees in \c paren. * @param paren the existing forest of parent trees * @param verbose if true the new child trees are displayed on stdout. * @return the vector of all incongruental trees with one more edge than * those of the \c paren list. */ vector<tree> gene(vector<tree> paren, bool verbose) { /* The result vector */ vector<tree> resul ; /* New trees are generated by spawning from each of the parent trees in turn, in a loop */ for(int t=0; t < int(paren.size()) ; t++) { /* The child trees of this particular parent are collected in an intermediate * forest \c chi. */ vector<tree> chi = paren[t].childs() ; for(int c=0 ; c < chi.size() ; c++) { #if 0 cout << "--------------\n" ; cout << paren[t] << endl ; cout << "generated\n" ; cout << chi[c] << endl ; cout << "--------------\n" ; #endif /* Each of this child trees is tested for congruency relative to the * existing full set of trees in the result vector. */ bool congru=false ; for(int pres=0; pres < int(resul.size()); pres++) { if ( resul[pres].iscongru(chi[c]) ) { congru= true ; break ; } } /* If a congruent version is not yet there, we add the new child tree, * and if the verbose option was used, show and enumerate it on stdout. */ if ( !congru) { if( verbose ) { cout << resul.size() << ":\n" ; cout << chi[c] << endl ; } resul.push_back(chi[c]) ; } } } return resul ; } /** The main program. * There is only one command line option, -v, which means that * an explicit list of each generatino of trees is plotted in ASCII art * on stdout. Otherwise, only the raw counts are reported. */ int main(int argc, char *argv[]) { /* We test for the presence of the -v option */ bool verbose =false ; char oc ; while ( (oc=getopt(argc,argv,"v")) != -1 ) { switch(oc) { case 'v' : verbose = true ; break ; case '?' : cerr << "Invalid command line option " << optopt << endl ; break ; } } /* The fundamental root tree (great grand root tree) is the one from (0,0) to (1,0), * with one edge and two nodes only. */ tree root ; edge strt(node(0,0),node(1,0)) ; root.addedge(strt) ; root.bbox[0]=0 ;root.bbox[1]=0 ; root.bbox[2]=1 ;root.bbox[3]=0 ; vector<tree> lvl1 = root.childs() ; #if 0 for(int t=0 ;t < int(lvl1.size()) ; t++) { cout << t << ":\n" ; for(int e=0; e < lvl1[t].edges.size() ; e++) { const edge & ee = lvl1[t].edges[e] ; cout << "from " << ee.nodes[0].xy[0] << " " << ee.nodes[0].xy[1] << " to " ; cout << ee.nodes[1].xy[0] << " " << ee.nodes[1].xy[1] << endl ; } } #endif cout << "total " << lvl1.size() << endl ; /* The program generates the first level generation from all * children of the root tree, the second level generatio from all children of * all trees in the first level etc. Obviously one can extend this further * beyond the hard-coded maximum of 8 generations of trees. */ vector<tree> lvl2 = gene(lvl1,verbose) ; cout << "total " << lvl2.size() << endl ; vector<tree> lvl3 = gene(lvl2,verbose) ; cout << "total " << lvl3.size() << endl ; vector<tree> lvl4 = gene(lvl3,verbose) ; cout << "total " << lvl4.size() << endl ; vector<tree> lvl5 = gene(lvl4,verbose) ; cout << "total " << lvl5.size() << endl ; vector<tree> lvl6 = gene(lvl5,verbose) ; cout << "total " << lvl6.size() << endl ; vector<tree> lvl7 = gene(lvl6,verbose) ; cout << "total " << lvl7.size() << endl ; vector<tree> lvl8 = gene(lvl7,verbose) ; cout << "total " << lvl8.size() << endl ; return 0 ; }