Conference rusure::math

    re .0
    
    C(n,r) represents the number of ways of selecting a set of r objects
    from a pool of n. Note that I use the word set here to emphasise that
    the order of the objects is unimportant, as is applicable to your
    problem.
    
    Hence if there are X voters and N of them voted for you then there are
    C(X,N) ways of selecting the set of voters that voted for you.
    
                n!
    C(n,r) = --------
             r!(n-r)!
    
    where k! = 1.2.3...k.
    
    If avoiding overflow is important then this is not the best way
    actually to calculate C(n,r). If n is really large so that an
    approximation is acceptable then there are some simple approximation
    formulae.
    
    I am not familiar with a good algorithm to generate all combinations.
    
    What sort of size for X and N are we talking?

The number of voters in this problem are typicaly less than 50, and the vote 
acceptance variable may be upto 48!  In this domain, there may be a 
computational overflow.  IS there any other way of simply solving this?

    re .2
    
    If you don't require an exact integer answer and for number of voters
    less than 50 then you could use G-float double precision to get a
    reasonable value for the number of combinations without having to be
    careful about overflow.
    
    Actually enumerating the combinations would be infeasible e.g. for 50
    voters and 25 supporters.

I assume that when your talking about computational overflow, you mean
intermediate overflow.  That is, in the computation of

		            n!
		C(n,k) = --------
		         k!(n-k)!

At least one of n!, k! or (n-k)! overflows, but the quotient doesn't.  If
that is the case, a reasonable accurate result can be obtained using the 
lgamma(x) funtion.  The gamma function is (for discussion purposes) a
generalization of factorials to non-integer values.  For integer values

		gamma(n+1) = n!

The lgamma function is the natural log of the gamma function.  lgamma(x) is
frequently used because gamma(x) is such a rapidly increasing function it 
overflows most finite representations and because usually the real quantity
of interest is a ratio of gamma functions, like C(n,k).  So one might compute
C(n,k) by

		C(n,k) = exp(log(n!/(k!(n-k)!)))
		       = exp(log(n!) - log(k!) - log((n-k)!))
		       = exp(lgamma(n+1) - lgamma(k+1) - lgamma(n-k))

There are very good implementations of the lgamma function on all of the DEC
Alpha platforms.

re .2

i have not done things like this for a while, so this is not promised 
to work, but ask lisp to evaluate it.

or use one of the extended precision packages mentioned in this conference.

Well, consider X=n!/(r!(n-r)!).  If n=50, the largest value for x is
probably near r=25.  For that case we have

	X=50*49*48...*25*24*23...*2 / (25*24*23...*2 * 25*24*23...*2)


The 25! cancels on top and bottom, giving

	X=50*49*48...*26 / 25*24*23...*2

The 50 cancels with the 25*2:

	X=49*48...*26 / 24*23...*3

27 cancels with 3 and 9.
28 cancels with 4 and 8.
30 cancels with 6 and 5.

	X=49*48*47*46*45*44*43*42*41*40*39*38*37*36*35*34*33*32*31*29*26 /
		24*23*22*21*20*19*18*17*16*15*14*13*12*11*10*7

Useful reductions now are 48/24, 46/23, 44/22, 42/21, 40/20, 38/19, 36/18,
	34/17, 32/16 :

	X=49*2*47*2*45*2*43*2*41*2*39*2*37*2*35*2*33*2*31*29*26 /
		15*14*13*12*11*10*7

Now we apply 39/13, 35/7, 33/11:

	X=49*2*47*2*45*2*43*2*41*2*3*2*37*2*5*2*3*2*31*29*26 /
		15*14*12*10

Now 3*5/15, 3*2*2/12:

	X=49*2*47*2*45*2*43*2*41*2*2*37*2*31*29*26 /
		14*10

Take one of the 7's from the 49 combined with a 2, and cancel with the 14,
then take a 5 from the 45 with another 2, to cancel with the 10:

	X=7*47*9*2*43*2*41*2*2*37*2*31*29*26

Combining the 2's:

	X=7*47*9*43*41*37*32*31*29*26


Not quite small enough for a vax (or even my alpha/vms dcl [why not???] )

/Eric

re .6
    
