diff --git a/Makefile b/Makefile
index a38e14e..35009e0 100644
--- a/Makefile
+++ b/Makefile
@@ -2,10 +2,9 @@
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
-SRC = *.c
 CFLAGS = -Wall -Ofast
 
 all: SOCgen SOCadmissible
 
-%: %.c
+%: src/%.c
 	gcc $(CFLAGS) -o $@ $^
diff --git a/README.md b/README.md
index af79e04..6c898d2 100644
--- a/README.md
+++ b/README.md
@@ -89,7 +89,7 @@ This integer is interpreted as a vector of bits, where each bit specifies the ab
 Since we consider simple directed graphs only (no self-loops), there are `n(n-1)` possible directed edges, where `n` is the number of nodes.
 This means that the largest number of nodes possible is limited by `n=8`.
 
-While the SOCs generated by `SOCgen` satisfy some degree-order (see function `isdegreeordered(...)` in [SOCgen.c](./SOCgen.c)), `SOCgen` does not perform graph-isomorphism tests, and may output multiple isomorphic graphs.
+While the SOCs generated by `SOCgen` satisfy some degree-order (see function `isdegreeordered(...)` in [SOCgen.c](./src/SOCgen.c)), `SOCgen` does not perform graph-isomorphism tests, and may output multiple isomorphic graphs.
 
 ## License
