/*

 A PROGRAM THAT COUNTS ROOTED AND UNROOTED MAPS BY GENUS, NUMBER OF EDGES AND NUMBER OF VERTICES
 
 Timothy R. Walsh
 
 A map is a 2-cell imbedding of a connected graph, loops and multiple edges allowed, in a surface.  
 This program only counts maps on compact, orientable surfaces with no boundary,
 which are characterized by a single non-negative integer, their genus.  
 A surface of genus g is like a doughnut with g holes.  
 A map is rooted if a half-edge is distinguished as the root.  
 Counting maps means counting equivalence classes of maps under orientation-preserving homeomorphism.  
 In the case of rooted maps, the homeomorphism preserves the root and, therefore, 
 preserves all the edges, vertices and faces, making rooted maps easier to count.  
 The method used to count unrooted maps requires the numbers of rooted maps; so these must be counted first.
 
 The method used to count rooted maps was very recently published in 
 [S. R. Carrell and G. Chapuy. Simple recurrence formulas to count maps on orientable surfaces. 
 ArXiv e-prints, 2014. Submitted on 19 March. URL http://arxiv.org/abs/1402.6300]. 
 It is vastly more efficient than the previously best method, published in [T. R. S. Walsh and A. Giorgetti. 
 Efficient enumeration of rooted maps of a given orientable genus by number of faces and vertices. 
 Ars Mathematica Contemporanea, 7:263-280, 2014], which was limited to genus 10.  
 The method used to count unrooted maps is a refinement published in 
 [V.A. Liskovets. A multivariate arithmetic function of combinatorial and topological significance. 
 Integers, 10:155-177, 2010] of the one published in 
 [A. Mednykh and R. Nedela. Enumeration of unrooted maps of a given genus. 
 J. Comb. Theory, Ser. B, 96(5):706-729, 2006].  
 In this last article, unrooted maps were counted by genus and number of edges alone.  
 The way of adding the number of vertices as a parameter was published in 
 [T. R. S. Walsh, A. Giorgetti, and A. Mednykh. 
 Enumeration of unrooted orientable maps of arbitrary genus by number of edges and vertices. 
 Discrete Mathematics, 312(17):2660-2671, 2012]; a more pedagogical exposition can be found in 
 [T. R. S. Walsh. Counting maps on doughnuts. Theor. Comput. Sci., 502:4-15, 2013].
 
 The program that follows this documentation counts rooted and unrooted maps by 
 orientable genus, number of edges and number of vertices.  
 It is written in C++ and needs two software packages to run it: 
 CLN to handle large integers and the C/C++ compiler XCODE to run CLN.  
 Both of these packages can be found on the internet - Google CLN and XCODE to find the web sites - 
 and downloaded for free, but they need to be installed on your computer 
 and you may need some help installing them.  
 Once you have them up and running, choose Command Line Tools C++ under New Project, 
 copy the program below and paste it into the window that contains the words "Hello World".  
 Press "build and run" to compile and run the program.
 
 My colleague Alain Giorgetti informed me that the program can also be compiled with g++ (under linex 64 bits).
 
 It could be that your computer's memory is too small to hold the large arrays in this program, 
 in which case you will see a message to that effect when you try to compile the program.  
 You can adjust the size of the arrays by changing the first of two defined constants with numerical values: 
 edgemax, which is currently defined as 100.  
 This number is the greatest number of edges of the maps that this program can count.  
 You can increase or decrease this number by changing it to another one.  
 The amount of memory that the arrays will occupy varies as the cube of this number.  
 The other constant with a numerical value is maxdiv, the last of the defined constants, 
 which is currently set to 20.  
 This is 1 greater than the maximum number of positive divisors of any integer up to and including 4*genusmax+2, 
 where genusmax is half of edgemax, with fractions rounded down.  
 The current value of maxdiv is big enough as long as edgemax is less than 180.  
 If you want to know the smallest value of maxdiv that will work for the value of edgemax that you choose, 
 look up highly composite numbers on Google.
 
 Once your program is up and running, you will see a few lines on your screen telling you what you can do, 
 but I'll repeat this information here.
 
 You will be prompted to enter the maximum number of edges in the maps you want to count - 
 any non-negative integer up to and including edgemax, or -1 to quit.  
 If you enter a larger integer, it will be changed to edgemax.  
 If you enter a negative integer, the program will terminate execution - 
 this is one way to terminate the program (the other one is to press the red octagon), 
 which will otherwise let you repeat the enumeration with different choices as often as you want.  
 Then you will be prompted to enter the maximum genus of the maps you want to count - 
 any non-negative integer up to and including half of the greatest number of edges you chose, 
 with fractions rounded down.  
 If you enter a greater integer than that, it will be changed to the greatest one you were allowed to enter; 
 if you enter a negative integer, it will be changed to zero.  
 Then you will see the words "counting rooted maps" 
 and be prompted to enter 1 to store the numbers in a file called 'rootx.txt' or 0 not to.  Whichever you choose, 
 the computer will begin counting rooted maps by genus, number of edges and number of vertices.  
 If you entered a large enough number of edges to make the program take a measurable amount of time to run, 
 you will see a progress report: starting genus 0, starting genus 1, etc., 
 to reassure you that the program is actually running.  When it has finished the calculations, 
 it will print on your screen the number of rooted maps of genus g with e edges and v vertices, 
 where g runs from 0 to the maximum you chose, e runs from 2g to the maximum you chose 
 and v runs from 1 to 1+e-2g for a map with one face.  
 If you chose to have the numbers entered into a file, the file can be found in the same folder as your program - 
 click on the folder 'build', then on any other folders that appear and one of them will contain that file.
 
 Then you will be prompted to enter 1 to count unrooted maps or 0 not to.  If you enter any integer except 0, 
 you will be prompted to enter 1 to store these numbers in a file called 'unrx.txt' or 0 not to, 
 and then the same sequence of events will take place as with rooted maps except for the name of the file.  
 The folder will also contain a file called 'RHx.txt', which the computer uses to count unrooted maps, 
 but you don't have to be concerned with it.
 
 When the computer has finished counting unrooted maps, or if you chose not to do so, 
 you will again be prompted to enter the maximum number of edges or -1 to quit; 
 so you can repeat the calculations with different choices until you decide to quit.
 
 If you'd like to receive tables of numbers or otherwise communicate with me, e-mail me at <walsh.timothy@uqam.ca>.
 
*/

#include "stdio.h"
#include "string.h"
#include <fstream>
#include <iostream>
#include <iomanip>
#include <cln/cln.h>
using namespace cln;
using std::cout;
using std::setw;
using std::endl;
#define edgemax 100 // maximum number of edges
#define emaxmax edgemax+1 // maximum number of edges + 1
#define emaxp1 edgemax+2 // maximum number of vertices + 1
#define genusmax (edgemax/2) // maximum genus
#define gmaxmax genusmax+2 // maximum genus + 2
#define term1 genusmax*(17+genusmax*(4*genusmax-15))/3
#define term2 genusmax*edgemax*(5-2*genusmax+edgemax)
#define term3 (term1+term2)/2
#define term4 ((edgemax-2*genusmax)*(edgemax-2*genusmax+1))/2
#define dimen2 term3+term4+edgemax-2*genusmax+1 // 1 + dimension of Rt2
#define rmax 2*gmaxmax+2 // dimension of array of parts of partition
#define lmax 4*gmaxmax+4 // dimension of array of phi
#define gqmax genusmax/2+1 // maximum genus of quotient map + 1
#define maxdiv 20 // maximum number of divisors of any integer <= 4*genusmax+2

void convertorb(int Lgr[3], int orb[rmax], int vf[maxdiv], int ol[maxdiv]);
cl_I readorb(int gl, int g, int Lgr[3], int orb[rmax]);
int FVcalc(int fv[maxdiv], int fvsum[1], int ol[maxdiv], int L);
int Ecalc(int e[2], int L);
void emult(int mne, int e[2], cl_I fact[1]);
void fvmult(int fv[maxdiv], int fvsum[1], cl_I fact[1], int vf[maxdiv]);
void firstvf(int vf[maxdiv], int v[maxdiv], int f[maxdiv], int vsum[1], int fsum[1]);
int nextvf(int vf[maxdiv], int v[maxdiv], int f[maxdiv], int vsum[1], int fsum[1]);
void process(int emax, int L, int g, int gl, cl_I A, int ol[maxdiv], int vf[maxdiv], cl_I Rt2[dimen2], cl_I Ut2[dimen2]);
int ind2(int g, int emax, int e, int f);
int root2(int gmax, int emax, cl_I Rt2[dimen2]);
void writeorbs(int L, int g, cl_I autos, int r, int d[rmax]);
int rh(int gl, int ipr, int gmax, int phi[lmax], int pr[rmax], int d[rmax],
	   int cg[gmaxmax], int s[gmaxmax], int gcd[gmaxmax], int q[lmax], int expos[lmax][gmaxmax]);
int partitions(int r, int t, int g, int ipr, int L, int pr[rmax], int d[rmax], int gcd[gmaxmax],
			   int q[lmax], int s[gmaxmax], int expos[lmax][gmaxmax]);
cl_I mhat2g(int g, int L, int ipr, int pr[rmax], int gcd[gmaxmax], int expos[lmax][gmaxmax]);
cl_I Jordan(int r, int g, int L, int ipr, int pr[rmax], int gcd[gmaxmax], int expos[lmax][gmaxmax]);
cl_I Ecal(int r, int L, int ipr, int pr[rmax], int s[gmaxmax],
		  int d[rmax], int gcd[gmaxmax], int expos[lmax][gmaxmax]);
int Harvey(int r, int d[rmax], int g, int ipr, int L, int pr[rmax],
		   int expos[lmax][gmaxmax], int gcd[gmaxmax], int s[gmaxmax]);
int makediv(int L, int ipr, int expos[lmax][gmaxmax], int q[lmax]);
void factorize(int gl, int ipr, int pr[rmax], int expos[lmax][gmaxmax]);
int phipr(int gl, int phi[lmax], long int phi2[lmax], int pr[rmax]);
int sz(cl_I x);
int isz(int x);
void pad(int d, cl_I x);

int bug;
std::fstream fp0, fp1, fp2;
char ascg,fname0[]="RHx.txt", fname1[]="rootx.txt", fname2[]="unrx.txt";

