/************************************************
  This file contains the definition of a class "partit_base"
  that represents the number of partitions which divide an integer number of "n"
  indistinguishable elements into "m" non-empty subsets. The member
  functions next() constructs the sequence of all possible partitions, returned
  as a vector<int> of size "m" upon repeated calls.

  Based on this, there is another class "partit_unrest" of unrestricted partitions
  of "n" into 1 to n non-empty subsets, which are constructed by working
  through all "partit_base" instantiations with 1 to n subsets.

  Two possible subtypes are defined: the standard one allows that terms ni
  occur more than once in the decomposition n= n1 + n2 + n3 + .. + nm,
  allowing for a decomposition 5=1+1+3, for example, the more restricted
  subtype includes only decompositions with all different ni, say 5=2+3
  but not 5=1+1+1+2.
 
  Both classes contain a member function size() that returns the total count
  of partitions to be created by calling next(). For the class partit_unrest, this
  equals that tabulation in chapter of the "Handbook of Mathematical Functions"
  by Milton Abramowitz and Irene A Stegun.

  It also contains a simple representation of Stirling numbers of the
  second kind, and an primitive implementation of a prime number class.

  An example of use is contained in the final main().

  See also
  http://www.theory.csc.uvic.ca/~cos/gen/nump.html

  Richard J Mathar, http://www.strw.leidenuniv.nl/~mathar, 27 Aug 2001
***********************************************************/

#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <iostream.h>
// <map>, <vector> as contained in the STL
#include <map>
#include <vector>

using namespace std ;

// Trivial implementation of the factorial. n!=1*2*3*...*n; 0! =1
double factorial(const int n)
{
	double resul=1. ;
	for(int k=1; k<=n; k++)
		resul *= k ;
	return resul ;
}

// Stirling Number of Second Kind. m<=n. Often very large and
// therefore returned and calculated as  a double, not an int. Returns -1 on invalid (n,m) parameter combinations.
// Equation in section C of M Abramowitz and I A Stegun chapter 24.1.4.
// For example Stirling2(7,3)=301. Stirling2(8,5)=1050.
double Stirling2(const int n, const int m)
{
	if( m==0)
		return !n ;
	if( m>n || n<1 || m < 1)
		return -1. ;
	if( m==1 || n==m)
		return 1. ;
	double resul=0. ;
	for( int k=1 ; k<=m ; k++)
	{
		const double tmpk= pow(double(k),double(n))/factorial(m-k)/factorial(k) ;
		if( (m-k) %2 )	// odd m-k
			resul -= tmpk ;
		else		// even m-k
			resul += tmpk ;
	}
	return resul ;
}

enum partit_type {
	SAME_SIZE ,	// same (repeated) numbers allowed in the partitioning, 5=1+1+3 counts as a partition
	DIFF_SIZE ,	// no number allowed more than once in the partition terms, 5=2+2+1 *not* a partition, but 5+2+3 is.
	INVALID		// not a valid partitioning (invalid parameter combination at instatiation time)
} ;


// compute an integer ceil(numer/denom) which is the smallest integer equal to or larger than the fraction numer/denom.
// It is based on the fact that C/C++ rounds towards zero for the integer division (both signs of the intermediate fraction).
int iceil(const int numer, const int denom)
{
	if( numer*denom > 0 )		// positive numer/denom
		return (numer+denom-1)/denom ;
	else				// negative numer/denom, or numer==0
		return numer/denom ;
}

			// any sufficiently negative number here....
