[RESOLVED] Binomial coefficiemt question

[bilm_ks]

Is there optimum method for the calculation C(n, r) = n!/(r!*(n-r)!)
e.g. since this is an integer can i avoid the division?

[zaza]

since this is an integer can i avoid the division?

Why would it being an integer mean you can avoid the division? nCr is what it is - there isn't a shortcut unless you want to program the coefficients in by hand or use Pascal's triangle to work it out. The former might be marginally quicker but very limited. The latter would be very lengthy, particularly for higher orders.

zaza

[jemidiah]

I'm not certain about this, but like most of the other fuctions that calculators calculate (sin, cos, sqrt, etc.), I would assume there is something that would make it easier to calculate (for a machine) than using the definition you've given (maybe some type of series, or even odd binary tricks may help). Sadly, I have no idea what that method is, sorry

[krtxmrtz]

Is there optimum method for the calculation C(n, r) = n!/(r!*(n-r)!)
e.g. since this is an integer can i avoid the division?

Well, unless the involved numbers are large I don't see why division should be a problem. Maybe you can decompose n and r into their prime factors and see if the division(s) can be remnoved altogether or at least simplified.
If n is large you should use Stirling's formula, which in its simplest form has no divisions... but more complicated stuff.

[bilm_ks]

Thanks all for your posts.

What i mean is {after the obvious n!/((n-r)!*r!) = ((n-r+1)*(n-r+2)*...*n)/r!} since this is an integer is there any way to remove more factors? (For zaza: C(n,r) = k no denominator). For krtxmrtz: n is in most cases small and n!~SQR(2*PI*n)*n^n*exp(-n) is an aproximate type.
Thanks again.

[zaza]

Yes, thank you for pointing out what an integer is. But just because you know that it will be an integer doesn't mean that you can skip the division. Is your problem that you want to evaluate the factorials and so you have to do 1x2x3x..., so you want to minimise this?

One possibility to cut this process down would be to evaluate n-r against r to see which is bigger. Then, since n > r, n! / r!(n-r!) could become:

If r > (n-r): (n!/r!) x (1/(n-r)!)

or vice versa if (n-r) > r.
Then you can evaluate the product from r to n in the first division, because the entire r! cancels.
An example:

13 C 6 = 13! / 7! 6!

-> (13! / 7!) * (1/6!)

-> (8 x 9 x 10 x 11 x 12 x 13) * (1/ (1 x 2 x 3 x 4 x 5 x 6))

I think that's about as simple as you're likely to get.
Mathematically, of course, the way you'd describe (8 x 9 x ...) algebraically is as I put in the previous line - and there's no getting around the fact that that is the mathematical way to describe it.

But in terms of evaluating the factorials by multiplication, it's easier than having to do all 3 of them separately.

Of course, another way would be to do a For Next loop from 1 to n, with a check to see if the counter matches n, r or n-r. If so, then copy the product thus far out to a variable. That way you only have to do 1 loop through multiplying.


zaza

[bilm_ks]

One possibility to cut this process down would be to evaluate n-r against r to see which is bigger.

zaza

If r > n-r -> n-1 muls and 1 div
if r < n-r -> n-1 muls and 1 div also

[krtxmrtz]

I'm not sure what your concern is exactly, efficiency in terms of speed, lines of code... I'd just like to draw your attention to the fact that the calculation of factorials can be ellegantly accomplished by a recursive subroutine. I posted something about this not too long ago, let me know if you're interested.

[MagicT]

Just for grins, I computed the combination formula using the most inefficient code I could muster:

Private Sub cmdGo_Click()
   Dim i As Long
   Dim w As New StopWatch
   Dim r, n, x As Integer
   Dim rsum, nsum, xsum As Double
   Dim comb As Double
   
   r = 50
   n = 150
   x = n - r
   w.StartClock
   rsum = 1
   xsum = 1
   nsum = 1
 
   For i = 1 To n
       nsum = nsum * i
   Next i
   For i = 1 To r
       rsum = rsum * i
    Next i
    For i = 1 To x
      xsum = xsum * i
    Next i
    comb = nsum / (rsum * xsum)
    Label2.Caption = comb
    w.StopClock Label1
   
End Sub

Stopwatch is a class I wrote that uses the GetTickCount API to count how many ticks occur between start and stop...I get zero ticks for any combination problem I run. (Note that these sums get bigger than doubles with relatively low values of r and n)
I'm running on an old 1500 mhz cpu.

So, speed is not an issue! Go for esthetics...

-MagicT

[sql_lall]

Is there optimum method for the calculation C(n, r) = n!/(r!*(n-r)!)
e.g. since this is an integer can i avoid the division?