int main (void)
{
	int gmax, emax, e, gl, f, g, emin, i, szmax, ipr, L, count, flq;
	cl_I sum, m2g, J, E, autos; 
	int phi[lmax], pr[rmax], expos[lmax][gmaxmax], q[lmax], d[rmax], gcd[gmaxmax], s[gmaxmax];
	int  cg[gmaxmax];
	int unroot, dg, dd, indmax, gmx, indg, indgl, vf[maxdiv], ol[maxdiv];
	long int phi2[lmax];
	cl_I R2, U2, Rt2[dimen2], Ut2[dimen2], A, usum[emaxmax];
	int Lgr[3], orb[rmax];
	printf("This program counts maps by genus, number of edges and number of vertices.\n");
	printf("You choose the maximum number of edges and the maximum genus.\n");
	printf("The program will count rooted maps within these limits.\n");
	printf("You can then choose to put these numbers in a file called 'rootx.txt'\n");
	printf("The program will print the numbers and, if you have so chosen, put them in the file.\n");
	printf("Then you can choose to count unrooted maps by number of vertices and edges.\n");
	printf("If you do, you can choose to put those numbers in a file called 'unrx.txt'\n");
	printf("The program will print the numbers and, if you have so chosen, put them in the file.\n");
	while (1)
	{
		bug=1;
		printf("Enter maximum number of edges <= %d, -1 to quit: ", edgemax);
		scanf("%d",&emax);
		if (emax < 0) break;
		if (emax > edgemax) emax = edgemax;
		printf("Enter maximum genus <= %d: ",emax/2);
		scanf("%d",&gmax);
		if (gmax < 0) gmax = 0;
		if (gmax >= emax/2) gmax = emax/2;
		printf("Counting rooted maps.\n");
		printf("Enter 1 to store these numbers in a file or 0 not to: ");
		scanf("%d",&flq);		
		if (flq)
		{
			fp1.open(fname1, std::fstream::out | std::fstream::trunc);
			if (!fp1.is_open())
			{
				printf("Error opening output file.\n");
				flq=0;
			}
		}
		printf("g = genus, e = number of edges, v = number of vertices.\n");
		szmax = root2(gmax, emax, Rt2);
		if (szmax < 14) szmax = 14;
		printf("   g   e   v ");
		for (i=1;i<=szmax-13;i++) printf(" ");
		printf("# rooted maps\n");
		printf("   0   0   1");
		for (i=1;i<=szmax;i++) printf(" ");
		printf("1\n");
		printf("   0   0 sum");
		for (i=1;i<=szmax;i++) printf(" ");
		printf("1\n\n");
		if (flq)
		{
			fp1<<"   g   e   v ";
			for (i=1;i<=szmax-11;i++) fp1<<" ";
			fp1<<"# rooted maps\n";
			fp1<<"   0   0   1";
			for (i=1;i<=szmax+2;i++) fp1<<" ";
			fp1<<"1\n";
			fp1<<"   0   0 sum";
			for (i=1;i<=szmax+2;i++) fp1<<" ";
			fp1<<"1\n\n";
		}
		for (gl=0;gl<=gmax;gl++)
		{
			emin = 2*gl;
			if (gl==0) emin = 1;
			for (e=emin;e<=emax;e++)
			{
				sum=0;
				for (f=1;f<=e+1-2*gl;f++)
				{
					printf("%4d%4d%4d ",gl,e,f);
					R2 = Rt2[ind2(gl,emax,e,f)];
					pad(szmax,R2);
					cout << R2;
					sum = sum + R2;
					printf("\n");						
					if (flq)
					{
						for (i=isz(gl);i<=3;i++) fp1<<" ";
						fp1<<gl;
						for (i=isz(e);i<=3;i++) fp1<<" ";
						fp1<<e;
						for (i=isz(f);i<=3;i++) fp1<<" ";
						fp1<<f;
						for (i=sz(R2);i<=szmax+2;i++) fp1<<" ";
						fp1<<R2;
						fp1<<"\n";
					}
				} /* end for f */
				printf("%4d%4d sum ",gl,e);
				pad(szmax,sum);
				cout << sum;
				printf("\n\n");
				if (flq)
				{
					for (i=isz(gl);i<=3;i++) fp1<<" ";
					fp1<<gl;
					for (i=isz(e);i<=3;i++) fp1<<" ";
					fp1<<e;
					fp1<<" sum  ";
					for (i=sz(sum);i<=szmax;i++) fp1<<" ";
					fp1<<sum;
					fp1<<"\n\n";
				}
			} // end for e
		} // end for gl
		if (flq) fp1.close();
		// counting unrooted maps
		printf("Enter 1 to count unrooted maps, 0 not to: ");
		scanf("%d",&unroot);
		if (!unroot) continue;
		printf("Enter 1 to store these numbers in a file or 0 not to: ");
		scanf("%d",&flq);		
		if (flq)
		{
			fp2.open(fname2, std::fstream::out | std::fstream::trunc);
			if (!fp2.is_open())
			{
				printf("Error opening output file.\n");
				flq=0;
			}
		}
		indmax=ind2(gmax,emax,emax,1+emax-2*gmax);
		for (i=0;i<=indmax;i++) Ut2[i]=Rt2[i];
		ipr=phipr(genusmax,phi,phi2,pr);
		factorize(genusmax,ipr,pr,expos);
		dd=0;
		for (gl=0;gl<=gmax;gl++)
		{
			if (gl>0)
			{
				ipr=phipr(gl,phi,phi2,pr);
				count = rh(gl, ipr, gmax, phi, pr, d, cg, s, gcd, q, expos);
			} // end if gl>0	
			printf("Starting genus %d.\n",gl);
			if (gl==0)
			{
				L=2;
				vf[0]=1;
				ol[0]=1;
				vf[1]=2;
				ol[1]=1;
				g=0;
				A=1;
				process(emax,L,g,gl,A,ol,vf,Rt2,Ut2);
				for (L=3; L<=emax; L++)
				{
					A=phi[L];
					process(emax,L,g,gl,A,ol,vf,Rt2,Ut2);
				} // end for L
			}	
			else
			{
				if (gl==1)
				{	
					for (L=2; L*2<=emax; L++)
						for (e=2;e*L<=emax;e++)
							for (f=1;f<=e-1;f++)
							{
								indg = ind2(1,emax,e,f);
								R2 = Rt2[indg];
								indgl = ind2(gl,emax,L*e,L*f);
								Ut2[indgl]=Ut2[indgl]+R2*phi2[L];
							} // end for f, for e and for L
				} // end if gl==1
				gmx=(gl+1)/2;
				for (g=gmx;g>=0;g--)
				{
					fp0.open(fname0, std::fstream::in);
					if (!fp0.is_open()) {printf("Error opening file RHx.txt %d\n",gl); return gl;}
					while (1)
					{
						A = readorb(gl,g,Lgr,orb);
						if (A<0) break;
						if (g!=Lgr[1]) continue;
						L = Lgr[0];
						convertorb(Lgr,orb,vf,ol);
						process(emax,L,g,gl,A,ol,vf,Rt2,Ut2);
					}
					fp0.close();
				} // end for g
			} // end else
			for (e=2*gl;e<=emax;e++)
			{
				usum[e]=0;
				for (f=1;f<=e+1-2*gl;f++)
				{
					indgl=ind2(gl,emax,e,f);
					U2=Ut2[indgl];
					if (e>0) Ut2[indgl]=exquo(U2,2*e); else Ut2[indgl]=1;
					usum[e]=usum[e]+Ut2[indgl];
				} // end for f
			} // end for e
			dg=sz(usum[emax]);
			if (dg>dd) dd=dg;
		} // end for gl		
		if (dd<14) dd=14;
		printf("   g   e   v ");
		for (i=1;i<=dd-13;i++) printf(" ");
		printf("# unrooted maps\n");
		if (flq)
		{
			fp2<<"   g   e   v ";
			for (i=1;i<=dd-13;i++) fp2<<" ";
			fp2<<"# unrooted maps\n";
		}
		for (gl=0;gl<=gmax;gl++)
		{
			emin = 2*gl;
			if (gl==0) emin = 1;
			for (e=emin;e<=emax;e++)
			{
				sum=0;
				for (f=1;f<=e+1-2*gl;f++)
				{
					printf("%4d%4d%4d   ",gl,e,f);
					U2 = Ut2[ind2(gl,emax,e,f)];
					pad(dd,U2);
					cout << U2;
					sum = sum + U2;
					printf("\n");						
					if (flq)
					{
						for (i=isz(gl);i<=3;i++) fp2<<" ";
						fp2<<gl;
						for (i=isz(e);i<=3;i++) fp2<<" ";
						fp2<<e;
						for (i=isz(f);i<=3;i++) fp2<<" ";
						fp2<<f;
						for (i=sz(U2);i<=dd+2;i++) fp2<<" ";
						fp2<<U2;
						fp2<<"\n";
					}
				} /* end for f */
				printf("%4d%4d sum ",gl,e);
				pad(dd+2,sum);
				cout << sum;
				printf("\n\n");
				if (flq)
				{
					for (i=isz(gl);i<=3;i++) fp2<<" ";
					fp2<<gl;
					for (i=isz(e);i<=3;i++) fp2<<" ";
					fp2<<e;
					fp2<<" sum  ";
					for (i=sz(sum);i<=dd;i++) fp2<<" ";
					fp2<<sum;
					fp2<<"\n\n";
				}
			} // end for e
		} // end for gl
		if (flq) fp2.close();
	} // end while
	printf("End of execution.");
    return 0;
}