#define PARTIT_BASE_NO_MAXCEIL -9999
// A class representing the number of partitions of 'n' into 'm' subsets
// without respect to order. The total count of all possible partitions is returned with the size() function.
// If n=7, m=3, SAME_SIZE, we get size()=4 decompositions of 7 into 3 subsets. 2+2+3=1+3+3=1+2+4=1+1+5
// If n=7, m=3, DIFF_SIZE, we get size()=1 decompositions of 7 into 3 subsets. 1+2+4: the others had repeated 2's, 3's or 1's
// If n=7, m=3, SAME_SIZE, maxnels=4 with the creator we get size()=3 decompositions of 7 into 3 subsets. 2+2+3=1+3+3=1+2+4
//	as the 1+1+5 decomposition had one element 5 that exceeded 'maxnels'.
// See M Abramowitz- I Stegun chapter 24.1.4.
class partit_base {
	private :
	partit_base *child ;	// pointer to the 'child' partitions with m-1 subsets
	int maxceil ;		// upper limit allowed for any number in the terms of the decomposition of *this

	public:
	int n ;		// number of elements
	int m ;		// number of subsets
	int maxel ;	// maximum number of elements in the subsets, 1<=maxel<=n
	enum partit_type kind; 	// SAME_SIZE or DIFF_SIZE

	// constructor
	// The optional third argument sets a limit to the maximum size of any of the subsets
	// (which otherwise - in the case of a negative default - is taken as the logical default of n-m+1.)
	partit_base(const int nelements, const int nsubsets, enum partit_type typ = SAME_SIZE,
			const int maxnels = PARTIT_BASE_NO_MAXCEIL)
			: n(nelements), m(nsubsets), maxceil(maxnels), kind(typ)
	{
		// The representation of all possible partitions is done recursively through
		// {n1,n2,n3,...nm}= loop over [maxel = ceil(n/m)...n-m+1] of {n1',n2',n3',..,nm-1,maxel}
		// where ceil() denotes the smallest integer greater or equal to the argument, and where the vector of the
		// n1+n1+n2+...+nm=n number of elements is sorted arbitrarily such that the max(ni), i=1..m, comes last,
		// and where {n1',n2',n3',...,nm-1'} is a "child"  partition with m-1 sets of elements (with a maxel'<=maxel
		// number in the final position).
		// ceil(n/m) can be represented by the integer division (n+m-1)/m here:
		// (n=40,m=7)->5,... -> 6; (n=42,m=7) -> 6 -> 6, (n=43,m=7)-> 6,... -> 7
		switch ( kind )
		{
		case SAME_SIZE :
			if( maxceil == PARTIT_BASE_NO_MAXCEIL || maxceil > n-m+1)	// effectively no constraint set
				maxceil= n-m+1 ;
				// compute the smallest integer maxel such that maxel+...+maxel+maxel>= n
			maxel= iceil(n,m) ;
			break ;
		case DIFF_SIZE :
					// 1+2+..+(m-1)=m(m-1)/2, minimum sum of m-1 elements
			if( maxceil == PARTIT_BASE_NO_MAXCEIL || maxceil > n-(m*m-m)/2 )
				maxceil= n-(m*m-m)/2 ;
				// compute the smallest integer maxel such that (maxel-m+1)+(maxel-m+2)+...+(maxel-1)+maxel>= n
				// which means m*maxel-m*(m-1)/2>=n which means maxel >= n/m +(m-1)/2 = [2n+m(m-1)]/(2m)
			maxel= iceil(2*n+m*(m-1),2*m) ;
			break ;
		case INVALID :
			break ;
		default :
			cerr << "partit_base(): unknown kind " << kind << endl ;
			exit(EXIT_FAILURE) ;
		}

		// detect some more invalid parameter combinations
		if( n <= 0 || m <=0 || maxceil <= 0  || maxel > maxceil )
			kind = INVALID ;
		
		switch ( kind )
		{
		case SAME_SIZE :
			child = (m > 1 && maxel >= 1)? new partit_base(n-maxel,m-1,kind,maxel) : 0 ;
			break ;
		case DIFF_SIZE :
			child = (m > 1 && maxel >= 2)? new partit_base(n-maxel,m-1,kind,maxel-1) : 0 ;
			break ;
		case INVALID :
			child = 0 ;
		}
	}

	// destructor
	~partit_base()
	{
		delete child ;
	}

