I have to find a recursive function in C to compute big Binomial Coefficients. e.g. 59C10 I've written the below code but takes too much time. Is there a better way to do it ?
long long nCk(long long n, long long k)
{
if(n == k || k == 0)
{
return 1;
}
else if(k == 1)
{
return n;
}
else
{
return nCk(n-1,k-1) nCk(n-1,k);
}
}
CodePudding user response：
First of all, I don't think recursion is necessarily the right way to do this.
However, you can make the calculation take less than a second if you use memoization, like Eric Postpischil suggested.
The program below uses a naive and unnecessarily large cache to calculate 59C10. It outputs 62828356305, which seems to be correct:
#include <stdio.h>
#include <stdint.h>
uint64_t cache[60][60];
uint64_t nCk(uint64_t n, uint64_t k)
{
if (k > n)
{
return 0;
}
else if (k == n || k == 0)
{
return 1;
}
else if (cache[n][k])
{
return cache[n][k];
}
else
{
return cache[n][k] = nCk(n-1, k-1) nCk(n-1, k);
}
}
int main()
{
printf("%lld\n", nCk(59, 10));
}
The cache is actually way bigger than it needs to be: if we restructure things, it would only need to hold two rows of the binomial triangle at a time (the previously-computed one and the one we are working on).
Edit: I know you already accepted my answer, but here is a more sophisticated answer that uses a much smaller cache, provided by the caller of the function:
#include <stdio.h>
#include <stdint.h>
// Calculates (n choose k) for all values of k from k_min to k_max inclusive.
// Stores the results in output[k_min] to output[k_max].
// The 'output' and 'extra' buffers must each be large enough to hold
// k_max 1 elements.
void nCkCore(uint64_t n, uint64_t k_min, uint64_t k_max,
uint64_t * output, uint64_t * extra)
{
// Take care of the base case.
if (n == 0)
{
for (uint64_t k = k_min; k <= k_max; k )
{
output[k] = k == 0;
}
return;
}
// Optimization: Take care of any trivial outputs and reduce k_max
// if possible, so there is less work to do in the deeper recursive calls.
while (k_max > n) { output[k_max--] = 0; }
// Populate 'extra' with the values we need from the preivous
// row in the binomial triangle.
uint64_t k2_min = k_min > 0 ? k_min - 1 : k_min;
nCkCore(n - 1, k2_min, k_max, extra, output);
// Populate 'output'.
for (uint64_t k = k_min; k <= k_max; k )
{
output[k] = k == 0 ? 1 : extra[k - 1] extra[k];
}
}
// Cache must be large enough to hold 2*(n 1) elements.
uint64_t nCk(uint64_t n, uint64_t k, uint64_t * cache)
{
if (k > n) { return 0; }
nCkCore(n, k, k, cache, cache (n 1));
return cache[k];
}
int main()
{
uint64_t cache[1024];
printf("%lld\n", nCk(59, 10, cache));
}
CodePudding user response：
Consider:
That can be directly implemented as:
/* The helper assumes all sanity checks have been made */
static long long nCk_helper(int n, int k) {
if (k == 0) return 1;
return nCk_helper(n - 1, k - 1) * n / k;
}
long long nCk(int n, int k) {
if (n < k || k < 0)
return 0; /* Or some error value */
/* Take the shorter of the possible computations */
if (k <= n - k)
return nCk_helper(n, k);
else
return nCk_helper(n, n - k);
}
CodePudding user response：
You can use "better" maths and (still) make it recursive. The famous:
n!
----------
k! (n-k)!
is actually a simple fraction e.g. (20 6) becomes:
20 19 18 17 16 15
-----------------
6 5 4 3 2 (1)
which is 20/6 times (19 5).
The function:
long double binom_rec(long double n, long double k) {
if (k == 1)
return n;
return (n / k) * binom_rec(n - 1, k - 1);
}
With only double the result for C(59 10) was _4.9997.
To get rid of the mass of numbers, a cancellation algorithm is needed. int is enough. Here an illustration for (59 10):
59 58 57 56 55 54 53 52 51 50
10 9 8 7 6 5 4 3 2 .
59 29 19 7 11 6 53 13 51 5 (After 1 pass)
. . . 7 6 . . . . .
59 29 19 . 11 . 53 13 51 5 (Second pass, done)
. . . . . . . . . .
This still overflows long ints above (60 30). So at this point one should switch to floats. At even higher pascal numbers the bottom row does not get empty, so one division is possible. As in C(590 100):
$ ./a.out 590 100
590 589 588 587 586 585 584 583 582 581 580 579 578 577 576 575 574 573 572 571 570 569 568 567 566 565 564 563 562 561 560 559 558 557 556 555 554 553 552 551 550 549 548 547 546 545 544 543 542 541 540 539 538 537 536 535 534 533 532 531 530 529 528 527 526 525 524 523 522 521 520 519 518 517 516 515 514 513 512 511 510 509 508 507 506 505 504 503 502 501 500 499 498 497 496 495 494 493 492 491
100 99 98 97 96 95 94 93 92 91 90 89 88 87 86 85 84 83 82 81 80 79 78 77 76 75 74 73 72 71 70 69 68 67 66 65 64 63 62 61 60 59 58 57 56 55 54 53 52 51 50 49 48 47 46 45 44 43 42 41 40 39 38 37 36 35 34 33 32 31 30 29 28 27 26 25 24 23 22 21 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 .
10 19 6 587 293 9 8 11 6 7 10 193 17 577 6 23 7 573 11 571 6 569 8 7 566 113 6 563 562 11 7 13 6 557 139 15 554 7 6 19 10 9 548 547 6 545 8 543 542 541 6 7 538 537 8 535 6 13 7 59 53 23 6 31 526 7 524 523 6 521 13 519 7 11 6 515 514 9 8 73 6 509 508 13 11 505 6 503 502 501 5 499 83 497 8 5 13 493 41 491
. . . . . . . . . . . . . . . . . . . . . . 78 . . . . . 72 . 70 69 . . 66 . . 63 . . 60 . . . 56 . 54 . . . 50 49 48 . . 45 44 . 42 . . . . . 36 35 . 33 32 . 30 . 28 27 26 . 24 . 22 21 20 19 18 . 16 15 14 13 . 11 . . 8 . . . . . . .
. . . 587 293 . . . . . . 193 17 577 . . . 191 . 571 . 569 . . 566 113 . 563 562 . . . . 557 139 . 554 . 2 19 . . 137 547 . 109 . 181 542 541 . . 538 537 . 535 . . . 59 53 23 3 31 526 . 131 523 . 521 . 519 . . 3 103 514 . 2 73 2 509 127 13 . 101 . 503 502 167 . 499 83 497 2 . 13 493 41 491
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 . . . . . . . . . . 16 3 14 . . . . . 8 . . . . . . .
. . . 587 293 . . . . . . 193 17 577 . . . 191 . 571 . 569 . . 566 113 . 563 562 . . . . 557 139 . 554 . . 19 . . 137 547 . 109 . 181 542 541 . . 538 179 . 535 . . . 59 53 23 . 31 526 . 131 523 . 521 . 519 . . . 103 514 . . 73 . 509 127 13 . 101 . 503 251 167 . 499 83 497 . . 13 493 41 491
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 . . . . . . . . . . . . 7 . . . . . 8 . . . . . . .
. . . 587 293 . . . . . . 193 17 577 . . . 191 . 571 . 569 . . 566 113 . 563 562 . . . . 557 139 . 554 . . 19 . . 137 547 . 109 . 181 542 541 . . 538 179 . 535 . . . 59 53 23 . 31 526 . 131 523 . 521 . 173 . . . 103 514 . . 73 . 509 127 13 . 101 . 503 251 167 . 499 83 71 . . 13 493 41 491
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 . . . . . . .
--> 141436416421996337338823235258671000825125009186271742189481844402257617684227520675461582943717912848062992491741184.000000 / 8.000000 -->
C(590 100) = 17679552052749542167352904407333875103140626148283967773685230550282202210528440084432697867964739106007874061467648.000
C(590 100) = 1.77e 115
(scroll for the /8) -------->
Fascinating that there cannot be any remainder in this division, mathematically. In other words: my cancellation algorithm is too lazy, it does not halve 4 of the even numbers with that last 8.
Speed does not seem to be the problem. But overflow is. A non-recursive approach is cleaner in this respect:
long double binom_val_canc(int n, int k)
Now you can input "only" ints. like C(123456789 1000000), and the output is "unlimited". In between should lie some integer cancellation i.e. a real prime factor cancelling to avoid any division.
return nCk_helper(n - 1, k - 1) * n / k;
rici's answer is nice. It avoids remainders by recursion, math, and the switched ordering: the integer division has to be applied to the product, not isolated n.
First step in C(59 10) is (50*51)/2, not (51/2)*50, as I did (and the formula suggests), necessitating floating point division.
So who says recursion is no good for binomial coefficients?