StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

PO
text
Body
copied!<p>From <a href="https://stackoverflow.com/questions/1838368/calculating-the-amount-of-combinations/1838732#1838732">Andreas' answer</a>:</p> <blockquote> <p>Here's an ancient algorithm which is exact and doesn't overflow unless the result is to big for a <code>long long</code></p> <pre><code>unsigned long long choose(unsigned long long n, unsigned long long k) { if (k > n) { return 0; } unsigned long long r = 1; for (unsigned long long d = 1; d <= k; ++d) { r *= n--; r /= d; } return r; } </code></pre> <p>This algorithm is also in Knuth's "The Art of Computer Programming, 3rd Edition, Volume 2: Seminumerical Algorithms" I think.</p> <p><strong>UPDATE:</strong> There's a small possibility that the algorithm will overflow on the line:</p> <pre><code>r *= n--; </code></pre> <p>for <strong>very</strong> large n. A naive upper bound is <code>sqrt(std::numeric_limits<long long>::max())</code> which means an <code>n</code> less than rougly 4,000,000,000.</p> </blockquote> <p>Consider n == 67 and k == 33. The above algorithm overflows with a 64 bit unsigned long long. And yet the correct answer is representable in 64 bits: 14,226,520,737,620,288,370. And the above algorithm is silent about its overflow, choose(67, 33) returns:</p> <p>8,829,174,638,479,413</p> <p>A believable but incorrect answer.</p> <p>However the above algorithm can be slightly modified to never overflow as long as the final answer is representable.</p> <p>The trick is in recognizing that at each iteration, the division r/d is exact. Temporarily rewriting:</p> <pre><code>r = r * n / d; --n; </code></pre> <p>For this to be exact, it means if you expanded r, n and d into their prime factorizations, then one could easily cancel out d, and be left with a modified value for n, call it t, and then the computation of r is simply:</p> <pre><code>// compute t from r, n and d r = r * t; --n; </code></pre> <p>A fast and easy way to do this is to find the greatest common divisor of r and d, call it g:</p> <pre><code>unsigned long long g = gcd(r, d); // now one can divide both r and d by g without truncation r /= g; unsigned long long d_temp = d / g; --n; </code></pre> <p>Now we can do the same thing with d_temp and n (find the greatest common divisor). However since we know a-priori that r * n / d is exact, then we also know that gcd(d_temp, n) == d_temp, and therefore we don't need to compute it. So we can divide n by d_temp:</p> <pre><code>unsigned long long g = gcd(r, d); // now one can divide both r and d by g without truncation r /= g; unsigned long long d_temp = d / g; // now one can divide n by d/g without truncation unsigned long long t = n / d_temp; r = r * t; --n; </code></pre> <p>Cleaning up:</p> <pre><code>unsigned long long gcd(unsigned long long x, unsigned long long y) { while (y != 0) { unsigned long long t = x % y; x = y; y = t; } return x; } unsigned long long choose(unsigned long long n, unsigned long long k) { if (k > n) throw std::invalid_argument("invalid argument in choose"); unsigned long long r = 1; for (unsigned long long d = 1; d <= k; ++d, --n) { unsigned long long g = gcd(r, d); r /= g; unsigned long long t = n / (d / g); if (r > std::numeric_limits<unsigned long long>::max() / t) throw std::overflow_error("overflow in choose"); r *= t; } return r; } </code></pre> <p>Now you can compute choose(67, 33) without overflow. And if you try choose(68, 33), you'll get an exception instead of a wrong answer.</p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload