Developers Club geek daily blog

1 year, 1 month ago
At a solution of problems of combination theory often there is a need for calculation binomial coefficients. Newton binomial, i.e. decompositionimage also uses binomial coefficients. For their calculation it is possible to use the formula expressing binomial coefficient through factorials:image or to use a recurrent formula:image From Newton binomial and a recurrent formula it is clear that binomial coefficients — integral numbers. On this example there was a wish to show that even at a solution of a simple task it is possible to step on a rake.
Before passing to writing of procedures of calculation, we will consider a so-called triangle of Pascal.
          1
       1     1
    1     2     1
  1    3     3     1
1   4     6     4     1

or it, but is a little in other type. In the left column of a line value n, is farther in a line of valueimage for k=0. n
 n          биномиальные коэффициенты
 0      1
 1      1      1
 2      1      2      1
 3      1      3      3      1
 4      1      4      6      4      1
 5      1      5     10     10      5      1
 6      1      6     15     20     15      6      1
 7      1      7     21     35     35     21      7      1
 8      1      8     28     56     70     56     28      8      1
 9      1      9     36     84    126    126     84     36      9      1
10      1     10     45    120    210    252    210    120     45     10      1
11      1     11     55    165    330    462    462    330    165     55     11      1
12      1     12     66    220    495    792    924    792    495    220     66     12      1
13      1     13     78    286    715   1287   1716   1716   1287    715    286     78     13      1
14      1     14     91    364   1001   2002   3003   3432   3003   2002   1001    364     91     14      1
15      1     15    105    455   1365   3003   5005   6435   6435   5005   3003   1365    455    105     15      1
16      1     16    120    560   1820   4368   8008  11440  12870  11440   8008   4368   1820    560    120     16      1


In full accordance with a recurrent formula of valueimage are equal 1 and any number is equal to the sum of the number costing over it and numbers "over it + a step to the left". For example, in 7y to a line number 21, and in 6y to a line numbers 15 and 6: 21=15+6. It is visible also that values in line are symmetric concerning the middle of a line, i.eimage. This property of symmetry of Newton binomial is relative an and b and it is visible in a factorial formula.
I imagewill also be lower for binomial coefficients to use representation of C (n, k) (it it is simpler to gather and the formula picture can be inserted not everywhere.

Calculation of binomial coefficients through a factorial formula


as binomial coefficients non-negative, we will use unsigned type in calculations.
// расчет факториала n
unsigned fakt(int n)
{
   for (unsigned r=1;n>0;--n)  
         r*=n;
   return r;
}
// расчет C(n,k)
unsinged bci(int n,int k)
{
   return fakt(n)/(fakt(k)*fakt(n-k));
}


Let's cause the bci (10,4) function — it will return 210 and this correct value of coefficient C (10,4). Means, the problem of calculation is solved? Yes, it is solved. But not absolutely. We did not answer a question: at what maximum values n, k the bci function will correctly work? Before beginning to look for the answer, we will agree that the unsigned int type of the 4th used by us byte and the maximum value is equal to 232-1=4'294'967'295. At what n, k C (n, k) will exceed it? Let's address Pascal's triangle. It is visible that the maximum values are reached in the middle of the line, i.e. at k=n/2. If n chetno, then is available one maximum value and if n nechetno, then their two. Exact value C (34,17) is equal 2333606220, and exact value C (35,17) is equal 4537567650, i.e. there is already more maximum unsigned int.
Let's write test procedure
void test()
{
    for (n=10;n<=35;++n) 
    printf("%u %u",n,bci(n,n/2);
 // для C++ можно использовать cout<<n<<" "<<bci(n,n/2)<<endl;
}

Let's start it and we will see
10 252
11 462
12 924
13 532
14 50
15 9
16 1
17 2
18 1
19 0
20 0
21 1
22 0
23 4
24 1
25 0
26 1
27 0
28 1
29 0
30 0
31 0
32 2
33 2
34 0
35 0


Value in the next line has to be approximately twice more, than in previous. Therefore the last correctly calculated coefficient (the cm a triangle is higher) is the C (12,6) Though unsigned int contains 4 billion, values less than 1000 are correctly calculated. Here those time, why so? It is all about our bci procedure, more precisely in line which at first calculates a large number in numerator, and divides it then into a large number in a denominator. For calculation of C (13,6) at first 13 are found!, and this number> 6 billion and it does not get into unsigned int.
How to optimize calculationimage? Very simply: let's reveal 13! also we will reduce numerator and a denominator on 7! As a result it will turn outimage. Let's program calculation for this formula
unsigned bci(int n,int k) 
{
	if (k>n/2) k=n-k; // возьмем минимальное из k, n-k.. В силу симметричность C(n,k)=C(n,n-k)
	if (k==1)  return n;
	if (k==0)  return 1;
	unsigned r=1;
	for (int i=1; i<=k;++i) {
		r*=n-k+i;
		r/=i;
	}
	return r;
}

:
and again we will test
10 252
11 462
12 924
13 1716
14 3432
15 6435
16 12870
17 24310
18 48620
19 92378
20 184756
21 352716
22 705432
23 1352078
24 2704156
25 5200300
26 10400600
27 20058300
28 40116600
29 77558760
30 155117520
31 14209041
32 28418082
33 39374192
34 78748384
35 79433695



Explicitly better, the error arose when calculating C (31,15). The reason is clear — the same overflow. At first we multiply on 31 (an op-pas — overflow, then we divide on 15). And what if to use a recursive formula? There only addition, overflow should not be.
Well, we try:
unsigned bcr(int n,int k) 
{
	if (k>n/2) k=n-k;
	if (k==1)  return n;
	if (k==0)  return 1;
	return bcr(n-1,k)+bcr(n-1,k-1);
}
void test()
{
    for (n=10;n<=35;++n) 
    printf("%u %u",n,bcr(n,n/2);
 // для C++ можно использовать cout<<n<<" "<<bcr(n,n/2)<<endl;
}


Test result
10 252
11 462
12 924
13 1716
14 3432
15 6435
16 12870
17 24310
18 48620
19 92378
20 184756
21 352716
22 705432
23 1352078
24 2704156
25 5200300
26 10400600
27 20058300
28 40116600
29 77558760
30 155117520
31 300540195
32 601080390
33 1166803110
34 2333606220
35 242600354


Everything that got into unsigned int, was counted correctly. Here only the line reckoned with n=34 about a minute. When calculating C (n, n/2) two recursive calls therefore calculation time exponential depends on n become. What to do — it turns out or it is inexact, or slowly. An output — in use of 64 bit variables.

The note by results of discussions: at the end of article the section where the simple and fast option "bcr with storing" of one of participants of discussion is given is added.

Use of 64 bit types for calculation of C (n, k)


Let's replace as bci unsigned int with unsigned long long and we will test in n=34 range. 68. n=34 is the maximum value which correctly is considered bci, and C (68,34) ~ 2.8*1019 does not get into unsigned long long ~ 1.84*1019 any more
unsigned long long bcl(int n,int k) 
{
	if (k>n/2) k=n-k;
	if (k==1)  return n;
	if (k==0)  return 1;
	unsigned long long r=1;
	for (int i=1; i<=k;++i) {
		r*=n-k+i;
		r/=i;
	}
	return r;
}

void test()
{
    for (n=34;n<=68;++n) 
    printf("%llu %llu",n,bcl(n,n/2));
 // для C++ можно использовать cout<<n<<" "<<bcl(n,n/2)<<endl;
}
 

Test result
34 2333606220
35 4537567650
36 9075135300
37 17672631900
38 35345263800
39 68923264410
40 137846528820
41 269128937220
42 538257874440
43 1052049481860
44 2104098963720
45 4116715363800
46 8233430727600
47 16123801841550
48 32247603683100
49 63205303218876
50 126410606437752
51 247959266474052
52 495918532948104
53 973469712824056
54 1946939425648112
55 3824345300380220
56 7648690600760440
57 15033633249770520
58 30067266499541040
59 59132290782430712
60 118264581564861424
61 232714176627630544
62 465428353255261088
63 321255810029051666
64 66050867754679844
65 454676336121653775
66 350360427585442349
67 23341572944240599
68 46683145888481198


We see that the error arises at n=63 for the same reason, as in bci. At first multiplication on 63 (and overflow), then division on 31.

Further increase of accuracy and calculation at n> 67


The error arose at n=63, and at n=68 the result does not get into unsigned64 any more. Therefore it is possible to tell "to n<=62 функция bcl считает правильно, дальнейшее увеличение точности требует либо int128 либо длинной арифметики». А если очень высокая точность не нужна, но хочется считать биномиальные коэффициенты при n=100...1000? Снова берем процедуру bci и меняем в ней типы unsigned int на double:
double bcd(int n,int k) 
{
	if (k>n/2) k=n-k; // возьмем минимальное из k, n-k.. В силу симметричности C(n,k)=C(n,n-k)
	if (k==1)  return n;
	if (k==0)  return 1;
	double r=1;
	for (int i=1; i<=k;++i) {
		r*=n-k+i;
		r/=i;
	}
	return ceil(r-0.2); // округлим до ближайшего целого, отбросив дробную часть
                       // комментарий изменен после обсуждений: такой способ использован,  чтобы расхождение с точным значением 
                       // было как можно позже. Где-то оно превышало 0.5 и простой round не годился 
}
void testd()
{
    for (n=50;n<=1000;n+=50) 
    printf("%d %.16e\n",n,bcd(n,n/2));
 // для C++ можно использовать cout<<n<<" "<<bcd(n,n/2)<<endl;
}


Test result
50 1.2641060643775200e+014
100 1.0089134454556417e+029
150 9.2826069736708704e+043
200 9.0548514656103225e+058
250 9.1208366928185793e+073
300 9.3759702772827310e+088
350 9.7744946171567713e+103
400 1.0295250013541435e+119
450 1.0929255500575370e+134
500 1.1674431578827770e+149
550 1.2533112137626624e+164
600 1.3510794199619429e+179
650 1.4615494992533863e+194
700 1.5857433585316801e+209
750 1.7248900341772600e+224
800 1.8804244186835327e+239
850 2.0539940413411323e+254
900 2.2474718820660189e+269
950 2.4629741379276902e+284
1000 2.7028824094543663e+299

Even for n=1000 it was succeeded to count! Overflow of double will happen at n=1030.
Discrepancy of calculation of bcd with exact value begins with n=57. It small — only 8. At n=67 a deviation 896.

For thrill-seekers and "Olympians"


In principle, for practical problems of accuracy of the bcd function it is enough, but in Olympiad tasks the tests "on the verge" often are given. I.e. theoretically can the task where the accuracy of double has not enough also of C will meet (n, k) gets into unsigned long long hardly. How to avoid overflow for such extreme cases? It is possible to use the recursive algorithm. But if he for n=34 considered minute, then for n=67 will consider years 100. It is possible to remember rasschtanny values (cm Addition after a publiukation). It is also possible to use a recursion not for all n and k but only for "enough big". Here calculation procedure which considers correctly for <=67. Иногда и для n>n67 at small k (for example, considers C (82,21) =1.83*1019).
unsigned long long bcl(int n,int k) 
{
	if (k>n/2) k=n-k;
	if (k==1)  return n;
	if (k==0)  return 1;
	unsigned long long r;
	if (n+k>=90) {
		// разрядности может не хватить, используем рекурсию
		r=bcl(n-1,k);
		r+=+bcl(n-1,k-1);
	}
	else {
		r=1;
		for (int i=1; i<=k;++i) {
			r*=n-k+i;
			r/=i;
		}
	}
	return r;
}

In some of Olympiad tasks I needed to calculate a lot of C (n, k) for n> 70, i.e. they obviously did not get into unsigned long long. Naturally, it was necessary to use "long arithmetics" (). For this task I wrote "a recursion with memory": the calculated coefficients were remembered in an array and the exponential growth of time of calculation was not.

Addition after the publication


At discussion options with storing of the calculated values are often mentioned. I have a code with dynamic memory allocation, but I did not bring him. For this moment here the simplest and effective chersanya from the comment: http://habrahabr.ru/post/274689/#comment_8730359http://habrahabr.ru/post/274689/#comment_8730359

unsigned bcr_cache[N_MAX][K_MAX] = {0};

unsigned bcr(int n,int k) 
{
	if (k>n/2) k=n-k;
	if (k==1)  return n;
	if (k==0)  return 1;
        if (bcr_cache[n][k] == 0)
            bcr_cache[n][k] = bcr(n-1,k)+bcr(n-1,k-1);
        return bcr_cache[n][k];
}


If in the program it is necessary to use [almost] all coefficients of a triangle of Pascal (to any level), then the provided recursive algorithm with storing — the fastest method. The similar code suits also for unsigned long long and even for long arithmetics (though there, probably, it is better to calculate and select dynamically required memory size). Specific N_MAX values can be such:
35 — will consider all coefficients of C (n, k),< 35 для unsigned int (32 бита)
n 68 — will consider all coefficients of C (n, k),< 68 для unsigned long long (64 бита)
n 200 — will consider coefficients of C (n, k), < 200 и некоторых k для unsigned long long (64 бита). Например, для С(82,21)=~1.83*10n19
K_MAX — it can be N_MAX/2+1, but no more than 35 as C (68,34) abroad unsigned long long.
For simplicity it is always possible to take K_MAX=35 and not to think "will enter — will not enter" (until we pass to numbers with digit capacity> 64 bits).

This article is a translation of the original post at habrahabr.ru/post/274689/
If you have any questions regarding the material covered in the article above, please, contact the original author of the post.
If you have any complaints about this article or you want this article to be deleted, please, drop an email here: sysmagazine.com@gmail.com.

We believe that the knowledge, which is available at the most popular Russian IT blog habrahabr.ru, should be accessed by everyone, even though it is poorly translated.
Shared knowledge makes the world better.
Best wishes.

comments powered by Disqus