Integer division, in this case, isn't a problem if you structure it right.
Say you want to find 9 C 4 (or, 9 C 5 even)
this is equal to: (9 * 8 * 7 * 6) / (4 * 3 * 2 * 1)
= ((((((((6)/1)*7)/2)*8)/3)*9)/4)
if you do it in that order, you should *always* be dealing with whole integers.
because you're basically working out 6C1, multipying by (7/2) to get 7C2 (also an integer), multiplying by (8/3) to get 8C3, then finally by (9/4) to get 9C4.
At each stage, so long as you multiply first, you'll always be dealing with integers.

[MarcoA]

Hello:

You can use the recurrence relation:

C(n,r) = C(n-1,r) + C(n-1,r-1)

Then we have the binomials coefficients with the Pascal triangle:

1
1 1..................................n=1
1 2 1...............................n=2
1 3 3 1............................n=3
1 4 6 4 1..........................n=4
1 5 10 10 5 1....................n=5
..................

You have no divisions!

Regards

[zaza]

If r > n-r -> n-1 muls and 1 div
if r < n-r -> n-1 muls and 1 div also

Oh my goodness, are you from Gibberland too? You speak fluent Gibberish...

[bilm_ks]

First of all i want to apologize for the foolish

If r > n-r -> n-1 muls and 1 div
if r < n-r -> n-1 muls and 1 div also

I ended up with two approaches according to your suggestions (and thank you all for this). The first one is based on zaza’s and sql_lall comments. The second uses the Pascal triangle as MarcoA suggested. I tested the following code

#include <windows.h>
#include <stdio.h>

// constants for testing
#define TEST_REP 0xfffff // 1 M loops
const int  TestN[] = {22,23,24,25,26,27,28};

// global array
int * g_pArr;

// using PascalTriangle
int PascalTriangle(int n, int r)
{
	// need an array with n+1 elements
	int * pArr = new int[n + 1];
	
	int res, prev, count;

	// init arr for 2
	pArr[0] = 1;
	pArr[1] = 2;
	pArr[2] = 1;
	
	// starting from 3 to n-1
	for (count = 3; count < n; count++)
	{
		// first = 1
		pArr[0] = prev = 1;

		for (int i = 1; i < count; i++)
		{
			// save pArr[i] for the next loop
			int tmp = pArr[i];
			pArr[i] += prev;     // C(count, i) = C(count-1, i)+C(count-1, i-1)
			prev = tmp;
		}

		// last = 1
		pArr[count] = 1;
	}

	// C(n,r)
	res = pArr[r] + pArr[r-1];

	// free mem
	delete [] pArr;

	return res;
}

// using PascalTriangle with a global array this time
int PascalTriangleExt(int n, int r)
{	
	int res, prev, count;

	// init arr for 2
	g_pArr[0] = 1;
	g_pArr[1] = 2;
	g_pArr[2] = 1;
	
	// starting from 3 going recursively to n-1
	for (count = 3; count < n; count++)
	{
		
		// first = 1
		g_pArr[0] = prev = 1;

		for (int i = 1; i < count; i++)
		{
			// save g_pArr[i] for the next loop
			int tmp = g_pArr[i];
			g_pArr[i] += prev; // C(count, i) = C(count-1, i)+C(count-1, i-1)
			prev = tmp;
		}

		// last = 1
		g_pArr[count] = 1;
	}

	// C(n,r)
	res = g_pArr[r] + g_pArr[r-1];

	return res;
}


// normal approach
int BinomialCoeff(int n, int r)
{
	int res, i, j;

	res = 1;
	
	if (r > (n - r))  
		// simplifying r!
		for (i = r + 1, j = 1 ; i <= n; i++, j++)
		{
			res *= i;
			res /= j;
		}
	else
		// simplifying (n-r)!
		for (i = n - r + 1, j = 1; i <= n; i++, j++)
		{
			res *= i;
			res /= j;
		}

	return res;
}

// test function
void Test(int (*func)(int,int), const char* name)
{
	unsigned int start, fin;

	start = GetTickCount();

	for (int i = 0; i < TEST_REP; i++)
		func(TestN[i%7], TestN[i%7]/2); // choose the hard case (r = n/2)

	fin = GetTickCount();

	printf("%s : %u ms\n\n", name, fin-start);
}

int main(void)
{
	const char* normal = "Normal";
	const char* pasc_tr = "Pascal triangle";
	const char* pasc_glob = "Pascal triangle (no allocation)";

	// mem alloc for global array
	g_pArr = new int[30];

	// testing each method
	Test(&PascalTriangle, pasc_tr);
	Test(&BinomialCoeff, normal);
	Test(&PascalTriangleExt, pasc_glob);

	return 0;
}

The results were much better for the normal approach (456 ms) than the Pascal triangle (1141 ms) or the Pascal triangle again but with no array allocation in the function body (750 ms)

Thanks all and especially zaza, sql_lall and marcoA