#include <stdlib.h>
#include <iostream>

using namespace std;

#include <math.h>

typedef int *ptr_int;
typedef ptr_int *ptr_ptr_int;

int n=4;

class list{
public:
  list *first;
  list *next;
  int k,number;
  int lcm;
  list()
  {
    next=NULL;
  }

  void new_list(int new_lcm,int k1,int n1);
  void check(int new_lcm,int k1,int n1);

  void show();

};

void list::show()
{
  cout << k << " " << number << " with lcm " << lcm << endl;
  if( next!=NULL )
    next->show();
}


void list::new_list(int new_lcm,int k1,int n1)
{
  next=new list;
  next->first=first;
  next->k=k1;
  next->number=n1;
  next->lcm=new_lcm;
}

void list::check(int new_lcm,int k1,int n1)
{
  
  if( new_lcm==lcm )
    {
      k=k1; number=n1;
    }
  else
    {
      if( next!=NULL )
	{
	  next->check(new_lcm,k1,n1);
	}
      else
	{
	  // create a new member in the list
	  new_list(new_lcm,k1,n1);
	}      
    }
}




int fact(int k)
{
  if( k==0 )
    return 1;
  else
    return k*fact(k-1); 
}

int binom(int n, int k)
{
  return fact(n)/fact(k)/fact(n-k);
}

int gcd(int k, int l)
{
  int aux;
  
  if( l>k )
    {
      aux=l; l=k; k=aux;
      return gcd(k,l);
    }
  if( k%l==0 )
    return l;
  else 
    return gcd(l,k%l);
  
}

int gcdset(int *exponent,ptr_ptr_int *b, int k, int count)
{
  int i;
  int current_gcd=1;
  for( i=0; i<=k; i++)
    {
      current_gcd*=exponent[b[i][k][count]];
    }

  for( i=0; i<=k; i++)
    {
      current_gcd=gcd(exponent[b[i][k][count]],current_gcd);
    }
  return current_gcd;
}





int gcdsetcomplement(int *exponent,ptr_ptr_int *b, int k, int count)
{
  int i,j,found;
  int current_gcd=1;

  for( i=0; i<n; i++)
    {
      current_gcd*=exponent[b[i][n-1][0]];
    }

  for( i=0; i<n; i++)
    {
      found=0; j=0;
      while( !found && j<=k )
	{
	  if( b[j][k][count]==b[i][n-1][0] )
	    found=1;
	  j++;
	}
      if( !found )
	{
	  current_gcd=gcd(exponent[b[i][n-1][0]],current_gcd);
	}
    }

  return current_gcd;
}



int lcm(int k, int l)
{
  return k*l/gcd(k,l);
}

int lcmset(int *exponent,ptr_ptr_int *b, int k, int count)
{
  int i;
  int current_lcm=1;
  for( i=0; i<=k; i++)
    {
      
      current_lcm=lcm(exponent[b[i][k][count]],current_lcm);
    }
  return current_lcm;
}

int prodset(int *exponent,ptr_ptr_int *b, int k, int count)
{
  int i;
  int prod=1;
  for( i=0; i<=k; i++)
    {
      prod*=exponent[b[i][k][count]];
    }
  return prod;
}

int setcheck(ptr_ptr_int *b,int k,int count,int i,int l)
{
  int j1,j2,found;


  for( j2=0; j2<=i; j2++ )
    {
      j1=0; found=0;
      while( !found && j1<=k )
	{
	  if( b[j2][i][l]==b[j1][k][count] )
	    {
	      found=1; 
	    }
	  j1++;
	}
      if( !found )
	return 0;
    }
  return 1;
}