void convertorb(int Lgr[3], int orb[rmax], int vf[maxdiv], int ol[maxdiv])
// expresses an orbifold signature as orbit lengths and multiplicities of branch points
{
	int j, k, old;
	
	old=0;
	k=0;
	for (j=0; j<maxdiv; j++) {vf[j]=0; ol[j]=0;}
	for (j=0; j<Lgr[2]; j++)
	{
		if (orb[j]!=old)
		{
			k++;
			old=orb[j];
			ol[k]=old;
		}
		vf[k]++;
	}
	ol[0]=k;
	vf[0]=k;
}

cl_I readorb(int gl, int g, int Lgr[3], int orb[rmax])
// reads a single orbifold signature from a file
{
	int j; cl_I A;
	
	fp0>>Lgr[0]; //period
	if (Lgr[0]<0) return -1;
	fp0>>Lgr[1]; //genus
	fp0>>A; //# automorphisms
	fp0>>Lgr[2]; //# branch points
	for (j=0;j<Lgr[2];j++) fp0>>orb[j];
	return A;
}
	

int FVcalc(int fv[maxdiv], int fvsum[1], int ol[maxdiv], int L)
// Calculates the number of vertices or faces in the lifted map.
{
	int i,s;
	
	s=0;
	for (i=1;i<=ol[0];i++) s=s+fv[i]*ol[i];
	return L*(fv[0]-fvsum[0])+s;
}

int Ecalc(int e[2], int L)
// Calculates the number of edges in the lifted map.
{
	return L*e[0]+L*e[1]/2;
}

void emult(int mne, int e[2], cl_I fact[1])
// Calculates the number of insertions of semiedges in the quotient map * the rooting ratio.
{
	int i, num, den;
	
	if (e[0]==mne)
	{
		num=2*e[0]+1; den=1; fact[0]=1;
		for (i=1;i<=e[1];i++)
		{
			fact[0]=exquo(fact[0]*num,den);
			num++;
			den++;
		}	
	}
	else fact[0]=exquo(fact[0]*(2*e[0]+e[1]-1)*(2*e[0]+e[1]),(2*e[0]-1)*(2*e[0]));
}

void fvmult(int fv[maxdiv], int fvsum[1], cl_I fact[1], int vf[maxdiv])
// Calculates the number of ways to distribute the branch points among the vertices or faces.
{
	int i, j, num, den;
	num=fv[0]; den=1; i=1;
	
	if (fv[0]<fvsum[0])
	{
		fact[0]=0;
		return;
	}
	else if (fv[0]==fvsum[0])
	{
		fact[0]=1;
		for (j=1; j<=fvsum[0]; j++)
		{
			fact[0]=exquo(fact[0]*num,den);
			num--;
			while (fv[i]==0) {i++; if (i>vf[0]) break;}
			if (i>vf[0]) break;
			if (den<fv[i]) den++;
			else
			{
				i++;
				den=1;
			}
		}
	}
	else fact[0]=exquo(fact[0]*fv[0],(fv[0]-fvsum[0]));
}

