/*
	Set Class Generator
    Copyright (c) 2004 Nicholas M. Collins. All rights reserved.
	Written 24th May 2004.
	http://www.cus.cam.ac.uk/~nc272/

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

Set Class Generation Code to find the k-set representatives in n-space, for n<=64
with respect to <T,I> permutation subgroup
 
CAUTION- may take a long time to run and generate a large text file for n over 30

note that the output n-bit number is represented with the lowest n bits of the ULL- this is for the convenience of 
binary operations for transposition and inversion as utilised in testing canonicity.

If reading back the output text file of codes, read the least significant bit as index 0 of an array of 0s and 1s representing
the pitch class set (see the printpcset() function below)
If you want the 8-sets in Z12, this code will produce the 4-sets. Get the 8-sets by complementation and replacement of representatives

I use globals and avoid some C++ stuff for efficiency
to make this even faster you should do your own stack implementation perhaps with pointers to the data 
avoiding vtables and copy construction

Tested on the MAC OS X Xcode compiler.
This code is hopefully cross platform given that it only uses standard libraries 
BUT the behaviour of unsigned long long should always be checked for your machine! 
*/

#include <iostream>
#include <fstream>
#include <stack>
#include <stdlib.h>
#include <time.h>

using std::endl;

#define ULL unsigned long long

int n=12;
int k=6;

int digits;
int ones;
int penultimate;
	
ULL mask[64], mask2[64];

//debug testing/status
ULL candconsid=0;
int stackmax=0;
ULL found=0;

//make an array of masks before you begin to save time
void makemasks() {
int i;

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

mask[i]=0x1;
mask[i]=mask[i]<<i;
mask2[i]=mask[i]-1;
}

}


//test code. also shows how to get a pcset array representation from the binary number
void printpcset(ULL pcset)
{
std::cout << pcset << "     ";

for(int i=0; i<64; ++i) {
int out=0;
if (pcset & (mask[i])) out=1;
std::cout << out; 
}
std::cout << '\n';
}


//it's actually the inverse (lowest bit most significant) that gets passed in, so construct input set to test from the inverse 
bool testrepresentative(ULL inverse) {
	++candconsid;

	//return false as early as possible
	//generate set class, <T, I> initially
	//n bits of input important
	//only have to do transpositions that keep a 1 in leading place
	//only possibilities have minimal zero runs for most sig.bits in binary rep (could refine that further)

	ULL input=0;

	//construct inverse, then can transpose this around without recalculating
	int i;

	for(i=0; i<n; ++i)
	{
	if((inverse & (mask[i]))>0)
	input += mask[n-i-1];
	}

	//printpcset(input);
	//printpcset(inverse);

	//test inverse code against original
	if(inverse>input) return false;

	//brute force
	for(i=1; i<n; ++i)
	{

	//printpcset(input >> i);
	//printpcset(input & (mask2[i]));
	//printpcset((input & (mask2[i])) << (n-i));

	ULL code1 = (input >> i) | ((input & (mask2[i])) << (n-i)); 

	//printpcset(code1);

	if(code1>input) return false;

	code1 = (inverse >> i) | ((inverse & (mask2[i])) << (n-i)); 

	//printpcset(code1);

	if(code1>input) return false;

	}

	return true;
}


struct Candidate{
public:
ULL pcset;
int level;
int lastzero;

Candidate() {pcset=0; level=-1; lastzero=0;}

Candidate(ULL set, int lvl, int lz){pcset=set; level=lvl; lastzero=lz;}

Candidate(const Candidate& c){pcset=c.pcset; level=c.level; lastzero=c.lastzero;}

Candidate& operator=(const Candidate& c){pcset=c.pcset; level=c.level; lastzero=c.lastzero; return *this;}

};

//this is probably inefficient- use of Candidate* might help 
std::stack<Candidate> stack;

//make output file, gloabl to avoid passing a reference around
std::ofstream to("output.txt");	


//only add if canonical, return number considered
void addcandidates(Candidate& c) {

	int lev=c.level;
	int last1=c.lastzero;
	
	int left= digits-1-last1;
	int tofit= ones-lev;
	
	if(left<tofit) return;		//nothing to do
	
	//generate in order of increasing code 	
	for(int i=0; i<(left-tofit+1);++i) {

		ULL output = c.pcset;
		int ind= i+last1+1;
		output= output + mask[ind];
		
		//add back missing one at least significant bit
		ULL test= output;
		
		//std::cout <<"check return of least sig bit" << '\n';
		//printpcset(test);
		test=(test<<1) +1; 
		//printpcset(test);
		
		if(testrepresentative(test))
		//will make a temporary copy- need to make reference or pointer based for efficiency
		stack.push(Candidate(output,lev+1,ind));
	}
	
}


