/* lpicks.c -- routines for assorting picks
 * 
 * int	n, N;
 * int	sum;
 * int	vals[N];
 * long	l, g, e;
 * l = lpicks(n, N, sum, vals);
 * g = gpicks(n, N, sum, vals);
 * e = epicks(n, N, sum, vals);
 *
 * lpicks returns the number of ways to pick n items from a total of N
 * such that their sum is less than sum.  The values are passed to lpicks
 * in vals.  gpicks and epicks are similar but return the
 * number of picks with sums > or = sum respectively.  
 * 
 * The return values from these routines can become very large, even for
 * reasonable values of N, and may overflow the integer representation.
 * In this case the following routines can be used:
 *
 * int		n, N;
 * int		sum;
 * int		vals[N];
 * double	ll, lg, le;
 * ll = llpicks(n, N, sum, vals);
 * lg = lgpicks(n, N, sum, vals);
 * le = lepicks(n, N, sum, vals);
 *
 * These return the natural logarithm of the same values as are returned
 * by the former routines.  If the number of picks satisfying the
 * condition is 0, -1.0 is returned.  Otherwise these routines will
 * never return a negative value.
 */
#include <stdio.h>
#include <math.h>
#include "common.h"
#include "logs.h"

#define	USELOGS	32			/* where pick can get big	*/

typedef struct AVF {			/* adjusted vals, frequencies	*/
    unsigned	val;
    unsigned	freq;
    unsigned	n;			/* # picked			*/
    unsigned	cf;			/* cumulative frequencies	*/
} AVF;

char		*ealloc();		/* allocate from the heap	*/
char		*erealloc();		/* reallocate space		*/
void		qsort();		/* quick sort			*/
long		pick();			/* number of distinct picks	*/
double		logpick();		/* log(pick)			*/
double		log();			/* natural logarithm		*/
double		exp();			/* e^x				*/
double		log1p();		/* log(1+x)			*/
double		fabs();			/* |x|				*/
void		dprint();		/* print debug info		*/
static void	massage();		/* put data in optimal form	*/
static int	cmp();			/* comparison for sort		*/
static int	argssame();		/* same args as last call?	*/
static void	legpicks();		/* collect all stats		*/
static void	lecount();		/* count picks <= ts		*/
static void	llegpicks();		/* callect logs of stats	*/
static void	llecount();		/* gets logs of counts		*/
static void	chkpick();		/* check that picks are OK	*/
static void	chkdiff();		/* check that increase right	*/

extern int	cflag;			/* check flag			*/
extern int	dflag;			/* debug flag			*/
extern int	indent;			/* indent for debug output	*/

static int	n = 0;			/* saved arguments		*/
static int	N = 0;
static int	sum = 0;
static int	*vals = NULL;

static unsigned	nv;			/* # vals			*/
static AVF	*avf;			/* adjusted vals, freqs		*/
static unsigned	ndv;			/* # distinct vals		*/
static unsigned	*cs;			/* cumulative sums		*/
static unsigned	ts;			/* target sum			*/
static long	lsum = 0L;		/* # < sum			*/
static long	esum = 0L;		/* # = sum			*/
static long	gsum = 0L;		/* # > sum			*/
static double	llsum = LOG0;		/* logs of picks		*/
static double	lesum = LOG0;
static double	lgsum = LOG0;
static double	csum = LOG0;		/* chkpick's sum		*/

long
lpicks(an, aN, asum, avals)
int		an;
int		aN;
int		asum;
int		*avals;
{
    if (!argssame(an, aN, asum, avals)) {
	legpicks();
    }
    return(lsum);
}

long
epicks(an, aN, asum, avals)
int		an;
int		aN;
int		asum;
int		*avals;
{
    if (!argssame(an, aN, asum, avals)) {
	legpicks();
    }
    return(esum);
}

long
gpicks(an, aN, asum, avals)
int		an;
int		aN;
int		asum;
int		*avals;
{
    if (!argssame(an, aN, asum, avals)) {
	legpicks();
    }
    return(gsum);
}