void firstvf(int vf[maxdiv], int v[maxdiv], int f[maxdiv], int vsum[1], int fsum[1])
// Calculates the initial numbers of branch points in the vertices and faces.
{
	int i;
	
	vsum[0] = 0; fsum[0]=0;
	for (i=1; i<=vf[0]; i++)
	{
		v[i]=0; f[i]=vf[i]; fsum[0]=fsum[0]+f[i];
	}
}

int nextvf(int vf[maxdiv], int v[maxdiv], int f[maxdiv], int vsum[1], int fsum[1])
// Updates the numbers of branch points in the vertices and faces.
{
	int i;
	
	for (i=vf[0]; i>0; i--)
	{
		if (f[i]>0)
		{
			f[i]--; fsum[0]--;
			v[i]++; vsum[0]++;
			break;
		}
		else
		{
			fsum[0]=fsum[0]+vf[i]-f[i];
			f[i]=vf[i];
			vsum[0]=vsum[0]-v[i];
			v[i]=0;
		}
	}
	return i;
}

void process(int emax, int L, int g, int gl, cl_I A, int ol[maxdiv], int vf[maxdiv], cl_I Rt2[dimen2], cl_I Ut2[dimen2])
// Calculates the contribution of all the rooted genus-g quotient maps and the A orbifolds of period L
// with signatures expressed by orbit lengths in ol and multiplicities in vf (vf[0]=ol[0]=# orbit lengths)
// to the number of unrooted maps with at most emax edges.
{
    int vf1,indu, indr, vfo,e[2],v[maxdiv],f[maxdiv],vsum[1],fsum[1],vgo,fgo,vmin,fmin,V,F,E,mne;
    cl_I vmul[1],fmul[1],emul[1];
	
	vfo=vf[1];
	if (vf[0]>0&&2*ol[1]==L) vf1=vf[1]; else vf1=0; // If there is a branch index 2, there may be semiedges.
	for (e[1]=0;e[1]<=vf1;e[1]++)  //e[0]=# edges of quotient map, e[1]=# semiedges
	{			
		firstvf(vf,v,f,vsum,fsum); //v and f contain the initial number of vertices and faces in the branch points
		vmul[0]=1;
		if (vsum[0]==0) vmin=1; else vmin=vsum[0];
		vgo=1; v[0]=vmin;
		while (vgo) //Number of vertices of quotient map increases until the lifted map has too many edges.
		{
			fvmult(v,vsum,vmul,vf);
		    V=FVcalc(v,vsum,ol,L);
			fmul[0]=1;
			if (fsum[0]==0) fmin=1; else fmin=fsum[0];
			mne=fmin+v[0]-2+2*g;
			fgo=1; f[0]=fmin;
			while (fgo) //Number of faces of quotient map increases until the lifted map has too many edges.
			{
				fvmult(f,fsum,fmul,vf);
				F=FVcalc(f,fsum,ol,L);
				e[0]=f[0]+v[0]-2+2*g;
				emult(mne,e,emul);
				E=Ecalc(e,L);
				if(E<=emax&&V>0&&F>0)
				{	
					indu = ind2(gl,emax,E,V);
					indr = ind2(g,emax,e[0],v[0]);
					Ut2[indu]=Ut2[indu]+Rt2[indr]*vmul[0]*fmul[0]*emul[0]*A;
					f[0]++;
				}
				else fgo=0;
			} // end while (fgo)
			if (f[0]>fmin) v[0]++; else vgo=0;
		} // end while (vgo)	
		while (1) // Updates the distribution of branch points among the faces and vertices.
		{
            if (nextvf(vf,v,f,vsum,fsum)<=0) break;
			vmul[0]=1;
			if (vsum[0]==0) vmin=1; else vmin=vsum[0];
			vgo=1; v[0]=vmin;
			while (vgo)
			{
				fvmult(v,vsum,vmul,vf);
				V=FVcalc(v,vsum,ol,L);
			    fmul[0]=1;
			    if (fsum[0]==0) fmin=1; else fmin=fsum[0];
				mne=fmin+v[0]-2+2*g;
				fgo=1; f[0]=fmin;
				while (fgo)
			    {
					fvmult(f,fsum,fmul,vf);
					F=FVcalc(f,fsum,ol,L);
					e[0]=f[0]+v[0]-2+2*g;
					emult(mne,e,emul);
					E=Ecalc(e,L);
					if(E<=emax&&V>0&&F>0)
					{	
						indu = ind2(gl,emax,E,V);
						indr = ind2(g,emax,e[0],v[0]);
						Ut2[indu]=Ut2[indu]+Rt2[indr]*vmul[0]*fmul[0]*emul[0]*A;
						f[0]++;
					}
					else fgo=0;
				} // end while (fgo)
				if (f[0]>fmin) v[0]++; else vgo=0;
			} // end while (vgo)						
		}  // end while (1) - updating the distribution of branch points among vertices and faces
		vf[1]--;
	}// end for (e[1]) = # semiedges.
	vf[1]=vfo;
} // end process

int ind2(int g, int emax, int e, int f)
{
	int g2, t0, t1, t2, t3;
	
	g2=2*g;
	t0 = (g*(17+g*(4*g-15))/3 + g*emax*(5-g2+emax))/2;
	t1 = ((e-g2)*(e-g2+1))/2;
	t2 = f-1;
	t3 = t0+t1+t2;
	return t3;
} // end ind2