//write to file if canonical, return number considered
void testcandidates(Candidate& c) {
	int lev=c.level;
	int last1=c.lastzero;
	
	int left= digits-1-last1;
	int tofit= ones-lev;
	
	if(left<tofit) return;		//nothing to do
	
	//generate in order of increasing code 	
	for(int i=0; i<(left-tofit+1);++i) {
		
		ULL output = c.pcset;
		int ind= i+last1+1;
		output= output + mask[ind];
		
		//add back missing one at least significant bit
		ULL test= output;
		test=(test<<1) +1; 
		
		if(testrepresentative(test)) {
			//will make a temporary copy? Need to make reference or pointer based for efficiency
			to << test << " "<<endl;
			//print candidate
			//printpcset(test);
			++found;
		}
	}

}


//uses double calculation with running divisor to avoid overloads 
ULL combinations(double m, double r) {
	int i;
	double product=1;
	
	for(i=0; i<r;++i)
	product *= (m-i)/(r-i);
	
	return (ULL)product;
}



int main (int argc, char * const argv[]) {

	//testing code
	//unsigned long long x;
	//ULL x= 0x1<<30; //shift operations will fail if you're not careful
	//0x1 <<31 assumes a long, must cast to ULL before start right shifting
	//std::cout << x <<" "<< (x<<20) <<" "<<(0x1<<31) <<" "<< (0x1<<31<<20) <<" "<< (0x1<<19)*x <<" "<< 8*x<< " " << sizeof(x) << '\n'; //endl  
	//for(int i=0;i<64;++i) {
	//printpcset(mask[i]);
	//printpcset(mask2[i]);
	//}
	//test it works for big spaces
	//n=64;
	//testrepresentative(mask2[50]);

	//time this!
	clock_t starttime;
	starttime=clock(); //time(NULL);

	int prog,level;
	Candidate pivot;	

	//get na dn k from inputs
	n=atoi(argv[1]);
	k=atoi(argv[2]);
		
	std::cout <<"finding the "<< k <<"-sets in " << n <<"-space"<< std::endl;

	//prepare masks, for optimisation
	makemasks();
		
	//floor of n/2
	int maxk= n>>1;
	
	//stupid cases
	if(k==0) {return 0;}
	if(k==n) {to << mask[n-1];  return 0;}
	if(k==1) {to << mask[0]; return 0; }
	if(k==(n-1)) {to << mask[n-2];  return 0;}
	
	if(k>maxk) k=n-k;
	
	//only need to deal with numbers one less than k and n since can assume always trailing one
	//this will be put back in before any representative test
	digits=n-1;
	ones=k-1;
	penultimate=ones-1;
	
	//make stack using list and add and pop
	Candidate start(0,0,-1);
	stack.push(start);
	
	prog=0;
	
	//really underestimates for small n, many cases where wrong but hard to tell how many without having the answer already...
	//estimate number to find as 1% of all possible pcsets, so use (choose k from n), and compare to found
	ULL comb= combinations(n,k)/100;
	ULL modtest= comb/25; //25 reports
	if(modtest<5) modtest=5;
	ULL nexttest=modtest;
	//std::cout << combinations(n,k)<< " " <<comb<<" "<<modtest<<" "<<nexttest<<" " << std::endl;
	
	while(!(stack.empty()))
	{
		pivot= stack.top();
		stack.pop();
		
		level= pivot.level;
		
		//take these progress reports out for faster performance
		if(found>nexttest) {
			//estimate progess
			nexttest+=modtest;
			std::cout <<(100*((double)found/comb))<<"%"<<endl;
		}
		
		if(level==penultimate) {
		//std::cout <<endl<< "pivot" << std::endl;
		//printpcset(pivot.pcset);
		//std::cout <<"from this pivot " << std::endl;

			testcandidates(pivot);
		}
		else
		addcandidates(pivot);
	
		if((stack.size())>stackmax) {stackmax=stack.size();}
		
	}

	std::cout <<"found "<< found <<"considered "<< candconsid <<" stackmax" << stackmax << std::endl;

	//difftime(time(NULL),starttime)
	std::cout <<"took "<< (((double)(clock()-starttime))/CLOCKS_PER_SEC) <<" seconds "<< std::endl;

	//output file is global and will close on return
	
    return 0;
}