/*

A026820.c A070289.c A213008.c
Copyright (C) June, 2012 Sergei Viznyuk <phystech@hotmail.com>

This program takes positive integers n and k as inputs.
For the given values of n and k it prints the pair of numbers M,L
where M is the number of partitions  of the given integer n into <=k parts,
and L is the number of distinct values of multinomial
coefficients for the partitions of the given integer n into <=k parts

The number of partitions M of n into <=k parts forms OEIS sequences A026820.

The number of distinct values of multinomial coefficients L
forms OEIS sequence A070289, when k=n.

The number of distinct values of multinomial coefficients L
forms OEIS sequence A213008.

The algorithm for finding partitions of a given integer n is implemented
from the following paper:
Constant Time Generation of Integer Partitions
YAMANAKA et al. IEICE Trans Fundamentals.2007; E90-A: 888-895 
The example MATLAB code used:
http://www.mathworks.com/matlabcentral/newsreader/view_thread/162379

This code is limited by LONG DOUBLE precision arithmetics and
accuracy of built-in functions in math library. The code
compares logarithms of the value of multinomial coefficient, and decides
if the values are the same or different to within 1/10000000 precision.
That is to say, for sufficiently large values of n the difference between
values of multinomial coefficients may be undetectable by this code and requires
infinite precision arithmetics.

This 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.

*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#ifdef __CYGWIN__
#define EXP exp
#define LOG log
#define LGAMMA lgamma
#define POW pow
#define FABS fabs
#else
#define EXP expl
#define LOG logl
#define LGAMMA lgammal
#define POW powl
#define FABS fabsl
#endif

typedef struct { long double E,D; long unsigned k,m[]; } wt;

static wt **s;
static long unsigned k,K,cr=0,inc=1000;
static long double dn,en;

static void cmpt(wt *a)
{ long unsigned i=0,j=0,k=a->k,*m=a->m;
  long double d=dn, e=en;
  while ( i<k )
    { if ( m[i] > 1 ) e-=LGAMMA(1.0L+m[i]);
      if ( m[i]!=m[j] ) { if ( i>j+1 ) d-=LGAMMA(1.0L+i-j); j=i; }
      i++;
    }
  if ( i>j+1 ) d-=LGAMMA(1.0L+i-j);
  if ( K>i+1 ) d-=LGAMMA(1.0L+K-i);
  a->E=e;
  a->D=d;
}

static int compar(wt **ap,wt **bp)
{ wt *a=*ap, *b=*bp;
  if ( !a->E ) cmpt(a);
  if ( !b->E ) cmpt(b);
  if ( a->E > b->E ) return 1;
  if ( a->E < b->E ) return -1;
  return 0;
}

static void children3(wt *a)
{ wt *b;
  long unsigned i=a->k,*m=a->m;
  if ( cr==inc )
    s=realloc(s,(inc+=inc)*sizeof(wt *));
  s[cr++]=a;
  if ( *m > m[1] )
    { if ( i < k )
	{ b=malloc(sizeof(wt)+(i+1)*sizeof(long unsigned));
	  memcpy(b,a,sizeof(wt)+i*sizeof(long unsigned));
	  (*b->m)--; b->m[i]=1; b->k++; children3(b);
	}
      if ( m[i-2] > m[i-1] && ( i > 2 || (m[i-2]-m[i-1]) > 1 ) )
	{ b=malloc(sizeof(wt)+i*sizeof(long unsigned));
	  memcpy(b,a,sizeof(wt)+i*sizeof(long unsigned));
	  (*b->m)--; b->m[i-1]++; children3(b);
	}
    }
}

int main(int argn,char *argc[])
{ wt *a;
  long double d,e,r,q,f,p=0.0L;
  long unsigned i,N,L=1;
  if ( argn!=3 )
    { fprintf(stderr,"Usage: %s <n> <k>\n",argc[0]); return -1; }
  N=atol(argc[1]);
  k=K=atol(argc[2]);
  if ( N < k ) k=N;
  if ( !k ) return 0;
  r=LOG(1.0L*K);
  q=EXP((1.0L-N)*r);
  f=N*r;
  en=LGAMMA(1.0L+N);
  dn=LGAMMA(1.0L+K);
  if ( k==1 ) goto ne;
  s=malloc(inc*sizeof(long unsigned *));
  a=malloc(sizeof(wt)+2*sizeof(long unsigned));
  a->E=0.0L; a->k=2; *a->m=N-1; a->m[1]=1;
  children3(a);
  if ( cr > 1 )
    qsort(s,cr,sizeof(long unsigned *),
    (int (*)(const void *,const void *))&compar);
  else
    compar(&a,&a);
  for (i=0;i<cr;i++)
    { a=s[i]; e=a->E; d=a->D;
      q+=EXP(d+e-f);
      if ( p<10000000.0L*FABS(p-e) )
	{ L++; p=e; }
      free(a);
    }
  free(s);
ne:
  fprintf(stdout,"%lu,%lu\n",cr+1,L);
  if ( q>1000000000000.0L*FABS(q-1.0L) ) return 0;
  fprintf(stderr,"*** Checksum error: q=%Le\n",q);
  return 1;
}