int root2(int gmax, int emax, cl_I Rt2[dimen2])
{
	int i, indmax, g, e, f, f1min, f1max, v, g1, g2, e1, e2, f1, f2, emin, szc, szmax;
	cl_I s, sum;
	indmax=ind2(gmax,emax,emax,1+emax-2*gmax);
	for (i=1;i<=indmax;i++)
		Rt2[i]=0;
    Rt2[0]=1;
    szmax = 0;
	for (g=0;g<=gmax;g++)
    {
		printf("Starting genus %d.\n",g);
		emin = 2*g;
    	if (g==0) emin = 1;
		for (e=emin;e<=emax;e++)
		{
			sum = 0;
			for (f=1;f<=e+1-2*g;f++)
			{
				s=0;
				if (e>2*g)
				{	
					v = e + 2*(1-g) - f;
					if (v>1) s = s + (4*e-2)*Rt2[ind2(g,emax,e-1,f)];
					if (f>1) s = s + (4*e-2)*Rt2[ind2(g,emax,e-1,f-1)];
				}
				if (g>0) s = s + (e-1)*(2*e-1)*(2*e-3)*Rt2[ind2(g-1,emax,e-2,f)];
				for (g1=0;g1<=g;g1++)
				{
					g2 = g - g1;
					for (e1=2*g1; e1<=e-2-2*g2; e1++)
					{
						e2 = e - 2 - e1;
						f1min = f - e2 - 1 + 2*g2;
						if (f1min<1) f1min=1;
						f1max = e1 + 1 - 2*g1;
						if (f1max > f-1) f1max=f-1;
						for (f1=f1min; f1<=f1max;f1++)
						{
							f2 = f - f1;
							s = s + 3*(2*e1+1)*(2*e2+1)*Rt2[ind2(g1,emax,e1,f1)]*Rt2[ind2(g2,emax,e2,f2)];
						} // end for f1
					} // end for e1
				} // end for g1
				s = exquo(s,(e+1));
				Rt2[ind2(g,emax,e,f)] = s;
				sum = sum + s;
			}// end for f
			szc = sz(sum);
			if (szc>szmax) szmax = szc;
		} // end for e
	} // end for g
	return szmax;
} // end root2

void writeorbs(int L, int g, cl_I autos, int r, int d[rmax])
// Stores the orbifold signatures in a file
{
	int n;
	fp0<<L;
	fp0<<" ";
	fp0<<g;
	fp0<<" ";
	fp0<<autos;
	fp0<<" ";
	fp0<<r;
	fp0<<" ";
	for (n=1;n<=r;n++)
	{
		fp0<<d[n];
		fp0<<" ";
	}
	fp0<<"\n";
}

int rh(int gl, int ipr, int gmax, int phi[lmax], int pr[rmax], int d[rmax],
	   int cg[gmaxmax], int s[gmaxmax], int gcd[gmaxmax], int q[lmax], int expos[lmax][gmaxmax])
/* generates all solutions to the Riemann-Hurwitz equation for given gamma=gl */
{
	int L,r,i,t,g,ndiv,count,lmx,gmx,rnum,rmn,rmx,npts;
	
	fp0.open(fname0, std::fstream::out | std::fstream::trunc);
	if (!fp0.is_open()) printf("Error opening output file.\n");
	count=0;
	for (i=0;i<=gmax;i++) cg[i]=0;
	lmx=4*gl+2;
	for (L=2;L<=lmx;L++)
	{
		if (phi[L]>2*gl) continue;
		if(gl!=1&&(gl-1)%L==0)
		{
			g=1+(gl-1)/L; r=0; t=0;
			npts=partitions(r, t, g, ipr, L, pr, d, gcd, q, s, expos);
			cg[g]=cg[g]+npts;
			count=count+npts;
		} /* end if gl */
		ndiv=makediv(L, ipr, expos, q);
		gmx=(gl-1+q[1])/L;
		for (g=gmx;g>=0;g--)
		{
			rnum=2*(gl+L-g*L-1);
			if (rnum%(L-1)==0) rmn=rnum/(L-1); else rmn=1+rnum/(L-1);
			if (rmn<2) rmn=2;
			rmx=rnum/(L-q[1]);
			if (rmx>2*gl+2) rmx=2*gl+2;
			for (r=rmn;r<=rmx;r++)
			{
				t=2-2*gl-L*(2-2*g-r);
				if (t<0) continue;
				npts=partitions(r, t, g, ipr, L, pr, d, gcd, q, s, expos);
				cg[g]=cg[g]+npts;
				count=count+npts;
			} /* end for r */
		} /* end for g */
	} /* end for L */
	fp0<<-1;
	fp0.close();
	return count;
} /* end rh */

int partitions(int r, int t, int g, int ipr, int L, int pr[rmax], int d[rmax], int gcd[gmaxmax],
			   int q[lmax], int s[gmaxmax], int expos[lmax][gmaxmax])