/* argssame -- return TRUE iff current args are same as from last call	*/
static int
argssame(an, aN, asum, avals)
int		an;
int		aN;
int		asum;
int		*avals;
{
    register int	*v;
    register int	i;
    int			same;

    v = (int *) ealloc(aN * sizeof(int));
    for(i=0; i<aN; i++) v[i] = avals[i];
    qsort(v, aN, sizeof(int), cmp);
    same = (an == n) && (aN == N) && (asum == sum) && (NULL != vals);
    if (same) {
	for(i=0; i<N; i++) {
	    if (v[i] != vals[i]) break;
	}
	same = (i >= N);
    }
    n = an;
    N = aN;
    sum = asum;
    if (NULL != vals) free(vals);
    vals = v;
    return(same);
}

static int
cmp(i1, i2)
int	*i1, *i2;
{
    return(
	(*i1 < *i2) ? -1 : (
	    (*i1 > *i2) ? 1 : 0
	)
    );
}

/* legpicks -- collect stats for picks
 * 
 * legpicks is called to determine the number of picks less than, equal
 * to, and greater than the specified sum.  The arguments n, N, sum,
 * and vals (the latter sorted from smallest to largest) are assumed to
 * have been copied to the static storage areas by argssame.  legpicks
 * leaves the number of picks less than, equal to, or greater than sum
 * in lsum, esum, and gsum.
 *
 * legpicks is mainly concerned with massaging the problem into the most
 * efficiently solvable problem.  It first adjusts the data values so
 * that all are non-negative by adding a number adj to them.  The target
 * sum, ts, becomes sum+n*adj.  Second, it arranges to use the smaller
 * of n and N-n.  That is, if sv is the sum of all the values, every
 * pick of n values that sums to less than ts coresponds to a pick of
 * N-n values that sums to >sv-ts.  If N-n is less than n, it pays to
 * solve this smaller problem instead.  nv is left equal to n or N-n,
 * whichever is smaller, nsw is TRUE if nv is N-n, and ts is set to
 * sv-ts if appropriate.  Finally, legpicks attempts to estimate
 * whether more picks will sum to greater than ts, or to less than ts.
 * In the first case, it inverts the problem by replacing each value v
 * with max-v, where max is the maximum value.  sv becomes max*N - sv,
 * and ts becomes max*nv - ts.  Since the actual counts are always made
 * by searching for just those picks that add up to <= ts, and then
 * computing gsum at the end, this inversion of the problem speeds up
 * the search.
 *
 * After massaging the problem statement, legpicks processes the data to
 * obtain two new arrays useful in gathering stats.  cs[i] is set equal
 * to Sum(v[j], {j, 0, i}).  Since the v[i] are sorted, cs[i] is the
 * minimum value a pick of i elements can add up to.  ndv is set equal
 * to the number of distinct values, and avf is set up as an array of
 * ndv value-frequency pairs.  These are sorted by value.
 *
 * Finally, legpicks calls lecount to count the number of picks that sum
 * up to less than ts (lsum), and the number that sum up equal to ts
 * (esum).  gsum is set to pick(N, n) - lsum - esum.  Depending on
 * previous rearrangements, gsum and lsum are switched to be appropriate
 * for the original problem.
 */