	// iterator; return the next partition in the sequence, or a vector of zero size() if all have been returned.
	// partit_base(n,m) is selected from set of all vectors of length m-1 of partit_base(n-maxel,m-1), where
	// ceil(n/m)<=maxel<=n-m-1.
	// This recursive scanning means for example: look at all partitions of n=40 and m=7 by looking (in sequence)
	// at the "child" partitions (n',m')=(34,6)+6, (33,6)+7, (32,6)+8, ...(6,6)+34
	vector<int> next()
	{
		vector<int> resul ;		// the resulting vector of size 'm' or 0
		if ( kind != SAME_SIZE && kind != DIFF_SIZE)
			return resul ;		// zero-length vector
		if( child)
			resul = child->next() ;		// recursive attempt to count forward in the sub-partition
		else if ( m==1 )		// the m=1 case means no further partitions into m-1 sets, and maxel=n
						// (vector of length 1) to be returned first, then no further results.
		{
			resul.clear() ;
			if( maxel++ <= maxceil)
				resul.push_back(n) ;
			return resul ;		// either vector with one element (equal to n=maxel) or empty vector
		}
		else
			return resul ;		// empty vector

		// now resul may contain a 0-size vector (which means that there are no more child partitions with the current
		// maxel) or a m-1-size vector (which indicates a regular child partitions has been returned).
		if ( resul.size() )	// successful delivery of child->next()
			resul.push_back(maxel) ;	// just append to have an m-size vector
		else			// attempt to increase maxel to its maximum
		{
			delete child ;
			child=0 ;
			maxel ++ ;	// next generation children
			if( maxel <= maxceil )	// still elements left in the (virtual) loop over maxel
			{
				if( kind == SAME_SIZE)
					child = (maxel >= 1) ? new partit_base(n-maxel,m-1,kind,maxel) : 0 ;
				else
					child = (maxel >= 2) ? new partit_base(n-maxel,m-1,kind,maxel-1) : 0 ;
				return next() ;	// recursive call with a new sub-partition
			}
		}
		return resul ;
	}

	// Number of non-zero solutions to the decomposition problem.Often very large and
	// therefore returned and calculated as  a double, not an int. Returns 0 on invalid (n,m) parameter combinations.
	double size () const
	{
		if( m==0) 	// ????
			return !n ;
		if( kind != SAME_SIZE && kind != DIFF_SIZE )
			return 0. ;
		if( m==1 || n==m)
			return 1. ;
		double resul=0. ;
			// use some local copies of 'maxel' and 'child' not to disturb the potential calls to next()
		int lastel ;
		switch (kind)		// code below same criteria as with the constructor before...
		{
		case SAME_SIZE:
			lastel = iceil(n,m) ;
			break ;
		case DIFF_SIZE:
			lastel= iceil(2*n+m*(m-1),2*m) ;
			break ;
		}
		for(  ; lastel <= maxceil ; lastel++)
		{
			partit_base *cntchil ;
			if( kind == SAME_SIZE)
				cntchil = (lastel >=1 ) ? new partit_base(n-lastel,m-1,kind,lastel) : 0 ;
			else
				cntchil = (lastel >= 2)? new partit_base(n-lastel,m-1,kind,lastel-1) : 0;
			if ( cntchil)
				resul +=  cntchil->size() ;
			delete cntchil ;
		}
		return resul ;
	}
} ;

// A class representing all partitions of 'n' into non-empty subsets
// without respect to order. The total count of all possible partitions is returned with the size() function.
// See M Abramowitz- I Stegun chapter 24.2.1.
// Example: 5=1+4=2+3=1+1+3=1+2+2=1+1+1+2=1+1+1+1+1 for n=5 and kind=SAME_SIZE; yields size()=7.
// Example: 5=1+4=2+3 for n=5 and kind=DIFF_SIZE; yields size()=2.
class partit_unrest {
	private :
	partit_base *child ;	// pointer to the 'child' partitions with m subsets

	public:
	int n ;		// number of elements
	int m ;		// number of subsets
	int maxm ;	// maximum number of elements in the subsets, 1<=maxel<=n
	enum partit_type kind; 

