/* 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; }