static void
legpicks()
{
    register int	i;
    register unsigned	u;
    register unsigned	*v;		/* scratch value array		*/
    int			adj;		/* value adjustment		*/
    unsigned		sv = 0;		/* summed values		*/
    int			nsw;		/* n has been switched		*/
    int			inv;		/* values have been inverted	*/

    indent = 0;
    dprint("n = %d, N = %d, sum = %d\n", n, N, sum);
    if (dflag) {
	for(i=0; i<N; i++) {
	    dprint("vals[%d] = %d\n", i, vals[i]);
	}
    }
    if (cflag) csum = LOG0;
    v = (unsigned *) ealloc(N * sizeof(unsigned));
    /*
     * adjust values to all >= 0
     */
    adj = -vals[0];
    for(i=0; i<N; i++) {
	v[i] = vals[i] + adj;
	sv += v[i];
    }
    if (sum + n*adj < 0) {
	lsum = esum = 0;
    }
    else {
	ts = sum + n*adj;		/* target sum */
	/*
	 * use the smaller of n and N-n (becomes nv)
	 */
	if (nsw = (N-n < n)) {
	    nv = N - n;
	    ts = sv - ts;
	}
	else {
	    nv = n;
	}
	/*
	 * invert the values if mostly smaller than ts/nv
	 */
	if (inv = ((nv * (v[N/2] + v[(N-1)/2])) < 2*ts)) {
	    for(i=0; i<N/2; i++) {	/* invert the values		*/
		u = v[i]; v[i] = v[N-1-i]; v[N-1-i] = u;
	    }
	    u = v[0];
	    for(i=0; i<N; i++) {
		v[i] = u - v[i];
	    }
	    sv = N*u - sv;
	    ts = nv*u - ts;
	}
	/*
	 * set up vals as ndv value-frequency pairs
	 */
	ndv = 1;
	for(i=1; i<N; i++) {
	    if (v[i] != v[i-1]) ndv++;
	}
	avf = (AVF *) ealloc(ndv * sizeof(AVF));
	u = 0;
	avf[0].val = v[0];
	avf[0].freq = 1;
	for(i=1; i<N; i++) {
	    if (v[i] != avf[u].val) {
		u++;
		avf[u].val = v[i];
		avf[u].freq = 1;
	    }
	    else {
		avf[u].freq++;
	    }
	}
	avf[0].cf = avf[0].freq;
	for(u=1; u<ndv; u++) {
	    avf[u].cf = avf[u].freq + avf[u-1].cf;
	}
	/*
	 * get cumulative sums
	 */
	for(i=1; i<nv; i++) {
	    v[i] += v[i-1];
	}
	cs = v = (unsigned *) erealloc(v, nv * sizeof(unsigned));
	/*
	 * Do it!
	 */
	if (dflag) {
	    dprint("nv = %d, N = %d, ndv = %d, ts = %d\n", nv, N, ndv, ts);
	    for(i=0; i<nv; i++) {
		dprint("cs[%d] = %d\n", i, cs[i]);
	    }
	    for(i=0; i<ndv; i++) {
		dprint("avf[%d].val = %d, .freq = %d, .cf = %d\n",
		    i, avf[i].val, avf[i].freq, avf[i].cf
		);
	    }
	}
	lecount(nv, ndv, ts, &lsum, &esum);
    }
    if ((nsw && inv) || (!nsw && !inv)) {
	gsum = pick(N, nv) - esum - lsum;
    }
    else {
	gsum = lsum;
	lsum = pick(N, nv) - esum - gsum;
    }
    if (cflag) printf("sum = e^%.15g\n", csum);
    free(avf);
    free(cs);
}

/* lecount -- count picks that sum to less than or equal a target
 * 
 * unsigned	cn;
 * unsigned	cndv;
 * unsigned	cts;
 * long		l, e;
 * lecount(cn, cndv, cts, &l, &e);
 *
 * lecount counts the number of picks of cn values that sum to < cts,
 * and that sum = cts, placing the results in l and e respectively.
 * The values are chosen from the first cndv value-frequency pairs in
 * avf.  lecount assumes the problem has been properly massaged by
 * legpicks beforehand.
 */