	// constructor
	partit_unrest(const int nelements, enum partit_type typ = SAME_SIZE) : n(nelements) , kind(typ)
	{
		// The representation of all possible partitions is a loop over the "restricted" partitions
		// with m=1,...,n subsets (terms in the sum). In the template case n=5: m=1 yields '5', m=2 yields 1+4 and 2+3,
		//	m=3 yields 1+1+3 and 1+2+2, m=4 yields 1+1+1+2, and m=5 yields 1+1+1+1+1
		m=1 ;
		switch( kind)
		{
		case SAME_SIZE :
			child = (n >= m) ? new partit_base(n,m,kind) : 0;
			maxm = n ;
			break ;
		case DIFF_SIZE :
			child = ( n >= (m*m+m)/2 ) ? new partit_base(n,m,kind) :0 ;
			maxm = -0.5+sqrt(0.25+2.*n) ;	// ensures that 1+2+...+m=m(m+1)/2 <= n
			break ;
		default:
			cerr << "partit_unrest(): unknown kind " << kind << endl ;
			exit(EXIT_FAILURE) ;
		}
	}

	// destructor
	~partit_unrest()
	{
		delete child ;
	}

	// Iterator function; return the next partition in the sequence, or a vector of zero size() if all have been returned.
	// partit_base(n,1) is selected first, then the elements with m=2, then with m=3 and so on.
	vector<int> next()
	{
		vector<int> resul ;		// the resulting vector of size 1..'n', or 0 if the subspace is exhausted
		if( child)
			resul = child->next() ;		// recursive attempt to count forward in the sub-partition
		// now resul may contain a 0-size vector (which means that there are no more child partitions with the current
		// m) or a m-size vector (which indicates a regular child partition has been returned).
		if ( resul.size() )	// successful delivery of child->next()
			return resul ;
		else			// attempt to increase m to its maximum
		{
			delete child ;
			child=0 ;
			m ++ ;	// next generation children
			if( m <= maxm)	// still elements left in the (virtual) loop over maxel
			{
				child = new partit_base(n,m,kind) ;
				return next() ;	// recursive call with a new sub-partition
			}
		}
		return resul ;
	}

	// Number of non-zero solutions to the decomposition problem. Often very large and
	// therefore returned and calculated as  a double, not an int. Returns -1 on invalid (n,m) parameter combinations.
	// Equals p(n) in the notation of Abramowitz, Stegun, Table 24.5,
	double size () const
	{
		double resul=0. ;
				// short lookup tables p(n) and q(n) for smallest 'n', starting at n=1, not n=0...
		const int p[]={1,2,3,5,7,11,15,22,30,42,56,77,101,135,176,231,297,385,490,627,792,1002,1255} ;
		const int q[]={1,1,2,2,3,4,5,6,8,10,12,15,18,22,27,32,38,46,54,64,76,89,104,122,142} ;
		if( kind == SAME_SIZE && n <= sizeof(p)/sizeof(int) )		// n  found in the lookup table
			resul= p[n-1] ;
		else if( kind == DIFF_SIZE && n <= sizeof(q)/sizeof(int) )		// n found in the lookup table
			resul= q[n-1] ;
		else if( kind == INVALID )
			;
		else
		{
			for( int subm = 1 ; subm <= maxm ; subm++)
			{
				partit_base *cntchil = new partit_base(n,subm,kind) ;
				resul +=  cntchil->size() ;
				delete cntchil ;
			}
		}
		return resul ;
	}
} ;
#undef PARTIT_BASE_NO_MAXCEIL

// for purposes of sorting map<int,int,less_int>
struct less_int {
	bool operator()(const int &i1, const int &i2) const
	{
		return i1 < i2 ;
	}
} ;