-[GPL-3.0-or-later](./LICENSES/GPL-3.0-or-later.txt)
\ No newline at end of file
+[GPL-3.0-or-later](./LICENSES/GPL-3.0-or-later.txt)
diff --git a/src/SOCadmissible.c b/src/SOCadmissible.c
new file mode 100644
index 0000000..424f94a
--- /dev/null
+++ b/src/SOCadmissible.c
@@ -0,0 +1,404 @@
+// SPDX-FileCopyrightText: 2023 Ämin Baumeler <amin@indyfac.ch> and Eleftherios-Ermis Tselentis <eleftheriosermis.tselentis@oeaw.ac.at>
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+/***
+ * Read graphs from  file.
+ * Input `filename': 	string
+ * Input `cnt':		pointer to int
+ * Input `maxdim':	pointer to int
+ * Output:		pointer to graph data
+ *
+ * This function sets `cnt' to the total number of graphs read from `filename',
+ * sets `maxdim' to the dimension (number of nodes) of the largest graph,
+ * and returns the graph data. Graph data is structured as follows:
+ * {dimension graph 1 | (0,0) | (0,1) | ... | (1,0) | (1,1) | ... | dimension graph 2 | ...}
+ * The ith graph is at position i*(maxdim*maxdim + 1).
+ ***/
+int* read_graphs(char* filename, int *cnt, int *maxdim) {
+	// Open file
+	FILE *fp;
+	fp = fopen(filename, "r");
+	if(fp == NULL)
+		return NULL;
+	// Count number of lines, allocate space and initiate counter
+	char ch;
+	int n = 0;
+	*cnt = 0;               // Total number of graphs
+	*maxdim = 0;            // Maximum number of vertices
+	while(!feof(fp)) {
+		ch = fgetc(fp);
+		switch(ch) {
+		case '\n':
+			(*cnt)++;
+			n = 0;
+			break;
+		case '{':
+		case '}':
+			*maxdim = (n > *maxdim) ? n : *maxdim;
+			n = 0;
+			break;
+		case '0':
+		case '1':
+			n++;
+			break;
+		case ',':
+		case ' ':
+		case 0xffffffff:
+			break;
+		default:
+			return NULL;	// Format error
+		}
+	}
+	rewind(fp);     // Return to start of file
+
+	// Allocate space for graph data
+	const int gdatasize = *maxdim * *maxdim + 1;          // Size allocation per graph
+	const int arraysize = (1+ *cnt) * gdatasize;          // Total space to allocate
+	int *graphs = (int*)malloc(sizeof(int)*arraysize);    // Allocation
+	if(!graphs) {
+		fprintf(stderr, "ERROR: Could not allocate enough memory to store the graphs.\n");
+		return NULL;
+	}
+	// Parse file
+	int col = 0;    // Column number
+	int row = 0;    // Row number
+	int i = 0;      // Current graph index
+	n = -1;         // Graph dimension
+	int j = 0;      // Running index
+	while ((ch = fgetc(fp)) != EOF) {
+		// New line encountered
+		if(ch == '\n') {
+			if(row != n) {
+				fprintf(stderr, "ERROR: File not properly formatted.\n");
+				return NULL;
+			}
+			graphs[i*gdatasize] = n;
+			i++;
+			j = 0;
+			col = 0;
+			row = 0;
+			n = -1;
+		// Matrix entry encountered
+		} else if(ch == '0' || ch == '1') {
+			int x = ch - '0';
+			graphs[i*gdatasize + 1 + j] = x;
+			col++;
+			j++;
+		// End-of-row or end-of-matrix encountered
+		} else if(ch == '}') {
+			if(n == -1) {
+				n = col;
+			} else if(col > 0 && n != col) {
+				fprintf(stderr, "ERROR: File not properly formatted.\n");
+				return NULL;
+			}
+			// End of row
+			if(col != 0) {
+				row++;
+				col = 0;
+			}
+		}
+	}
+	// Close file
+	fclose(fp);
+	return graphs;
+}
+
+/***
+ * Increase intervention counter by one.
+ * Returns 0 if we run out of interventions, and 1 otherwise.
+ ***/
+int nextintervention(int n, int *interventionslen, int *interventionidx) {
+	for(int i=0; i<n; i++) {
+		if(interventionidx[i] < interventionslen[i] - 1) {
+			interventionidx[i]++;
+			for(int j=0; j<i; j++)
+				interventionidx[j] = 0;
+			return 1;
+		}
+	}
+	return 0;
+}
+
+/***
+ * Early check
+ * With this, we might find some fixed-points entries without the need of entering any recursion.
+ * If `party' is not in the range of a parent's intervention, then `party' will definitely receive a 0.
+ * If `party' has no parents, then `party' will definitely receive a 1.
+ * If the fixed-point value cannot be inferred, return -1, otherwise return the fixed-point entry.
+ ***/
+int alphapre(int party, const int *parents, const int *parentslen, const int *interventions, const int *interventionidx, int maxn) {
+	if(parentslen[party] == 0)
+		return 1;
+	const int idatasize = (maxn+1)*(maxn+1)*2;		// Max nr of interventions/party
+	for(int pidx=0; pidx<parentslen[party]; pidx++) {
+		const int parent = parents[maxn*party + pidx];
+		const int f0 = interventions[idatasize*parent + 2*interventionidx[parent] + 0];
+		const int f1 = interventions[idatasize*parent + 2*interventionidx[parent] + 1];
+		if(f0 != party && f1 != party)
+			return 0;
+	}
+	return -1;
+}
+
+/***
+ * Recursively compute the fixed point
+ * This methods reuses already computed fixed points.
+ * `fp'    pointer to fixed-point array
+ * `fplen' the first `fplen' parties have already computed the fixed points
+ *
+ * Make sure to run alphapre before running this.
+ ***/
+int alpha(int *path, int pathlen, const int *parents, const int *parentslen, const int *interventions, const int *interventionidx, const int *fp, int n, int maxn) {
+	const int party = path[pathlen-1];
+	for(int i=0;i<pathlen-1;i++) {
+			if(party == path[i]) {
+				return 0;
+		}
+	}
+	const int idatasize = (maxn+1)*(maxn+1)*2;		// Max nr of interventions/party
+	// Iterate over all parents
+	// If one parent does not vote for `party', immediattely return a zero.
+	// If all parents vote for `party' --- which can only be inferred after we queried all parents --- return a one.
+	for(int pidx=0; pidx<parentslen[party]; pidx++) {
+		const int parent = parents[maxn*party + pidx];
+		const int f0 = interventions[idatasize*parent + 2*interventionidx[parent] + 0];
+		const int f1 = interventions[idatasize*parent + 2*interventionidx[parent] + 1];
+		if(f0 == party && f1 == party)
+			continue;
+		// Re-use already computed fixed points, else enter recursion
+		int val;
+		if(fp[parent] != -1) {
+			val = fp[parent];
+		} else {
+			path[pathlen] = parent;
+			val = alpha(path, pathlen+1, parents, parentslen, interventions, interventionidx, fp, n, maxn);
+		}
+		// Evaluate the intervention on the input
+		const int muOfVal = (val == 0) ? f0 : f1;
+		if(muOfVal != party)
+			return 0;
+	}
+	return 1;
+}
+
+/***
+ * ceil(log(x, base=10))
+ * Helper function for pretty printing.
+ ***/
+int ceillog10(int x) {
+	int res = 0;
+	while(x > 0) {
+		x /= 10;
+		res++;
+	}
+	return res;
+}
+
+/***
+ * Verify the admissibility of the graphs
+ * Go through all graphs and test them
+ * `startidx': Index from which to start, assumed to be non-negative
+ * `len':      Number of graphs to test, assumed to be non-negative and such that startidx+len does not exceed the length of the `graphs' array
+ * `maxn':     Max. number of nodes
+ * `graphs':   Pointer to graph structure
+ *
+ * Return 0 if test succeeds, else return graphidx+1 of failed graph (this is the line number).
+ * Return -1 on error.
+ ***/
+int test_graphs(int startidx, int len, int maxn, const int *graphs) {
+	const int gdatasize = maxn * maxn + 1;        // Size allocation per graph
+	const int idatasize = (maxn+1)*(maxn+1)*2;    // Max number of interventions per party
+	int range[maxn+1];                            // Range of intervention
+	int interventions[maxn*idatasize];            // Hold all interventions per party {f_0(0), f_0(1), g_0(0), g_0(1), ... , f_1(0), f_1(1), ...}
+	int interventionslen[maxn];                   // Number of interventions per party
+	int interventionidx[maxn];                    // Current intervention executed
+	int parents[maxn*maxn];                       // List all parents per party
+	int parentslen[maxn];                         // Number of parents per party
+	int children[maxn*maxn];                      // List all children per party
+	int childrenlen[maxn];                        // Numer of children per party
+	const int padlen = ceillog10(len);
+	time_t t0 = time(NULL);
+	const time_t tstart = t0;
+	for(int graphidx=startidx; graphidx<startidx+len; graphidx++) {
+		time_t t1 = time(NULL);
+		if(t1 > t0) {
+			const int count = graphidx - startidx + 1;
+			const float percentage = 100*(float)count/(float)len;
+			const int deltat = t1-tstart;
+			const float rate  = (float)count/(float)deltat;
+			fprintf(stderr, "\r%6.2f%%    %*d/%d (%*.2f graphs/s in %d seconds; current line %d)",percentage,padlen,count,len,padlen+3,rate,deltat,graphidx+1);
+			t0 = t1;
+		}
+		// Test graph with index `graphidx'
+		// `n' is the order of this graph
+		const int n = graphs[graphidx*gdatasize];
+		// Get parents and children
+		// Adj matrix is such that (row,col) = 1 <=> edge from row to col
+		for(int i=0; i<n; i++) {
+			parentslen[i] = 0;
+			childrenlen[i] = 0;
+		}
+		int idx = 0;
+		for(int row=0; row<n; row++) {
+			for(int col=0; col<n; col++) {
+				if(graphs[graphidx*gdatasize+1+idx] == 1) {
+					// Edge from row to col
+					children[maxn*row + childrenlen[row]] = col;
+					childrenlen[row]++;
+					parents[maxn*col + parentslen[col]] = row;
+					parentslen[col]++;
+				}
+				idx++;
+			}
+		}
+		// Fill in interventions
+		for(int party=0; party<n; party++) {
+			// This is a small trick to speed-up processing:
+			// If the party `party' has two or more children, ignore the `discard' intervention (-1) (tmp=0),
+			// else make use of the `discard' internvetion (tmp=1).
+			// Populate range of `party's interventions
+			// The range has childrenlen[party]+tmp entries
+			range[0] = -1;
+			const int tmp = (childrenlen[party]<=1) ? 1 : 0;
+			for(int j=0; j<childrenlen[party]; j++) {
+				range[j+tmp] = children[maxn*party+j];
+			}
+			// Enumerate, store, and count all possible interventions
+			// Each intervention is a tuple (mu(0), mu(1)), x corresponds to the first entry, y to the second
+			int cnt = 0;
+			for(int x=0; x<childrenlen[party]+tmp; x++) {
+				for(int y=0; y<childrenlen[party]+tmp; y++) {
+					interventions[party*idatasize+2*cnt] = range[x];
+					interventions[party*idatasize+2*cnt+1] = range[y];
+					cnt++;
+				}
+				interventionslen[party] = cnt;
+			}
+		}
+		// Loop over intervetions and verify fixed-point expression
+		int fp[n];                                       // Array to store alpha(party) values
+		int afterfp[n];                                  // \omega(\mu(fp))
+		int path[n+1];                                   // Current path
+		for(int i=0; i<n; i++) interventionidx[i] = 0;	 // Initialize intervention index to the first intervention
+		int running = 1;  // Flag to see whether we are still verifying, or we verified all interventions
+		while(running) {
+			// Reset fixed-point entries
+			for(int party=0; party<n; party++)
+				fp[party] = -1;
+			// Fill out fixed-points that are trivially computed
+			for(int party=0; party<n; party++)
+				alphapre(party, parents, parentslen, interventions, interventionidx, maxn);
+			// Invoke recursive function
+			for(int party=0; party<n; party++) {
+				path[0] = party;
+				const int val = alpha(path, 1, parents, parentslen, interventions, interventionidx, fp, n, maxn);
+				fp[party] = val;
+			}
+
+			// Compute \omega(\mu(fp))
+			for(int i=0; i<n; i++)
+				afterfp[i] = -1;
+			// Loop over parties
+			for(int party=0; party<n; party++) {
+				// Loop over parents of this party
+				for(int pidx=0; pidx<parentslen[party]; pidx++) {
+					const int parent = parents[maxn*party + pidx];
+					const int f0 = interventions[idatasize*parent + 2*interventionidx[parent] + 0];
+					const int f1 = interventions[idatasize*parent + 2*interventionidx[parent] + 1];
+					// If `parent' never votes for `party', this value is 0.
+					if(f0 != party && f1 != party) {
+						afterfp[party] = 0;
+						break;
+					}
+					// If `parent' always votes for `party', check the other parents.
+					if(f0 == party && f1 == party)
+						continue;
+					// Evaluate intervention on the input (fp[parent]) 
+					const int muOfVal = (fp[parent] == 0) ? f0 : f1;
+					// If `parent' does not vote for `party', this value is 0, else check the other parents.
+					if(muOfVal != party) {
+						afterfp[party] = 0;
+						break;
+					}
+
+				}
+				// All parents voted `party'
+				if(afterfp[party] == -1)
+					afterfp[party] = 1;
+			}
+			// Compare the output of the recursive alpha function (fp) with \omega(\mu(fp)) (afterfp)
+			// If these two arrays differ, return the line number of the graph we falsified
+			for(int party=0;party<n;party++)
+	 			if(fp[party] != afterfp[party])
+					return graphidx+1;
+
+			// Next intervention
+			running = nextintervention(n, interventionslen, interventionidx);
+		}
+	}
+	const int deltat = time(NULL) - tstart;
+	const float rate  = (deltat == 0) ? len : (float)len/(float)deltat;
+	fprintf(stderr, "\r%6.2f%%    %*d/%d (%*.2f graphs/s in %d seconds; current line %d)\n",100.0,padlen,len,len,padlen+3,rate,deltat,startidx+len);
+	return 0;
+}
+
+int main(int argc, char *argv[]) {
+	// Parse command-line arguments
+	int START = 0;
+	int NUM = 0;
+	if(argc >= 3)
+		START = atoi(argv[2]);
+	if(argc >= 4) {
+		if(argv[3][0] == '+')
+			NUM = atoi(argv[3]+1);
+		else
+			NUM = atoi(argv[3])-START+1;
+	}
+	if(!(argc >= 2 && (argc < 3 || START > 0) && (argc < 4 || NUM > 0) && argc <= 4)) {
+		fprintf(stderr, "Usage: %s <filename> [<startline> [<endline> | +<count>]]\n", argv[0]);
+		fprintf(stderr, "       <filename>     File name with adjacency matrices of simple directed graphs\n");
+		fprintf(stderr, "       <startline>    Verify graphs starting from line `startline'\n");
+		fprintf(stderr, "       <endline>      Verify graphs up to and including line `endline'\n");
+		fprintf(stderr, "       +<count>       Verify `count' number of graphs\n");
+		fprintf(stderr, "\n");
+		fprintf(stderr, "[FILE FORMAT]\n");
+		fprintf(stderr, "       Each line in `filename' must contain the adjacency matrix of a simple directed graph in the format\n");
+		fprintf(stderr, "         {{a00,a01,...},{a10,a11,...},...}  where aij=1 if and only if the graph has the edge i -> j\n");
+		fprintf(stderr, "       The file `filename' may contain graphs with different order (number of vertices)\n");
+		fprintf(stderr, "\n");
+		fprintf(stderr, "This program verifies the admissibility of simple directed graphs.\n");
+		return -1;
+	}
+
+	// Read Graphs from File
+	int graphcount = 0;
+	int maxdimension = -1;
+	const int *graphs = read_graphs(argv[1], &graphcount, &maxdimension);
+	if(graphs == NULL) {
+		fprintf(stderr, "ERROR: Could not read graphs from file `%s'.\n", argv[1]);
+		return -1;
+	}
+	const int startidx = START ? START-1 : 0;		// Index of graph to start verifying
+	const int tot = NUM ? NUM : graphcount-startidx;	// Number of graphs to be verified
+	if(startidx >= graphcount || tot <= 0 || tot+startidx-1 >= graphcount) {
+		fprintf(stderr, "ERROR: Lines to be verified are out of range\n");
+		return -1;
+	}
+
+	// Test `tot' graphs, starting from index `startidx'
+	printf("Verifying the admissibility of %d graphs in the file `%s' (line %d to line %d)\n", tot, argv[1], startidx+1, startidx+tot);
+	const int result = test_graphs(startidx, tot, maxdimension, graphs);
+	if (result == 0)
+		printf("These graphs are admissible\n");
+	else if(result == -1)
+		fprintf(stderr, "ERROR: Something went wrong\n");
+	else 
+		printf("The function alpha does not represent the fixed point, or the graph on line %d is inadmissible\n", result);
+	return result;
+}
diff --git a/src/SOCgen.c b/src/SOCgen.c
new file mode 100644
index 0000000..6e53bcc
--- /dev/null
+++ b/src/SOCgen.c
@@ -0,0 +1,519 @@
+// SPDX-FileCopyrightText: 2023 Ämin Baumeler <amin@indyfac.ch> and Eleftherios-Ermis Tselentis <eleftheriosermis.tselentis@oeaw.ac.at>
+//
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+
+/***
+ * Print the graph's adjacency matrix
+ ***/
+void dumpgraph(int n, const int *children, const int *childrenlen) {
+	int adj[n*n];
+	for(int i=0; i<n*n; i++) adj[i] = 0;
+	for(int a=0; a<n; a++) {
+		for(int bidx=0; bidx<childrenlen[a]; bidx++) {
+			const int b = children[a*n+bidx];
+			adj[a*n+b] = 1;
+		}
+	}
+	// Print
+	printf("{");
+	for(int a=0; a<n; a++) {
+		if(a) printf(",");
+		printf("{");
+		printf("%d", adj[a*n+0]);
+		for(int b=1; b<n; b++) printf(",%d", adj[a*n+b]);
+		printf("}");
+	}
+	printf("}\n");
+}
+
+/***
+ * Block size of our datastructure.
+ ***/
+int blocksize(int n) {
+	int res=(n+1)*n;
+	for(int i=0; i<n; i++)
+		res *= n-i;
+	return res;
+}
+
+/***
+ * Recrusive bi-directional graph traversal, used to check connectivity.
+ ***/
+void traverse(int vertex, int *visited, int n, const int *parents, const int *parentslen, const int *children, const int *childrenlen) {
+	visited[vertex] = 1;
+	for(int i=0; i<parentslen[vertex]; i++) {
+		int u = parents[n*vertex + i];
+		if(visited[u] == 0)
+			traverse(u, visited, n, parents, parentslen, children, childrenlen);
+	}
+	for(int i=0; i<childrenlen[vertex]; i++) {
+		int u = children[n*vertex + i];
+		if(visited[u] == 0)
+			traverse(u, visited, n, parents, parentslen, children, childrenlen);
+	}
+}
+
+/***
+ * Return 1 if the graph is connected, and 0 otherwise.
+ ***/
+int isconnected(int n, const int *parents, const int *parentslen, const int *children, const int *childrenlen) {
+	// Early check: Is there a node with 0 total degree?
+	for(int i=0; i<n; i++) if(parentslen[i]+childrenlen[i]==0) return 0;
+	// Check connectivity
+	int visited[n];
+	for(int i=0; i<n; i++) visited[i] = 0;
+	// Start traversal
+	traverse(0, visited, n, parents, parentslen, children, childrenlen);
+	for(int i=0; i<n; i++) if(!visited[i]) return 0;
+	return 1;
+}
+
+/***
+ * Interpret graph number `graph' as a simple digraph, and populate the arrays `parents', `parentslen', `children', `childrenlen'
+ * The number `graph' is understood as the adjacency matrix with the diagonal removed (no self loops).
+ ***/
+void graphnrtolists(int n, unsigned long long int graph, int *parents, int *parentslen, int *children, int *childrenlen) {
+	for(int i=0; i<n; i++) {
+		parentslen[i] = 0;
+		childrenlen[i] = 0;
+	}
+	int row = 0;
+	int col = 1;    // We cannot start at (0,0); it does not exist
+	unsigned long long g = graph;
+	for(int i=0; i<n*(n-1); i++) {
+		int val = g%2;
+		g = g/2;
+		if(val) {
+			children[row*n+childrenlen[row]] = col;
+			childrenlen[row]++;
+			parents[col*n+parentslen[col]] = row;
+			parentslen[col]++;
+		}
+		col++;
+		if(col==row)
+			col++;
+		if(col==n) {
+			col = 0;
+			row++;
+		}
+
+	}
+}
+
+/***
+ * Find the tail in `path' of length `pathlen' that includes a cycle, and save the canonical representation.
+ * The canonical representation is the listing of the nodes of the cycle starting from the smallest one.
+ ***/
+void canonical_cycle(int n, int *cycle, int *cyclelen, const int *path, int pathlen) {
+	int startidx;
+	int pathhascycle = 0;
+	const int last = path[pathlen-1];                       // This is the last node
+	for(int k=0; k<pathlen-1; k++) {
+		if(path[k]==last) {             // We have found the first occurance of node `last' on the path
+			// The cycle the following part of path: path[startidx], path[startidx+1], ..., path[pathlen-2], and therefore has length pathlen-startidx-1
+			startidx = k;
+			(*cyclelen) = pathlen-startidx-1;
+			pathhascycle = 1;
+			break;
+		}
+	}
+	if(!pathhascycle) {                     // Sanity check. This should NEVER happen!
+		fprintf(stderr, "ERROR: You provided a path to `canonical_cycle' that does NOT contain a cycle as its tail! I'm returning now.\n");
+		return;
+	}
+	// Generate canonical representation
+	// First, find smallest node on cycle
+	int smallestidx = -1;
+	int smallest = n;
+	for(int k=startidx; k<pathlen-1; k++) {
+		if(smallest > path[k]) {
+			smallest = path[k];
+			smallestidx = k;
+		}
+	}
+	// Now, copy cycle to `cycle' starting from the smallest
+	int j=0;                                        // This is the index on where to write on the array `cycle'
+	for(int i=smallestidx;; i++) {
+		if(j==*cyclelen) break;                 // Stop filling in, if we reached the end
+		if(i==pathlen-1) i = startidx;          // If we reach the end of the path, jump to the start
+		cycle[j] = path[i];                     // Write at position `j' the element from the cycle on the path at position `i`'
+		j++;                                    // Increase index to write on `cycle'
+	}
+}
+
+/***
+ * Return 1 if the cycle is not in `cycles', otherwise 0.
+ ***/
+int cycleisnew(int n, const int *c, int clen, const int *cycles, const int *cyclescnt) {
+	const int cnt = cyclescnt[clen];
+	const int bs = blocksize(n);
+	for(int i=0; i<cnt; i++) {                      // Iterate over all cycles of the same length
+		int j;
+		const int startidx = clen*bs + i*clen;  // We start here
+		for(j=0; j<clen; j++)                   // Go position by position
+			if(cycles[startidx+j] != c[j])
+				break;                  // They differ at position j
+		if(j==clen)                             // We didn't break
+			return 0;
+	}
+	return 1;
+}
+
+/***
+ * Save the given cycle `c' of length `clen' in `cycles'.
+ ***/
+void savecycle(int n, const int *c, int clen, int *cycles, int *cyclescnt) {
+	const int m = cyclescnt[clen];                          // There are already `m' cycles of length `clen'
+	const int bs = blocksize(n);
+	const int startidx = clen*bs + m*clen;                  // We start writing here
+	for(int i=0; i<clen; i++) cycles[startidx+i] = c[i];
+	cyclescnt[clen]++;                                      // Increase the counter that holds the number of cycles of length `clen'
+}
+
+
+/***
+ * Recursive depth-first-search to find cycles in the graph specified by `children' and `childrenlen'.
+ * The cycles are stored in `cycles' and `cyclescnt'.
+ * Parameter `useless':     n-array where useless nodes are marked
+ * Parameter `path':        Array of nodes that represent a path on the graph
+ * Parameter `pathlen':     Lenght of `path'
+ * Parameter `visited':     n-array to mark visited nodes
+ * Parameter `n':           Number of nodes
+ * Parameter `found':	Pointer to integer to count the number of cycles found
+ ***/
+void dfs(const int *useless, int *path, int pathlen, int *visited, int n, const int *children, const int *childrenlen, int *cycles, int *cyclescnt, unsigned int *found) {
+	int cycle[n];                           // Temporary store for cycle
+	int cyclelen;                           // Length of temporarilly stored cycle
+	const int vertex = path[pathlen-1];     // We go on duing a dfs from node `vertex'
+	for(int i=0; i<childrenlen[vertex]; i++) {
+		const int u = children[n*vertex + i];           // Enter node `u'
+		if(useless[u]) continue;                        // This one is useless
+		if(visited[u] == 0) {                           // Not visited yet
+			visited[u] = 1;                         // Mark as visited
+			path[pathlen] = u;                      // Append to path, and re-enter recursion...
+			dfs(useless, path, pathlen+1, visited, n, children, childrenlen, cycles, cyclescnt, found);
+			visited[u] = 0;                         // Undo marking, we are going to take another route
+		} else {
+			// We have reached a node on the path, i.e., we have found a cycle!
+			// The cycle is somewhere along the path
+			// path[0], path[1], ... path[pathlen-1]
+			// More specifically, node `u' is somewhere there.
+			// To detect this, we first place `u' on the path.
+			// The method `canonical_cycle' will then use this information to find and return the canonical form of the cycle.
+			path[pathlen] = u;
+			canonical_cycle(n, cycle, &cyclelen, path, pathlen+1);          // Find cycle and generate the canonical representation.
+			if(cycleisnew(n, cycle, cyclelen, cycles, cyclescnt)) {         // Is it new? If yes, save
+				(*found)++;                                             // Increase the number of found cycles by 1
+				savecycle(n, cycle, cyclelen, cycles, cyclescnt);       // Save cycle
+			}
+			// else dont do anything..
+		}
+	}
+	return;
+}
+
+/***
+ * Find cycles in the `n'-nodes graph specified by the arrays `children' and `childrenlen', and return the number of cycles found.
+ * The cycles are stored in `cycles' and `cyclescnt'.
+ * This uses a recursive depth-first-search.
+ ***/
+int find_cycles(int *cycles, int *cyclescnt, int n, const int *children, const int *childrenlen) {
+	int visited[n];
+	int useless[n];                 // Some nodes are useless, mark them as such
+	int path[n+1];                  // We use n+1 here because we will "close" the path to form a cycle, i.e., one node will be repeated.
+	unsigned int found = 0;         // Store total number of cycles found
+	// Initialize cyclescnt to zero
+	for(int i=0; i<n+1; i++)
+		cyclescnt[i] = 0;
+	// Mark nodes without children as useless
+	for(int i=0; i<n; i++) {
+		useless[i] = (childrenlen[i]==0) ? 1 : 0;
+	}
+	// Iterate over nodes
+	for(int vertex=0; vertex<n; vertex++) {
+		if(useless[vertex]) continue;
+		for(int i=0; i<n; i++) visited[i] = 0;  // Mark all as univisited
+		path[0] = vertex;                       // We start at `vertex'
+		visited[vertex] = 1;                    // Mark `vertex' as visited. If we revisit a node, we found a cycle.
+		// Enter depth-first-search recursion
+		dfs(useless, path, 1, visited, n, children, childrenlen, cycles, cyclescnt, &found);
+		// We processed this. Every cycle that goes through `vertex' has been found.
+		useless[vertex] = 1;
+	}
+	return found;
+}
+
+/***
+ * gissoc tests wheter the graph provided is a SOC or not.
+ * Parameter `n':		Number of nodes
+ * Parameter `parents':	Pointer to list of parents per node
+ * Parameter `parentslen':	Pointer to list of number of parents per node
+ * Parameter `children':	Pointer to list of children per node
+ * Parameter `childrenlen':	Pointer to list of number of children per node
+ * Parameter `cycles':          Pointer to list of cycles in the graph
+ * Parameter `cyclescnt':       Pointer to list of cycles count
+ *
+ * This function retruns 0 if the graph is NOT a SOC, and 1 otherwise.
+ *
+ * A SOC (Siblings-on-Cycles) is a simple digraph where every cycle has siblings, i.e., nodes with common parents.
+ ***/
+int gissoc(int n, const int *parents, const int *parentslen, const int *children, const int *childrenlen, int *cycles, int *cyclescnt) {
+	// We store the cycles in groups of the cycle lenghts (length of cycle is the number of edges, or equivalentelly, the number of distinct nodes)
+	// Since the graph is a simple digraph, the smallest cycles have length 2.
+	// cycles[0], cycles[1] is the first cycle of length 2
+	// cycles[2], cycles[3] is the second cycle of length 2
+	// etc.
+	// There are at most (n choose 2) cycles of length 2.
+	// There are at most (n choose n/2) cycles of length n/2. This is upper bounded by 2^n/sqrt(pi*n/2) < 2^n; this is the largest binomial coefficient.
+	// Since generally n will be small, we just use this bound for all groups.
+	// So, the cycles with length k are saved at and after cycles[k*(2^n)].
+	// Each block has size 2^n.
+	/* Early and saftey check for self loops. SOCs dont have self-loops, and we should never generate them. Anyhow, test */
+	for(int i=0; i<n; i++) {
+		int num = childrenlen[i];
+		for(int cidx=0; cidx<num; cidx++)
+			if(i == children[n*i+cidx]) {
+				fprintf(stderr, "ERROR: THIS GRAPH HAS A SELF LOOP (node %d)\n", i);
+				return 0;
+			}
+	}
+	// Check SOC property
+	const int bs = blocksize(n);
+	// Iterate over cycle length; min length is 2, max length is n
+	// then iterate over cycle of length `len'
+	// then iterate over nodes of that cycle.
+	// In that iteration process, mark all parents.
+	// If a marked parent is marked again, we have found a common parent.
+	for(int len=2; len<n+1; len++) {
+		const int cnt = cyclescnt[len];             // Number of cycles of length `len'
+		for(int i=0; i<cnt; i++) {                  // Iterate over cycles
+			const int startidx = len*bs+len*i;  // The cycle we look at (length=`len', index=`i') is saved starting from index `startidx'
+			// Initialize marks
+			int marked[n];
+			for(int k=0; k<n; k++) marked[k] = 0;
+			// Flag: Have we found a common parent?
+			int commonparentfound = 0;
+			// Iterate over cycle
+			for(int k=0; k<len; k++) {
+				const int vertex = cycles[startidx+k];
+				for(int pidx=0; pidx<parentslen[vertex]; pidx++) {
+					const int parent = parents[n*vertex+pidx];
+					if(marked[parent]) commonparentfound = 1;    // The node `parent' is a common parent
+					else marked[parent] = 1;                     // Mark the parent node
+					if(commonparentfound) break;                 // We have found a common parent, break out of the for loop over `pidx'
+				}
+				if(commonparentfound) break;                         // We have found a common parent, break out of the for loop over `k'
+			}
+			if(!commonparentfound)                          // On cycle `i' we did NOT find a common parent, this is not a SOC
+				return 0;
+		}
+	}
+	return 1;                                                       // If all test pass, we have found a SOC
+}
+
+/***
+ * Number of grandparents, and number of grandchildren
+ ***/
+int grandparentslen(int node, const int *parents, const int *parentslen) {
+	int tot = 0;
+	for(int i=0; i<parentslen[node]; i++)
+		tot += parentslen[i];
+	return tot;
+}
+int grandchildrenlen(int node, const int *children, const int *childrenlen) {
+	int tot = 0;
+	for(int i=0; i<childrenlen[node]; i++)
+		tot += childrenlen[i];
+	return tot;
+}
+/***
+ * Speed-up by ignoring all graphs that do not satisfy some ordering of out- and in-degrees.
+ * If the graph does not satisfy this `canonical' form, then another ismorphic to it (it's simply a permutation of the nodes)
+ * will be degree-ordered.
+ ***/
+int isdegreeordered(int n, const int *parents, const int *parentslen, const int *children, const int *childrenlen) {
+	// Decreasingling order by in-degree, then out-degree, then ...
+	for(int i=0; i<n-1; i++) {
+		if(parentslen[i] > parentslen[i+1]) {
+			return 0;
+		} else if(parentslen[i] == parentslen[i+1]) {
+			if(childrenlen[i] > childrenlen[i+1]) {
+				return 0;
+			} else if(childrenlen[i] == childrenlen[i+1]) {
+				if(grandparentslen(i,parents,parentslen) > grandparentslen(i+1,parents,parentslen)) {
+					return 0;
+				} else if(grandparentslen(i,parents,parentslen) == grandparentslen(i+1,parents,parentslen)) {
+					if(grandchildrenlen(i,children,childrenlen) > grandchildrenlen(i+1,children,childrenlen)) {
+						return 0;
+					}
+				}
+			}
+		}
+	}
+	return 1;
+}
+
+/* Random nubers */
+#if RAND_MAX/256 >= 0xFFFFFFFFFFFFFF
+#define LOOP_COUNT 1
+#elif RAND_MAX/256 >= 0xFFFFFF
+#define LOOP_COUNT 2
+#elif RAND_MAX/256 >= 0x3FFFF
+#define LOOP_COUNT 3
+#elif RAND_MAX/256 >= 0x1FF
+#define LOOP_COUNT 4
+#else
+#define LOOP_COUNT 5
+#endif
+
+u_int64_t rand_uint64(unsigned long long max) {
+	u_int64_t r = 0;
+	for (int i=LOOP_COUNT; i > 0; i--) {
+		r = r*(RAND_MAX + (u_int64_t)1) + rand();
+	}
+	return r%max;
+}
+/* End of random numbers */
+
+int main(int argc, char *argv[]) {
+	// Parse command-line arguments
+	int n = -1;
+	int NONDAGONLY = 0;
+	int NOSINK = 0;
+	int NOSOURCE = 0;
+	int RANDOM = 0;
+	int UNKOPTION = 0;
+	for(int i=1; i<argc; i++) {
+		if(strcmp(argv[i], "-n") == 0 && i+1 < argc) {
+			n = atoi(argv[i+1]);
+			i++;
+		} else if(strcmp(argv[i], "-r") == 0 && i+1 < argc) {
+			RANDOM = atoi(argv[i+1]);
+			i++;
+		} else if(strcmp(argv[i], "-c") == 0)
+			NONDAGONLY = 1;
+		else if(strcmp(argv[i], "--no-sink") == 0)
+			NOSINK = 1;
+		else if(strcmp(argv[i], "--no-source") == 0)
+			NOSOURCE = 1;
+		else {
+			UNKOPTION = 1;
+			break;
+		}
+	}
+	if(UNKOPTION || n<=1) {
+		fprintf(stderr, "Usage: %s -n <order> [-r <num> ] [FILTER ...]\n", argv[0]);
+		fprintf(stderr, "       -n <order>     Generate SOCs with `order' connected nodes\n");
+		fprintf(stderr, "       -r <num>       Pick directed graphs at random, and exit after having found `num' SOCs\n");
+		fprintf(stderr, "\n");
+		fprintf(stderr, "[FILTER]              Consider only simple directed graphs ...\n");
+		fprintf(stderr, "       -c             ... that are cyclic (i.e., not DAGs)\n");
+		fprintf(stderr, "       --no-sink      ... without sink nodes (this logically implies -c)\n");
+		fprintf(stderr, "       --no-source    ... without source nodes (also this logically implies -c)\n");
+		fprintf(stderr, "\n");
+		fprintf(stderr, "This program prints the found SOCs as adjacency matrices to stdout.\n");
+		fprintf(stderr, "To exclude (some) of the isomorphic SOCs, it uses a degree-order filter.\n");
+		return -1;
+	}
+	// EO Parse command-line arguments
+	// Setup datastructures
+	int parents[n*n];       // Each node i can have at most n-1 parents. These are listed at parents[n*i+k] for 0<=k<n
+	int parentslen[n];      // Number of parents of node i is stored in parentslen[i]
+	int children[n*n];      // Each node i can have at most n-1 children. These are listed at parents[n*i+k] for 0<=k<n
+	int childrenlen[n];     // Number of children of node i is stored in childrenlen[i]
+	const int bs = blocksize(n);
+	int *cycles = (int*)malloc(sizeof(int)*(n+2)*bs);
+	if(!cycles) {
+		fprintf(stderr, "Failed to allocate memory to store cycles (tried to allocate %lu bytes)\n", sizeof(int)*(n+2)*bs);
+		return -1;
+	}
+	int cyclescnt[n+1];      // cyclescnt[k] stores the number of cycles with length k
+	// EO Setup datastructures
+	// Initiate graph enumeration
+	const int m = n*(n-1);
+	if(m>64) {
+		fprintf(stderr, "Too many graphs to enumarate with an unsiged long long integer (number of node pairs = %d)\n", m);
+		return -1;
+	}
+	const unsigned long long max = 1L << m;		// Largest `graphnumber'
+	const int padlen = (int)((float)m/3.322)+1;	// Convert log 2 to log 10
+	int len = 0;
+	time_t t0 = time(NULL);
+	time_t t = time(NULL);
+	srand((unsigned) time(&t));
+	// EO Initiate graph enumeration
+	fprintf(stderr, "Generating SOCs with %d nodes\n", n);
+	if(RANDOM) fprintf(stderr, "  Picking graphs at random\n");
+	if(NONDAGONLY) fprintf(stderr, "  Filter: Omitting DAGs\n");
+	if(NOSINK) fprintf(stderr, "  Filter: Omitting graphs with sink nodes\n");
+	if(NOSOURCE) fprintf(stderr, "  Filter: Omitting graphs with source nodes\n");
+	unsigned long long graphnumber = 0;	// We will always omit the graph with `graphnumber' 0; it is unintersting anyway
+	unsigned long long graphschecked = 0;
+	for(;graphnumber < max && (!RANDOM || len < RANDOM);) {
+		const time_t now = time(NULL);
+		// Break if we exhaustively check all graphs or if we found enough SOCs at random
+		//if(graphnumber >= max || (RANDOM && len >= RANDOM))
+		//	break;
+		// Print status every second
+		if(now > t) {
+			t = now;
+			const int deltat = now - t0;
+			const float rate = (float)graphschecked/(float)deltat;
+			const float SOCrate = (float)len/(float)deltat;
+			const float percentage = (RANDOM==0) ? 100*(float)graphnumber/(float)max : 100*(float)len/(float)RANDOM;
+			const float ETC = (RANDOM==0) ? ((float)(max - graphschecked)/(rate))/3600 : ((float)(RANDOM - len)/(SOCrate))/3600;
+			fprintf(stderr, "\r%6.2f%%    %*llu/%llu (%i SOCs found, %d seconds, %*.2f graphs/s, %*.2f SOCs/s, %*.2fh estimated time left)",percentage,padlen,graphnumber,max,len,deltat,padlen+3,rate,padlen+3,SOCrate,padlen,ETC);
+			fflush(stderr);
+		}
+		// Convert graph index `grpahnumber' to the lists parents, children, parentslen, childrenlen
+		graphnrtolists(n, graphnumber, parents, parentslen, children, childrenlen);
+		// Increase checked counter and prepare for next iteration
+		graphschecked++;
+		if(RANDOM)
+			graphnumber = rand_uint64(max);
+		else
+			graphnumber++;
+		// EO Increase checked counter and prepare for next iteration
+		// If enabled, compare against degree options
+		if(NOSINK || NOSOURCE) {
+			// Get some degree properties
+			int minindegree = n;
+			int minoutdegree = n;
+			for(int node=0; node<n; node++) {
+				if(minindegree > parentslen[node])
+					minindegree = parentslen[node];
+				if(minoutdegree > childrenlen[node])
+					minoutdegree = childrenlen[node];
+			}
+			if(NOSINK && minoutdegree == 0) continue;
+			if(NOSOURCE && minindegree == 0) continue;
+		}
+		// Check whether this graph is canonical or not; continue if graph has not proper ordering of degress (there is an isomorphic graphs to this one that has been or will be checked)
+		if(!isdegreeordered(n, parents, parentslen, children, childrenlen)) continue;
+		// Ignore graphs that are not connected
+		if(!isconnected(n, parents, parentslen, children, childrenlen)) continue;
+		// Find cycles
+		const int num_cycles = find_cycles(cycles, cyclescnt, n, children, childrenlen);
+		// If enabled, ignore DAGs
+		if(NONDAGONLY && num_cycles > 0) continue;
+		// Test the Siblings-On-Cycles property
+		// A DAG is trivially a SOC
+		if(num_cycles > 0 && !gissoc(n, parents, parentslen, children, childrenlen, cycles, cyclescnt)) continue;
+		// We have found a SOC
+		len++;
+		dumpgraph(n, children, childrenlen);
+	}
+	const int deltat = time(NULL) - t0;
+	const float rate = (deltat == 0) ? graphschecked : (float)graphschecked/(float)deltat;
+	const float SOCrate = (deltat == 0) ? len : (float)len/(float)deltat;
+	fprintf(stderr, "\r%6.2f%%    %*llu/%llu (%i SOCs found, %d seconds, %*.2f graphs/s, %*.2f SOCs/s, %*.2fh estimated time left)",100.00,padlen,graphnumber,max,len,deltat,padlen+3,rate,padlen+3,SOCrate,padlen,0.00);
+	fprintf(stderr, "\nFound %d SOCs in %d seconds\n", len,deltat);
+	// Free manually allocated memory
+	free(cycles);
+	return 0;
+}