static void
lecount(cn, cndv, cts, lp, ep)
unsigned	cn;
unsigned	cndv;
unsigned	cts;
long		*lp;
long		*ep;
{
    register int	i;
    register int	v, f;
    int			ccn, ccts, ccndv;
    long		l, e;
    long		ls, es;
    long		pfi;		/* pick(f, i)			*/
    double		icsum = csum;	/* initial csum			*/

    indent++;
    *lp = *ep = 0L;
    /*
     * make sure inputs valid
     */
    if (cndv <= 0) error("Can't Happen");
    /*
     * prune some easily handled cases
     */
    if (
	(cts < cs[cn-1]) ||		/* any sum of cn is > cts	*/
	(cn > avf[cndv-1].cf)		/* there aren't cn vals		*/
    ) {
	dprint("Can't be done\n");
	goto leave;
    }
    /*
     * none to pick: sum is zero, and there's one way to pick nothing.
     */
    if (0 == cn) {
	if (cts == 0) *ep = 1L;		/* 1 works...			*/
	else *lp = 1L;			/* ...trust me			*/
	if (cflag) chkpick();
	goto leave;
    }
    /*
     * Only one value available.  Because cts < cs[cn-1], we know we
     * can pick cn of this value and be <= cts.
     */
    if (1 == cndv) {			/* only one val left		*/
	pfi = pick(avf[0].freq, cn);	/* # ways to pick cn from freq	*/
	if (cs[cn-1] == cts) *ep = pfi;	/* exactly equal		*/
	else *lp = pfi;			/* less				*/
	if (cflag) {
	    avf[0].n = cn;
	    chkpick();
	    avf[0].n = 0;
	}
	goto leave;
    }
    /*
     * Only one to be picked.  Find the # of vals < cts.
     */
    if (1 == cn) {
	for(i=0; i<cndv; i++) {
	    if ((v = avf[i].val) >= cts) break;
	    if (cflag) {
		avf[i].n = 1;
		chkpick();
		avf[i].n = 0;
	    }
	}
	*lp = (i > 0) ? avf[i-1].cf : 0L;
	if (v == cts) {			/* found an equal value		*/
	    *ep = avf[i].freq;		/* that many are equal		*/
	}
	goto leave;
    }
    /*
     * Both cn and cndv are >1.  Try to pick values from one of the
     * avf[i], and call lecount recursively to find compatible lower
     * picks.
     */
    ls = es = 0L;			/* accumulators for sums	*/
    for(ccndv=cndv; ccndv-- > 1;) {
	f = avf[ccndv].freq;		/* # we can pick		*/
	v = avf[ccndv].val;		/* what each one costs		*/
	ccn = cn;			/* current number to pick	*/
	ccts = cts;			/* current target sum		*/
	pfi = 1L;			/* pick(f, i)			*/
	for(i=1; i<=f; i++) {		/* i = # of this val to use	*/
	    if (ccts < v) break;	/* we can't afford any more	*/
	    ccn--;			/* OK: Buy one of these		*/
	    ccts -= v;			/* ...and pay for it		*/
	    pfi *= f+1 - i;		/* how many ways we can do it	*/
	    pfi /= i;
	    if (cflag) avf[ccndv].n = i;
					/* now find lower picks that fit*/
	    if (0 == ccn) {		/* (lecount doesn't like 0's)	*/
		if (ccts == 0) es += pfi;
		else ls += pfi;
		if (cflag) chkpick();
    	        if (cflag) chkdiff(cndv, icsum, csum, logint(ls), logint(es));
		break;			/* no more can be taken		*/
	    }
	    lecount(ccn, ccndv, ccts, &l, &e);
	    if (dflag) dprint("%ld * %ld = %ld < %d + %d * %d = %d\n",
		l, pfi, pfi * l, ccts, v, i, cts
	    );
	    if (dflag) dprint("%ld * %ld = %ld == %d + %d * %d = %d\n",
		e, pfi, pfi * e, ccts, v, i, cts
	    );
	    ls += pfi * l;		/* sum up			*/
	    es += pfi * e;
    	    if (cflag) chkdiff(cndv, icsum, csum, logint(ls), logint(es));
	    dprint("Total %ld <, %ld == %d\n",
		ls, es, cts
	    );
	}
	if (cflag) avf[ccndv].n = 0;
    }
    lecount(cn, 1, cts, &l, &e);	/* handle cndv == 1		*/
    *lp = ls += l;
    *ep = es += e;

leave:
    if (cflag) chkdiff(cndv, icsum, csum, logint(*lp), logint(*ep));
    if (dflag) dprint("%ld <, %ld == %d, %d vals chosen from 1st %d \n",
	*lp, *ep, cts, cn, cndv
    );
    indent--;
}