/* Generates in inverse lexicographical order all the partitions d[1..r] of t with r parts
 from the set q, expressed as a strictly decreasing sequence of integers ending with 1. */
{
	int i, n, k, count, j[rmax];
	cl_I E, J, m2g, autos;
	
	if (t==0&&r==0)
	{
		autos=Jordan(r, g, L, ipr, pr, gcd, expos);
		writeorbs(L, g, autos, r, d);
		return 1; /* empty sequence */
	} /* end if t */
	if (r>t||q[1]*r<t) return 0; /* no sequences */
	if (r==t) /* unique sequence of 1 */
	{
		for (n=1;n<=r;n++) d[n]=1;
		k = Harvey(r, d, g, ipr, L, pr, expos, gcd, s);
		if (k>0)
		{
			E=Ecal(r, L, ipr, pr, s, d, gcd, expos);
			J=Jordan(r, g, L, ipr, pr, gcd, expos);
			m2g=mhat2g(g, L, ipr, pr, gcd, expos);
			autos = E*J*m2g;
			writeorbs(L, g, autos, r, d);
			return 1;
		} /* end if k */
	} /* end if r */
	if (q[1]*r==t) /* unique sequence of largest part */
	{
		for (n=1;n<=r;n++) d[n]=q[1];
		k = Harvey(r, d, g, ipr, L, pr, expos, gcd, s);
		if (k>0)
		{
			E=Ecal(r, L, ipr, pr, s, d, gcd, expos);
			J=Jordan(r, g, L, ipr, pr, gcd, expos);
			m2g=mhat2g(g, L, ipr, pr, gcd, expos);
			autos = E*J*m2g;
			writeorbs(L, g, autos, r, d);
			return 1;
		} /*  end if k */
	} /* end if q */
	count=0; i=1; j[1]=1;
	while (1)
	{
		d[i]=q[j[i]]; t=t-d[i]; /* assigning part #i */
		if (d[i]==1||t>d[i]*(r-i)) /* d[i] is too small; backtrack */
		{
			t=t+d[i]; i--;
			if (i<=0) break; /* All the partitions have been generated. */
			t=t+d[i]; j[i]++; /* decrease d[i] */
		} /* end backtrack */
		else /* d[i] isn't too small. */
		{
			if (t>r-i) {i++; j[i]=j[i-1];} /* advance */
			else /* t<=r-i */
			{
				if (t==r-i) /* We have a partition. */
				{
					for (k=i+1;k<=r;k++) d[k]=1;
					k = Harvey(r, d, g, ipr, L, pr, expos, gcd, s);
					if (k>0)
					{
						count++;
						E=Ecal(r, L, ipr, pr, s, d, gcd, expos);
						J=Jordan(r, g, L, ipr, pr, gcd, expos);
						m2g=mhat2g(g, L, ipr, pr, gcd, expos);
						autos = E*J*m2g;
						writeorbs(L, g, autos, r, d);
					} /* end if k */
				} /* end if t */
				t=t+d[i]; j[i]++; /*decrease d[i]*/
			} /* end else (t<=r-i) */
		} /* end else (d[i] isn't too small.) */
	}  /* end while */
	return count;
} /* end partitions */

cl_I mhat2g(int g, int L, int ipr, int pr[rmax], int gcd[gmaxmax], int expos[lmax][gmaxmax])
/* calculates lcm(m[1],...,m[r])^2g */
{
	int i,k,ar;
	cl_I m,m2g;
	
	if (g==0) return 1;
	m=1;
	for (i=0;i<=ipr;i++)
	{
		ar=expos[L][i]-gcd[i];
		if (ar>0)
			for (k=1;k<=ar;k++) m=m*pr[i]; /* m=lcm(m[1],...,m[r]) */
	} /* end for i */
	m=m*m; m2g=m;
	for (k=1;k<g;k++) m2g=m2g*m;
	return m2g;
} /* end mhat2g */

cl_I Jordan(int r, int g, int L, int ipr, int pr[rmax], int gcd[gmaxmax], int expos[lmax][gmaxmax])
/* calculates phi(2g,gcd(d[1],...,d[r]) */
{
	int i,p2,e;
	cl_I p2g,J;
	
	if (g==0) return 1;
	J=1;
	for (i=0;i<=ipr;i++) /* pr[i] is the ith prime */
	{
		if (r==0) gcd[i]=expos[L][i];
		if (gcd[i]>0) /* the ith prime divides gcd */
		{
			p2=pr[i]*pr[i]; p2g=p2;
			for (e=1;e<g;e++) p2g=p2g*p2; /* p2g=pr[i]^2g */
			for (e=1;e<gcd[i];e++) J=J*p2g;
			J=J*(p2g-1);
		} /* end if */
	} /* end for */
	return J;
} /* end Jordan */

cl_I Ecal(int r, int L, int ipr, int pr[rmax], int s[gmaxmax],
		  int d[rmax], int gcd[gmaxmax], int expos[lmax][gmaxmax])
/* calculates E(m[1],...,m[r]) */
{
	int i,j,k,rp,vp,rpd,vpd;
	cl_I pw,E;
	
	if (r==0) return 1;
	E=1;
	for (i=0;i<=ipr;i++)
	{
		rpd=expos[L][i]-gcd[i];
		if (rpd>0) /* ith prime divides lcm(m[1],...,m[r]) */
		{
			pw=1;
			for (k=1;k<s[i];k++) pw=pw*(pr[i]-1);
			if (s[i]%2==0) pw=pw+1; else pw=pw-1;
			pw=exquo(pw,pr[i]);
			E=E*pw;
			rp=0;
			vp=1;
			for (j=1;j<=r;j++)
			{
				vpd=expos[L][i]-expos[d[j]][i];
				if (vpd>0) { rp++;vp=vp+vpd-1; } /* ith prime divides m[j] */
			} /* end for j */
			vp=vp-rpd;
			for (k=1;k<=vp;k++) E=E*pr[i];
			for (k=1;k<=rp-s[i]+1;k++) E=E*(pr[i]-1);
		} /* end if */
	}
	return E;
} /* end Ecalc */

int Harvey(int r, int d[rmax], int g, int ipr, int L, int pr[rmax],
		   int expos[lmax][gmaxmax], int gcd[gmaxmax], int s[gmaxmax])