>Well, consider X=n!/(r!(n-r)!).  If n=50, the largest value for x is
>probably near r=25.  For that case we have
    
    For a given even n C(n,r) has a maximum at r=n/2. (Think Pascal's triangle).
    
    My calculator gives C(50,25) as around 1.3E14.
    
>Not quite small enough for a vax (or even my alpha/vms dcl [why not???] )
    
    Hence as you observe too large for 32-bit integer arithmetic. Alpha
    VMS is initially a straight port of VAX VMS hence DCL can't do 64-bit
    arithmetic even on Alpha. Maybe later... Followups in VMSNOTES. (-:

For small N, it's not too hard to take advantage of the cancellation idea
in a program:

$ type comb.for
        integer*8 n,k,i,m
        accept *,n,k
        if (k.gt.n-k) k=n-k
        m = 1
        do i = 1,k
          m = m*(n-i+1)/i
        enddo
        print *,m
        end
$ fort/check=overflow comb
$ link comb
$ run comb
50,25
      126410606437752
$ run comb
60,30
   118264581564861424
$ run comb
70,21
   385439532530137800
$ run comb
70,22
%SYSTEM-F-HPARITH, high performance arithmetic trap, Imask=00010000, Fmask=00000000, summary=40, PC=8048C798, PS=0000001B
-SYSTEM-F-INTOVF, arithmetic trap, integer overflow at PC=8048C798, PS=0000001B
%TRACE-F-NOMSG, Message number 0009804C
 Image Name   Module Name     Routine Name    Line Number  rel PC      abs PC
                                                        0 8048C798    8048C798
                                                        0 82FBDA50    82FBDA50

That's quite a compression of all the cancelation stuff I discussed !  I'm
still trying to figure out how your program works.  Thanks.

/Eric

The heart of the program (there's always one line that's the heart of
any computer program right?) is:

     integer*8 n,k,i,m
        accept *,n,k
        if (k.gt.n-k) k=n-k
        m = 1
        do i = 1,k
==>          m = m*(n-i+1)/i
        enddo
        print *,m
        end



Let's suppose we're doing (50,25).  The generated value is:

	1*50/1*49/2*48/3*47/4...*26/25

This value is computed in order from left to right, so we need to convince
ourselves (well, I need to convince myself!) that none of those divisions
will ever be non-integer.

Certainly the /2 is o.k. since either n or n-1 is even.

If I were being sloppy, I could convince myself that /3 is o.k. since n or
n-1 or n-2 is divisible by 3.  But I'm not quite that sloppy, so I need
to worry that perhaps whichever of those 3 was divisible by 3 no longer is,
once we did the /2.  But on further thought, I realize that the /2 can't
"damage" the divisible-by-3-ness of a number.

But what about the /4 ?  Certainly either n or n-1 or n-2 or n-3 is divisible
by 4, but suppose it was the n or n-1 that was divisible by 4.  One of those
was already divided by 2, so can we be certain that what's left is divisible
by 4 ?

Then of course there's the whole question of why the entire computation
just happens to be equivalent to (n,k), even if all the divisibility concerns
are satisfied.

I guess this slow mathematician here will need to think about all of this
a bit more while the whole thing must be obvious to the rest of you...

Thanks.

/Eric

o.k. I'm not yet convinced that all those immediate operations will be
exactly integers, but I think I've figured out why if they are integers
that the expression is equivalent to (n,k).  We have

	50*49*48...*1
	--------------------------------
	25*24*23...*1    * 25*24*23...*1

The first set cancels, leaving us with

	50*49*48...*26
	--------------
	25*24*23...*1

which we can reverse the bottom of:

	50*49*48...*26
	--------------
	1*2*3...*25

which can of course be written as

	1*50/1*49/2*48/3...*26/25

but this doesn't answer the first question yet about guaranteeing that
each intermediate result is an exact integer...

/Eric

> But what about the /4 ?  Certainly either n or n-1 or n-2 or n-3 is divisible
> by 4, but suppose it was the n or n-1 that was divisible by 4.  One of those
> was already divided by 2, so can we be certain that what's left is divisible
> by 4 ?
    
    Because no matter which one was divisible by 4 you still have another
    one divisible by 2.  E.g., if n-1 is divisible by 4 then n-3 is
    divisible by 2 (but not 4).
    
    Similarly when trying to divide by 6: n, n-1, n-2, n-3, n-4, n-5;
    Whichever one is divisible by six (say it's n-4) leaves you another
    that's divisible by 3 (in this case n-1), one that's divisible by 4
    (either n-2 or n) and one that's divisible by 2 (the "4" candidate that
    wasn't divisible by 4). And so on.
    
      John

    Re .9, .10, .11:
    
    At step i, the value of m is (n,i).  The number of combinations is
    known to be an integer (particularly because it can be defined
    inductively by adding previous values, which are also integers).
    
    
    				-- edp
    
    
Public key fingerprint:  8e ad 63 61 ba 0c 26 86  32 0a 7d 28 db e7 6f 75
To find PGP, read note 2688.4 in Humane::IBMPC_Shareware.

re .14
    
>   depending on whether rem(a,c) + rem(b,c) < or > 1.
    
    depending on whether rem(a,c) + rem(b,c) < or >= 1.
    
    This is also a somewhat non-obvious use of the name "rem". Given the
    names floor() and ceil() perhaps the right function name is height(). (-:
    
    Equivalently expressed...
    
    depending on whether a mod c + b mod c < or >= c.