double
llpicks(an, aN, asum, avals)
int		an;
int		aN;
int		asum;
int		*avals;
{
    if (!argssame(an, aN, asum, avals)) {
	llegpicks();
    }
    return(llsum);
}

double
lepicks(an, aN, asum, avals)
int		an;
int		aN;
int		asum;
int		*avals;
{
    if (!argssame(an, aN, asum, avals)) {
	llegpicks();
    }
    return(lesum);
}

double
lgpicks(an, aN, asum, avals)
int		an;
int		aN;
int		asum;
int		*avals;
{
    if (!argssame(an, aN, asum, avals)) {
	llegpicks();
    }
    return(lgsum);
}

/* llegpicks -- legpicks, but use logs 					*/
static void
llegpicks()
{
    register int	i;
    register unsigned	u;
    register unsigned	*v;		/* scratch value array		*/
    int			adj;		/* value adjustment		*/
    unsigned		sv = 0;		/* summed values		*/
    int			nsw;		/* n has been switched		*/
    int			inv;		/* values have been inverted	*/

    indent = 0;
    dprint("n = %d, N = %d, sum = %d\n", n, N, sum);
    if (dflag) {
	for(i=0; i<N; i++) {
	    dprint("vals[%d] = %d\n", i, vals[i]);
	}
    }
    if (cflag) csum = LOG0;
    v = (unsigned *) ealloc(N * sizeof(unsigned));
    /*
     * adjust values to all >= 0
     */
    adj = -vals[0];
    for(i=0; i<N; i++) {
	v[i] = vals[i] + adj;
	sv += v[i];
    }
    if (sum + n*adj < 0) {
	llsum = lesum = LOG0;
    }
    else {
	ts = sum + n*adj;		/* target sum			*/
	/*
	 * use the smaller of n and N-n (becomes nv)
	 */
	if (nsw = (N-n < n)) {
	    nv = N - n;
	    ts = sv - ts;
	}
	else {
	    nv = n;
	}
	/*
	 * invert the values if mostly smaller than ts/nv
	 */
	if (inv = ((nv * (v[N/2] + v[(N-1)/2])) < 2*ts)) {
	    for(i=0; i<N/2; i++) {	/* invert the values		*/
		u = v[i]; v[i] = v[N-1-i]; v[N-1-i] = u;
	    }
	    u = v[0];
	    for(i=0; i<N; i++) {
		v[i] = u - v[i];
	    }
	    sv = N*u - sv;
	    ts = nv*u - ts;
	}
	/*
	 * set up vals as ndv value-frequency pairs
	 */
	ndv = 1;
	for(i=1; i<N; i++) {
	    if (v[i] != v[i-1]) ndv++;
	}
	avf = (AVF *) ealloc(ndv * sizeof(AVF));
	u = 0;
	avf[0].val = v[0];
	avf[0].freq = 1;
	for(i=1; i<N; i++) {
	    if (v[i] != avf[u].val) {
		u++;
		avf[u].val = v[i];
		avf[u].freq = 1;
	    }
	    else {
		avf[u].freq++;
	    }
	}
	avf[0].cf = avf[0].freq;
	for(u=1; u<ndv; u++) {
	    avf[u].cf = avf[u].freq + avf[u-1].cf;
	}
	/*
	 * get cumulative sums
	 */
	for(i=1; i<nv; i++) {
	    v[i] += v[i-1];
	}
	cs = v = (unsigned *) erealloc(v, nv * sizeof(unsigned));
	/*
	 * Do it!
	 */
	if (dflag) {
	    dprint("nv = %d, N = %d, ndv = %d, ts = %d\n", nv, N, ndv, ts);
	    for(i=0; i<nv; i++) {
		dprint("cs[%d] = %d\n", i, cs[i]);
	    }
	    for(i=0; i<ndv; i++) {
		dprint("avf[%d].val = %d, .freq = %d, .cf = %d\n",
		    i, avf[i].val, avf[i].freq, avf[i].cf
		);
	    }
	}
	llecount(nv, ndv, ts, &llsum, &lesum);
    }
    if ((nsw && inv) || (!nsw && !inv)) {
	lgsum = logdif(logdif(logpick(N, nv), lesum), llsum);
    }
    else {
	lgsum = llsum;
	llsum = logdif(logdif(logpick(N, nv), lesum), lgsum);
    }
    if (cflag) printf("sum = e^%g\n", csum);
    free(avf);
    free(cs);
}