/* tests whether (d[1],...,d[r]) satisfies Harvey's condition */
{
	int i,j,k;
	
	if (g==0)
		for (i=0;i<=ipr;i++)
			if (expos[L][i]>0) /* the ith prime pr[i] divides L. */
			{
				k=0;
				for (j=1;j<=r;j++)
					if (expos[d[j]][i]==0) /* the ith prime doesn't divide d[j] */
					{ k=1; break; }
				if (k==0) /* the ith prime divides all the d[j]s */
				{
					return 0;
				} /* end if expos and for j and if expos and for i and if g */
			} /* end if expos and for i and if g */
	gcd[0]=expos[L][0];
	if (expos[L][0]>0) /* L is even */
	{
		for (j=1;j<=r;j++) /* computing the exponent of 2 in gcd(d[1],...,d[r]) */
			if (expos[d[j]][0]<gcd[0]) gcd[0]=expos[d[j]][0];
		k=0;
		for (j=1;j<=r;j++) /* counting how many d[j] have the same exponent of 2 */
			if (gcd[0]==expos[d[j]][0]) k++;
		s[0]=k;
		if (k%2==1&&expos[L][0]>gcd[0]) return 0; /* s[0] is odd and L/gcd([1],...,d[r]) is even */
	} /* end outer if expos */
	for (i=0;i<=ipr;i++) /* calculating the prime power decomposition of gcd(d[1],...,d[r]) */
	{
		gcd[i]=expos[L][i];
		s[i]=0;
		if (expos[L][i]>0) /* the ith prime divides L */
		{
			for (j=1;j<=r;j++) /* computing the exponent of the ith prime in gcd(d[1],...,d[r]) */
				if (expos[d[j]][i]<gcd[i]) gcd[i]=expos[d[j]][i];
			k=0;
			for (j=1;j<=r;j++) /* counting how many d[j] have the same exponent of the ith prime */
				if (gcd[i]==expos[d[j]][i]) k++;
			s[i]=k;
			if (k==1)
			{
				return 0;
			} /* end if k */
		} /* end expos */
	} /* end for i */
	return 1;
} /* end Harvey */

int makediv(int L, int ipr, int expos[lmax][gmaxmax], int q[lmax])
/* q contains the positive divisors of L in decreasing order */
{
	int ndiv,i,j,f;
	
	ndiv=1;
	for (i=0;i<=ipr;i++) ndiv=ndiv*(expos[L][i]+1); /* ndiv = #divisors of L */
	i=0; j=ndiv-1;
	f=1;
	while (i<=j)
	{
		if (L%f==0)
		{
			q[j]=f; j--;
			q[i]=L/f; i++;
		} /* end if */
		f++;
	} /* end while */
	return ndiv;
} /* end makediv */

void factorize(int gl, int ipr, int pr[rmax], int expos[lmax][gmaxmax])
/* creates a table expos of exponents of prime factors <= 2gl+2 of +ve integers <=4gl+2 */
{
	int i,j,ip,p,lmx,prpwr[lmax];
	
	for (i=0;i<lmax;i++)
	{
		prpwr[i]=-1;
		for (j=0;j<gmaxmax;j++) expos[i][j]=0;
	} /* end for i */
	ip=0; p=2; lmx=4*gl+2;
	for (i=2;i<=lmx;i++)
	{
		if (i==p) /* i is a prime */
		{
			prpwr[i]=ip;
			j=i;
			while (j<=lmx)
			{
				prpwr[j]=ip;
				j=j*i;
			} /* end while j */
			ip++;
			if (ip>ipr) p=0; else p=pr[ip];
		} /* end if i */
	} /* end for i */
	for (i=2;i<=lmx;i++)
	{
		if (prpwr[i]>=0)
		{
			j=i; ip=prpwr[i];
			while (j<=lmx)
			{
				expos[j][ip]++;
				j=j+i;
			} /* end while j */
		} /* end while */
	} /* end for i */
} /* end factorize */

int phipr(int gl, int phi[lmax], long int phi2[lmax], int pr[rmax])
/* creates a table phi of the Euler function and a list pr of primes */
{
	int imax,pmax,i,j,ipr;
	
	imax=4*gl+5; pmax=2*gl+3;
	for (i=0;i<=imax;i++) {phi[i]=i; phi2[i]=i*i;}
	ipr=-1;
	for (i=2;i<=imax;i++)
	{
		if (phi[i]==i) /* i is a prime */
		{
			if (i<=pmax)
			{
				ipr++;
				pr[ipr]=i;
			} /* end if i<=pmax */
			j=i;
			while (j<=imax)
			{
				phi[j]=phi[j]-phi[j]/i;
				phi2[j]=phi2[j]-phi2[j]/(i*i);
				j=j+i;
			} /* end while j */
		} /* end if */
	} /* end for */
	return ipr;
} /* end phipr */

int sz(cl_I x)
/* Computes the number of characters occupied by an integer x of type cl_I */
{
	int s; cl_I y;
	
	y=x;
	if (y==0) return 1;
	if (y<0)
	{
		s=2;
		y=-y;
	}
	else s=1;
	while (y>=10)
	{
		y=exquo(y-mod(y,10),10);
		s++;
	}
	return s;
}

int isz(int x)
/* Computes the number of characters occupied by an integer x of type int */
{
	int s, y;
	
	y=x;
	if (y==0) return 1;
	if (y<0)
	{
		s=2;
		y=-y;
	}
	else s=1;
	while (y>=10)
	{
		y=y/10;
		s++;
	}
	return s;
}

void pad(int d, cl_I x)
/* prints blanks to the left of a big integer x so that it occupies at least d spaces */
{
	int i;
	for (i=sz(x);i<d;i++) printf(" ");
}