// prime number class: contains the first prime numbers in a sorted vector
class prime : public vector<int> {
	private:
	// Extend the current list until all primes less or equal to 'maxpr' are covered. This uses the sieve
	// of Erastothenes to extend the list of primes from the currently stored numbers up to twice the maximum
	// of these. The restriction to the twice the maximum stems from the simple way of constructing the new
	// numbers: a full list of integers is produced in a map<int,int,less_int> first, then integer multiples  of
	// the numbers already in the prime number list are removed from the map<>. If map<> would contain larger numbers
	// than twice the ones known in the present prime number list, those could be integer multiples of prime numbers
	// in the map<> but not in the present list, which would evade being eliminated....
	void sieve()
	{
		int sz= size() ;		// count of known primes
		int c= (*this)[sz-1] ;		// largest prime known (=in the list)
		const int maxpr = 2*c ;		// largest number in the new list of prime number candidates
		map<int,int,less_int> candid ;		// candidates to be inserted: key=the prime number, value=don't care
				// start with a set of candidates : all integers in the  range from the currently known
				// largest prime up to 'maxpr'
		for(c= (*this)[sz-1]+1 ; c <= maxpr ; c++)
			if( c % 2)		// pre-select  and omit even numbers (those are not primes)
				candid.insert(map<int,int,less_int>::value_type(c,c)) ;
		c= (*this)[sz-1] ;		// largest prime known (=in the list)
		// sieve of Erastothenes to get rid of the keys in the map that are not primes
		for(int pr=0; pr< sz ; pr++)
		{
			const int strtpr = (*this)[pr] ;	// multiples of the prime 'strtpr' to be elimiated from 'candid'
			int mul=c/strtpr ,			// could start with mul=1, but need only to test for numbers > c
			    tst ;
			while( (tst = mul*strtpr) <= maxpr )	// test all multiples of numbers in the set of known primes
			{
				candid.erase(tst) ;		// no need to check whether present or not...
				mul++ ;				// to next higher multiple of strtpr
			}
		}
		// re-insert the survivors of the sieve into the vector
		map<int,int,less_int>::iterator cptr=candid.begin() ;
		while( cptr != candid.end() )
		{
			push_back( (*cptr).first ) ;
			cptr ++ ;
		}
	}
	public:

	// default constructor; start with a small, complete set of known primes
	prime()
	{
		static const int lookup[]={2,3,5,7,11,13,17,19,23,29,31,37,41,43} ;
		clear() ;
		for(int i=0; i < sizeof(lookup)/sizeof(int) ; i++)
			push_back(lookup[i]) ;
	}

	// grow until all primes up to at least 'maxpr' are in the list.
	void grow(const int maxpr)
	{
		for(;;)
		{
			int sz=size() ;		// number of known primes
			int c = (*this)[sz-1] ;	// largest known prime
			if( c < maxpr)
				sieve() ;	// grow list until covering the list of primes 2,3,5,7,11...,<= 2*c
			else
				break ;
		}
	}

	// simple test whether an integer (input 'number') is prime: returns true if prime.
	bool isprime(const int number)
	{
		// Compile all primes up to about sqrt(number). This ensures that testing (number|prime_in_the_list)
		// for all prime_in_the_list gives a definite answer
		grow( int(sqrt(number)) ) ;
		const int sz=size() ;		// number of known primes
		const int c = (*this)[sz-1] ;	// largest known prime
		if( number <= c)
		{
			vector<int>::iterator i = find(begin(),end(),number) ;	// search for number in current list of primes
			if ( i != end() )	// found it: is a prime
				return true ;
			else			// didn't find it: not a prime
				return false ;
		}
		else		// 'number' not covered by the number range of the present prime number list
		{
			for(int i= 0 ; i < sz ; i++)
				if( (number % (*this)[i]) == 0 )	// number is divisible by a member of the list
									// (but itself not in the list): not a prime
					return false ;
			return true ;
		}
	}