/* llecount -- count picks that sum to less than or equal a target
 * 
 * unsigned	cn;
 * unsigned	cndv;
 * unsigned	cts;
 * double	ll, le;
 * llecount(cn, cndv, cts, &ll, &le);
 *
 * llecount is just like lecount, except that it returns in ll and le
 * the natural logs of the desired numbers.  llecount is less efficient
 * than lecount because it does floating point arithmetic and evaluates
 * transcental functions, but it uses lecount whenever numbers drop low
 * enough, so it is hoped that it will not be disastrously slow.
 */
static void
llecount(cn, cndv, cts, llp, lep)
unsigned	cn;
unsigned	cndv;
unsigned	cts;
double		*llp;
double		*lep;
{
    register int	i;
    register int	v, f;
    int			ccn, ccts, ccndv;
    long		l, e;
    double		ll, le;
    double		lls, les;
    double		lpfi;		/* logpick(f, i)		*/
    double		icsum = csum;	/* initial csum			*/

    indent++;
    *llp = *lep = LOG0;
    /*
     * make sure inputs valid
     */
    if (cndv <= 0) error("Can't Happen");
    /*
     * prune some easily handled cases
     */
    if (
	(cts < cs[cn-1]) ||		/* any sum of cn is > cts	*/
	(cn > avf[cndv-1].cf)		/* there aren't cn vals		*/
    ) {
	dprint("Can't be done\n");
	goto leave;
    }
    /*
     * none to pick: sum is zero, and there's one way to pick nothing.
     */
    if (0 == cn) {
	if (cts == 0) *lep = LOG1;	/* 1 works...			*/
	else *llp = LOG1;		/* ...trust me			*/
	if (cflag) chkpick();
	goto leave;
    }
    /*
     * Only one value available.  Because cts < cs[cn-1], we know we
     * can pick cn of this value and be <= cts.
     */
    if (1 == cndv) {			/* only one val left		*/
					/* # ways to pick cn from freq	*/
	lpfi = logpick(avf[0].freq, cn);
					/* exactly equal		*/
	if (cs[cn-1] == cts) *lep = lpfi;
	else *llp = lpfi;		/* less				*/
	if (cflag) {
	    avf[0].n = cn;
	    chkpick();
	    avf[0].n = 0;
	}
	goto leave;
    }
    /*
     * Only one to be picked.  Find the # of vals < cts.
     */
    if (1 == cn) {
	for(i=0; i<cndv; i++) {
	    if ((v = avf[i].val) >= cts) break;
	    if (cflag) {
		avf[i].n = 1;
		chkpick();
		avf[i].n = 0;
	    }
	}
	*llp = (i > 0) ? logint(avf[i-1].cf) : LOG0;
	if (v == cts) {			/* found an equal value		*/
	    *lep = logint(avf[i].freq);	/* that many are equal		*/
					/* ...and all before are less	*/
	}
	goto leave;
    }
    /*
     * Both cn and cndv are >1.  Try to pick values from one of the
     * avf[i], and call lecount recursively to find compatible lower
     * picks.
     */
    lls = les = LOG0;			/* accumulators for sums	*/
    for(ccndv=cndv; ccndv-- > 1;) {
	f = avf[ccndv].freq;		/* # we can pick		*/
	v = avf[ccndv].val;		/* what each one costs		*/
	ccn = cn;			/* current number to pick	*/
	ccts = cts;			/* current target sum		*/
	lpfi = LOG1;			/* pick(f, i)			*/
	for(i=1; i<=f; i++) {		/* i = # of this val to use	*/
	    if (ccts < v) break;	/* we can't afford any more	*/
	    ccn--;			/* OK: Buy one of these		*/
	    ccts -= v;			/* ...and pay for it		*/
					/* how many ways we can do it	*/
	    lpfi += logint(f+1 - i) - logint(i);
	    if (lpfi < LOG1) lpfi = LOG1;
					/* now find lower picks that fit*/
	    if (cflag) avf[ccndv].n = i;
	    if (0 == ccn) {		/* (lecount doesn't like 0's)	*/
		if (ccts == 0) les = logsum(les, lpfi);
		else lls = logsum(lls, lpfi);
		if (cflag) chkpick();
    	        if (cflag) chkdiff(cndv, icsum, csum, lls, les);
		break;			/* no more can be taken		*/
	    }
	    if (logpick(avf[ccndv].cf, ccn) < LOGMAXL) {
					/* safe to switch to integers	*/
	        lecount(ccn, ccndv, ccts, &l, &e);
		ll = logint(l);
		le = logint(e);
	    }
	    else {			/* too big: have to use logs	*/
		llecount(ccn, ccndv, ccts, &ll, &le);
	    }
	    if (dflag) dprint("e^%g * e^%g = e^%g < %d + %d * %d = %d\n",
		ll, lpfi, logmul(lpfi, ll), ccts, v, i, cts
	    );
	    if (dflag) dprint("e^%g * e^%g = e^%g == %d + %d * %d = %d\n",
		le, lpfi, logmul(lpfi, le), ccts, v, i, cts
	    );
	    lls = logsum(lls, logmul(lpfi, ll));
	    les = logsum(les, logmul(lpfi, le));
    	    if (cflag) chkdiff(cndv, icsum, csum, lls, les);
	}
	if (cflag) avf[ccndv].n = 0;
    }
    lecount(cn, 1, cts, &l, &e);	/* handle cndv == 1		*/
    *llp = lls = logsum(lls, logint(l));
    *lep = les = logsum(les, logint(e));

leave:
    if (cflag) chkdiff(cndv, icsum, csum, *llp, *lep);
    if (dflag) dprint("e^%g <, e^%g == %d, %d vals chosen from 1st %d \n",
	*llp, *lep, cts, cn, cndv
    );
    indent--;
}

