Maximize (L1+L2+..)!/L1!L2!..

Under restriction of Li>1 and L1+2L2+3L3..=N

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define NMAX 300
int N;
double lfactab[NMAX];

// log(n!)
double logfac(int n)
{
    double res=0;
    int i;
    for(i=2;i<=n;i++)res += log(i*1.0);
    return res;
}

// tabulize log(n!)
void init_lfac()
{
    int i;
    for(i=1;i<NMAX;i++)lfactab[i]=logfac(i);
}

// log((\sum_Lk)!/L1!L2!...)
double comb(int *nk, int dm)
{
    int i;
    int sum=0, sum2=0;
    for(i=1;i<=dm;i++)
    {
	sum+=nk[i];
#ifdef DEBUG
	sum2+=nk[i]*i;
#endif
    }
#ifdef DEBUG
    if(sum2!=N)
    {
	printf("#invalid\n");
	exit(1);
    }
#endif
    double c=lfactab[sum];
    for(i=1;i<=dm;i++)
	if(nk[i]>=2) c-=lfactab[nk[i]];
    return c;
}

// best values
int nrec=0;
#define MAXREC 10
int nkmax[MAXREC][NMAX+1];
double cmax=0;

// check if it is a new record
void check(int *nk, int dm)
{
    int i;
    double c=comb(nk, dm);
    if(c>cmax+1e-5)
    {
	cmax=c;
	nrec=0;
	for(i=0;i<=dm;i++)nkmax[nrec][i]=nk[i];
	for(i=dm+1;i<=N;i++)nkmax[nrec][i]=0;
    }
    else if(fabs(c-cmax)<1e-5)
    {
	nrec++;
	for(i=0;i<=dm;i++)nkmax[nrec][i]=nk[i];
	for(i=dm+1;i<=N;i++)nkmax[nrec][i]=0;
    }
#ifdef DEBUG
    printf("###:");
    for(i=1;i<=dm;i++)printf(" %d ",nk[i]);
    printf(" C=%f \n", exp(c));
#endif
}

// recursive loop
void rloop(int *nk, int lidx, int dm, int rest)
{
    int i;
#ifdef DEBUG
    if(rest<0)
    {
	printf("Err idx %d rest %d\n", lidx,rest);
	exit(1);
    }
    //debug
    printf("#rloop lidx=%d Rest %d [",lidx, rest);
    for(i=1;i<lidx;i++)printf("%d,",nk[i]);
    printf(" * ");
    for(i=lidx+1;i<=dm;i++)printf(" - ");
    printf("]\n");
#endif

    if(rest==0)
    {
	check(nk,lidx-1);
	return;
    }

    int nmax = rest/lidx;
    int ii;
    if(nmax==0 )
    {
	//rest!=0 no way to adj
	return;
    }
    if(nmax==1)
    {
	//no more to fit
	if(rest==lidx)
	{
	    // just fit
	    nk[lidx]=1;
	    check(nk,lidx);
	}
	return;
	//done
    }
    if(lidx==dm) return;

    for(ii=0;ii<=nmax;ii++)
    {
	nk[lidx]=ii;
        rloop(nk, lidx+1, dm, rest-ii*lidx);
    }
}

int main(int argc, char **argv)
{

    init_lfac();
    int i,j,k;
    //N=atoi(argv[1]);
    int nk[NMAX+1];
#define PMAX 20
    int justval[PMAX];
    for(i=0;i<PMAX;i++)
    {
	int nmax=i+1;
	int sum=0;
	int fac=1;
	for(j=nmax;j>=1;j--)
	{
	    sum += fac*j;
	    fac*=2;
	}
	justval[i]=sum;
	printf("# just %d = %d\n", i+1, justval[i]);
    }

    for(N=1;N<=100;N++)
    {
	int dm=0;
	while(justval[dm]<N)dm++;
	dm++;

	nk[1]=N;
	for(i=2;i<N;i++)nk[i]=0;
	check(nk,dm);
	rloop(nk, 1, dm, N);

	for(j=0;j<=nrec;j++)
	{
	    printf("%d  ",N);
	    int nmax=dm;
	    if(nkmax[j][dm]==0) nmax--;
	    for(i=1;i<=nmax;i++)printf("%d ",nkmax[j][i]);
	    printf(" #C=%f \n", exp(cmax));
	    fflush(stdout);
	}

    }
}