	// decompose an integer into prime factors; does n't yet work for n=1, as 1 is not in the list of primes.
	// n= key0^value0 * key1^value1 * key2^value2...
	map<int,int,less_int> decomp(int n)
	{
		grow(n) ;	// increase list of known prime numbers if necessary
		int prindx =0 ;	// index into the prime number vector table
		map<int,int,less_int>resul ;
		while( n >= 2)
		{
			int pwr = 0 ;	// power to be computed
			while( (n % (*this)[prindx]) == 0 ) 	// n divisible by the current prime
			{
				pwr++ ;
				n /= (*this)[prindx] ;
			}
			if( pwr)	// [prindx] was a (multiple) divisor: store in the result table
				resul.insert( map<int,int,less_int>::value_type((*this)[prindx],pwr)) ;
			prindx++ ;	// move over to next larger prime
			// if ( (*this)[prindx] > n) break ;
		}
		return resul ;
	}
} ;

// Test function: example of an application main()
ostream & operator<<(ostream &os, const vector<int> & somevect)
{
        copy(somevect.begin(),somevect.end(),ostream_iterator<int,char>(os," ")) ;
        return os ;
}

int main(int argc, char *argv[])
{
	if(argc != 3)
	{
		cerr << "usage: " << argv[0] << " no_elements no_subsets\n" ;
		exit(EXIT_FAILURE) ;
	}
	int n=atoi(argv[1]) ;
	int m=atoi(argv[2]) ;
	partit_base part(n,m) ;	// create/initialize partitioning of n elements into m subsets, type SAME_SIZE
	double count=part.size() ;
	cout << "restricted partitions type " << part.kind << endl ;
	cout << "expected number of elements: size(" << n <<"," << m << ")=" << count << endl ;
	vector<int> decomp;
	count=0. ;
	do {
		decomp=part.next() ;
		if( decomp.size() )
		{
			count++ ;
			cout << decomp << endl ;
		}
	} while ( decomp.size() ) ;
	cout << "constructed number of elements: " << count << endl ;

	partit_base part2(n,m,DIFF_SIZE) ;	// create/initialize partitioning of n elements into m subsets
	count=part2.size() ;
	cout << "restricted partitions type " << part2.kind << endl ;
	cout << "expected number of elements: size(" << n <<"," << m << ")=" << count << endl ;
	count=0. ;
	do {
		decomp=part2.next() ;
		if( decomp.size() )
		{
			count++ ;
			cout << decomp << endl ;
		}
	} while ( decomp.size() ) ;
	cout << "constructed number of elements: " << count << endl ;

	partit_unrest part3(n) ;	// create/initialize partitioning of n elements into any number of subsets, type SAME_SIZE
	cout << "unrestricted partitions type " << part3.kind << endl ;
	count=part3.size() ;
	cout << "expected number of elements: size(" << n << ")=" << count << endl ;
	count=0. ;
	do {
		decomp=part3.next() ;
		if( decomp.size() )
		{
			count++ ;
			cout << decomp << endl ;
		}
	} while ( decomp.size() ) ;
	cout << "constructed number of elements: " << count << endl ;

	partit_unrest part4(n,DIFF_SIZE) ;	// create/initialize partitioning of n elements into any number of subsets
	cout << "unrestricted partitions type " << part4.kind << endl ;
	count=part4.size() ;
	cout << "expected number of elements: size(" << n << ")=" << count << endl ;
	count=0. ;
	do {
		decomp=part4.next() ;
		if( decomp.size() )
		{
			count++ ;
			cout << decomp << endl ;
		}
	} while ( decomp.size() ) ;
	cout << "constructed number of elements: " << count << endl ;

	n= 200 ;
	cout << "primes up to " << n << endl ;
	prime prs ;
	prs.grow(n) ;		// collect list of primes up to 'n'
	cout << prs << endl ;	// print them
	n = 80 ;
	cout << n << " a prime: " << prs.isprime(n) << endl ;
	n = 1451 ;
	cout << n << " a prime: " << prs.isprime(n) << endl ;
	cout << prs << endl ;	// print them, containing e.g. 1451, 1453, 1459, 1471

	return 0 ;		// to please SunOS shells
}