/* chkpick -- check that a pick actually is valid			*/
static void
chkpick()
{
    register int	i;
    register int	cn = 0, s = 0;
    double		np = LOG1;

    for(i=0; i<ndv; i++) {
	if (i != 0) printf(" + ");
	cn += avf[i].n;
	s += avf[i].val * avf[i].n;
	np = logmul(np, logint(pick(avf[i].freq, avf[i].n)));
	printf("%d*%d", avf[i].n, avf[i].val);
    }
    csum = logsum(csum, np);
    if (s == ts) {
	printf(" = %d = %d\n", s, ts);
    }
    else {
	printf(" = %d < %d\n", s, ts);
    }
    if (cn != n) {
	dprint("n not correct\n");
    }
    if (s > ts) {
	dprint("sum too big\n");
    }
}

#define	DIFTOL	(1e-10)
/* chkdiff -- check that sums increased as they should			*/
static void
chkdiff(cndv, icsum, csum, ll, le)
int	cndv;
double	icsum, csum, ll, le;
{
    register int	i;
    double		np = LOG1;
    double		new;

    for(i=cndv; i<ndv; i++) {
	np = logmul(np, logint(pick(avf[i].freq, avf[i].n)));
    }
    new = logmul(np, logsum(ll, le));
    if (DIFTOL < fabs(csum - logsum(icsum, new))) {
	printf("%.12g + %.12g <> %.12g\n",
	    exp(icsum), exp(new), exp(csum)
	);
    }
}