main()
{
  cout << "Dimension (enter number of exponents): ";
  cin >> n;

  int i,count,k,l;
  int loop;
  int *a,*last;
  ptr_int *kappa,*lcms,*tkappa,*C;
  ptr_ptr_int *b;
  int *exponent;
  int sign,rmax,d;
  int found,maxlcm;
  list *orbits,*first;
  int *chrank;
  int mindeg,mindeg2,maxdeg,maxdeg2,period,degree;

  
 
  a=new int[n];
  last=new int[n];
  kappa=new ptr_int[n];
  tkappa=new ptr_int[n];
  lcms=new ptr_int[n];
  C=new ptr_int[n];

  b=new ptr_ptr_int[n];
  exponent=new int[n];
  
  orbits=new list;
  


  for( i=0; i<n; i++ )
    {
      cout << "exponent a_" << i << ": ";
      cin >> exponent[i]; 
    }
  
 
  for( l=0; l<n; l++ )
    {
      b[l]=new ptr_int[n];
      for( i=0; i<n; i++ )
	{
	  (b[l])[i]=new int[binom(n,i+1)];
	}
      kappa[l]=new int[binom(n,l+1)];
      tkappa[l]=new int[binom(n,l+1)];
      lcms[l]=new int[binom(n,l+1)];
      C[l]=new int[binom(n,l+1)];
    }
  
  
  
  for( i=0; i<n; i++ )
    {
      last[i]=binom(n,i+1)-1;
      for( k=0; k<=i; k++ )
	{
	  ((b[k])[i])[last[i]]=(n-1)-i+k;
	}
      for( k=0; k<=i; k++ )
	{
	  ((b[k])[i])[0]=k;
	}
    }
 
  // subsets
  for( i=0; i<n; i++ )
    {
      
      for( count=1; count<binom(n,i+1)-1; count++ )
	{
	  loop=1;
	  k=i;    
	  for( l=0; l<i; l++ )
	    {
	      ((b[l])[i])[count]=((b[l])[i])[count-1];
	    }
	  ((b[k])[i])[count]=((b[k])[i])[count-1]+1;
	  while( loop )
	    {
	      if( ((b[k])[i])[count]==((b[k])[i])[last[i]]+1 && k!=0 )
		{
		  ((b[k-1])[i])[count]=((b[k-1])[i])[count-1]+1;
		  k--;
		}
	      else
		loop=0;
	    }
	  
	  //reset
	  k++;
	  while( k<=i )
	    {
	      ((b[k])[i])[count]=((b[k-1])[i])[count]+1;	    
	      k++;
	    }
	    
	}
    }


  // compute the lcms for Randell's algorithm
  for( k=0; k<n; k++ )
    {
      for( count=0; count<binom(n,k+1); count++ )
	{
	  lcms[k][count]=lcmset(exponent,b,k,count);
	}
    }



  //  determine kappa
  rmax=0;
  for( k=0; k<n; k++ )
    {
      for( count=0; count<binom(n,k+1); count++ )
	{

	  sign=(k+1)%2;
	 if( sign==0 )
	   sign=+1;
	 else
	   sign=-1; 
	  kappa[k][count]=sign;
	 	  
	  // sum over subsets with i+1 elements
	  for(i=0; i<k; i++)
	    {
	      // sets with i+1 elements
	      for( l=0; l<binom(n,i+1); l++ )
		{
		  //check whether set is a subset of b[-][k][count]
		  if( setcheck(b,k,count,i,l) )
		    {
		      sign=(k-i)%2;
		      if( sign==0 )
			sign=+1;
		      else
			sign=-1;
		      kappa[k][count]+=sign*prodset(exponent,b,i,l)/lcmset(exponent,b,i,l);
		      
		    }
		}	      
	    }
	  // set b[-][k][count] always occurs as a subset of itself
	  kappa[k][count]+=prodset(exponent,b,k,count)/lcmset(exponent,b,k,count);

	  
	  sign=(n+1-k)%2;
	  if( sign==0 )
	    tkappa[k][count]=0;
	  else
	    tkappa[k][count]=kappa[k][count];

	  if( tkappa[k][count]>rmax )
	    rmax=tkappa[k][count];
  
	    
	}
    }

  //n odd implies rmax is at least 1 (empty set contributes in that case)
  if( n%2!=0 )
    {
      if( rmax==0 )
	rmax=1;
    }


  //cout << "C's" << endl;
  for( k=0; k<n; k++ )
    {
      for( count=0; count<binom(n,k+1); count++ )
	{
	  C[k][count]=gcdsetcomplement(exponent,b,k,count)/gcdset(exponent,b,n-1,0);
	  //cout << "   " << k << " " << count << ":    " << C[k][count] << "   ";
	 
	  // subsets with i+1 elements
	  for(i=0; i<k; i++)
	    {
	      // sets with i+1 elements
	      for( l=0; l<binom(n,i+1); l++ )
		{
		  //check whether set is a subset of b[-][k][count]
		  if( setcheck(b,k,count,i,l) )
		    {
		      C[k][count]/=C[i][l];
		      //	      cout << C[k][count] << "   ";
		    }
		}	      
	    }
	  //cout << endl;
	  
	    
	}
    }


  


  // compute Milnor number
  d=1;
  for( i=0; i<n; i++ )
    {
       d=d*(exponent[i]-1);
    }
  
  cout << "The milnor number of the natural filling of Sigma is " << d << endl;
  cout << "The rank of H_" << n-2 << "(Sigma;Z) is " << kappa[n-1][0] << endl;
  
  if( rmax>0 )
    {
  cout << "The torsion part of H_" << n-2 << "(Sigma;Z) consists of groups with order " << endl; 
  for( i=0; i<rmax; i++ )
    {
      d=1;
      for( k=0; k<n; k++ )
	{
	  
	    for( count=0; count<binom(n,k+1); count++ )
	    {
	      if( tkappa[k][count]>=i+1 )
		{
		  // cout << "bijdrage bij" << k << " " << count << endl;
		  d*=C[k][count];
		}
	    }
	}
      //empty set counts once if n is odd, tkappa=0	  
      if( n%2!=0 )
	d*=gcdset(exponent,b,n-1,0);
      if( d!=1 )
	cout << "d_" << i << "=" << d << endl;
    }
    cout << "There are no other torsion subgroups" << endl;
     }
   else
  {
    cout << "H_" << n-2 << "(Sigma;Z) is torsion-free" << endl;
  }

  cout << endl;

  //lcm determines orbit type
  orbits->lcm=lcms[1][0];
  orbits->k=1; orbits->number=0;
  for( k=1; k<n; k++ )
    {
      for( count=0; count<binom(n,k+1); count++ )
	{
	  // checks whether this lcm has already been found and updates structure
	  orbits->check(lcms[k][count],k,count);
	}
    }
  //orbits->show();

  first=orbits;

  
  //determine minimal and maximal degree in the first period of ch (rough estimates to get bounds for a table)
  //get period as well;
  d=1;
  mindeg=-2-2*d; maxdeg=4*n-6-2*d;
  period=-2*lcms[n-1][0];
  for( i=0; i<n; i++ )
    {
      mindeg+=2*((int) d/exponent[i]); maxdeg+=2*((int) d/exponent[i]);
      period+=2*(lcms[n-1][0]/exponent[i]);
    }




  for( d=2; d<=lcms[n-1][0]; d++ )
    {
      mindeg2=2*n-2-2*d; maxdeg2=4*n-6-2*d;
      for( i=0; i<n; i++ )
	{
	  mindeg2+=2*((int) d/exponent[i]); maxdeg2+=2*((int) d/exponent[i]);
	}
      if( mindeg2<mindeg )
	mindeg=mindeg2;
      if( maxdeg2>maxdeg )
	maxdeg=maxdeg2;
    }
   //cout << "mindeg/maxdeg: " << mindeg << " " << maxdeg << endl;
  mindeg-=10;
  maxdeg+=10;
  //normalize maxdeg
  maxdeg=maxdeg-mindeg;
  chrank=new int[maxdeg+1];
  //cout << "mindeg/maxdeg: " << mindeg << " " << maxdeg << endl;
  for( i=0; i<=maxdeg; i++ )
    {
      chrank[i]=0;
    }
  // consider time pi/2*d
  for( d=1; d<=lcms[n-1][0]; d++ )
    {
      //find largest orbit space which period divides d*pi/2
      k=0; count=0; found=0; maxlcm=1;

      while( orbits!=NULL )
	{
	  if( d%orbits->lcm==0 )
	    {
	      if( found )
		{
		  if( orbits->lcm>maxlcm )
		    {
		     
		      k=orbits->k;
		      count=orbits->number;
		      maxlcm=orbits->lcm;
		    }
		}
	      else
		{
		  
		  found=1;
		  k=orbits->k;
		  count=orbits->number;
		  maxlcm=orbits->lcm;
		}
	      	      
	    }
	  orbits=orbits->next;
	}
      
      //if d is divisible by one of the lcms, we get a contribution for homology
      if( found )
	{
	  //	  cout << d << ": " << k << " " << count << endl;
	  //increase the rank of the appropriate chgroups
	  degree=0;
	  for( i=0; i<n; i++ )
	    {
	      if( d%exponent[i]==0 )
		degree+=2*d/exponent[i];
	      else
		degree+=2*((int) d/exponent[i])+1;
	    }
	  //k+1 is the number of exponents: 2(k+1)-4 is dim orbit space
	  degree=degree-2*d+n-4-(k+1)+2;
	  for( i=0; i<=2*(k+1)-4; i+=2)
	    {
//	      cout << "Degree " << degree << " " << k << " " << count << endl;
	      chrank[degree+i-mindeg]+=1;
	    }
	  chrank[degree+(k+1)-2-mindeg]+=kappa[k][count];
	}
      orbits=first;

    }


  
  cout << "The mean Euler characteristic is ";
  count=0;
  for( i=0; i<=maxdeg; i++ )
    {
      if( (abs(i+mindeg))%2==0 )
	sign=1;
      else
	sign=-1;
      count+=chrank[i]*sign;
    }
  cout << (float) count/((float) fabs(period) ) << "=" << count << "/" << abs(period) << endl;


  cout << "Press enter to list some contact homology groups" << endl;
  cin.ignore();
  cin.ignore();



  // now list some ranks, simple version which is not always correct
  cout << "Contact homology is given by (only correct if period!=0 and no generators in degree -1,0 or 1" << endl << "The period is " << period << endl;
  for( i=0; i<=maxdeg; i++ )
    {
      count=chrank[i];
      k=i-period;
      while( k<=maxdeg && k>=0 && period!=0 )
	{
	  count+=chrank[k];
	  k-=period;
	}
      cout << "In degree " << i+mindeg << " ch has rank " << count << endl; 
    }


  cout << endl;

  
  delete chrank;
  
}