Kirjoittaja: jlaire
Kirjoitettu: 25.10.2015 – 25.10.2015
Tagit: algoritmit, matematiikka, koodi näytille, vinkki
Tämä vinkki esittelee erilaisia toteutuksia Pascalin kolmiosta tuttujen binomikertoimien laskemiseen. Niille on käyttöä esimerkiksi kisaohjelmoinnissa.
Koodissa keskitytään toimivuuteen ja selkeyteen eikä juurikaan tehokkuuteen. Koodit käyttävät C++11:n uusia ominaisuuksia, mikä täytyy huomioida koodin kääntämisessä ainakin vanhemmilla kääntäjillä.
Aloitetaan yhtälöstä C(n,k) = n!/(k!(n-k)!), joka on helppo muuttaa koodiksi:
template <typename Int>
Int factorial(Int n) {
Int result = 1;
for (Int i = 1; i <= n; ++i) {
result *= i;
}
return result;
}
template <typename Int>
Int binomial_with_factorial(Int n, Int k) {
return factorial(n) / (factorial(k) * factorial(n - k));
}Kertoma kuitenkin kasvaa huomattavasti binomikertoimia nopeammin: C(21,10) = 352716, mutta 21! > 264, joten funktio palauttaa tähän tapaukseen väärän vastauksen, vaikka käytettäisiin 64-bittisiä kokonaislukuja. Tämä on hyvä pitää mielessä, vaikka kielessä olisi sisäänrakennettuna mielivaltaisen tarkkuuden kokonaisluvut; niiden käyttäminen turhaan voi hidastaa koodia huomattavasti.
Seuraava toteutus käyttää dynaamista ohjelmointia, joka pohjautuu sääntöön C(n,k) = C(n-1,k) + C(n-1,k-1). Tämän koodin suorituksessa väliarvot eivät koskaan ole lopputulosta suurempia, joten ylivuodot eivät ole ongelma.
#include <map>
template <typename Int>
Int binomial_with_dp(Int n, Int k) {
static std::map<std::pair<Int, Int>, Int> cache{{{0, 0}, 1}};
if (k < 0 || k > n) {
return 0;
}
std::pair<Int, Int> key{n, k};
if (!cache[key]) {
cache[key] = binomial_with_dp<Int>(n - 1, k) + binomial_with_dp<Int>(n - 1, k - 1);
}
return cache[key];
}Tämä lienee useimpiin käyttötapauksiin ihan hyvä ratkaisu. Joskus halutaan kuitenkin välttyä pienempien arvojen tallentamiselta. Palataan takaisin kertomakaavaan, mutta esitetään se seuraavassa muodossa:
n-k+1 n-k+2 n-k+k
C(n,k) = ----- * ----- * ... * -----
1 2 kVäliarvot pysyvät paljon pienempinä, kun kerto- ja jakolaskuja tehdään vuorotellen. Hyödyntämällä sääntöä C(n,k) = C(n,n-k) voidaan varmistaa, että tapaukset kuten C(63,62) toimivat yhtä hyvin kuin pienet k:n arvot.
#include <algorithm>
template <typename Int>
bool normalize_k(Int n, Int& k) {
if (k < 0 || k > n) {
return false;
}
k = std::min<Int>(k, n - k);
return true;
}
template <typename Int>
Int binomial_with_loop(Int n, Int k) {
if (!normalize_k(n, k)) {
return 0;
}
Int result = 1;
for (Int i = 1, a = n - k + 1; i <= k; ++a, ++i) {
result *= a;
result /= i;
}
return result;
}Esimerkiksi arvon C(63,28) = 629308289804197437 laskeminen menee oikein, mikä on huomattava parannus kertomatoteutukseen! Parametreilla (n,k) = (63,29) tulee kuitenkin väärä tulos. Syy tähän on se, että silmukan sisällä tehdään kertolasku ensin. Tämä on yleensä hyvä käytäntö kokonaislukujen kanssa, koska se välttää pyöritysvirheet jakolaskussa. Se voi kuitenkin aiheuttaa ylivuodon ja siten väärän lopputuloksen, kuten tässä. Järjestyksen muuttaminen ei suoraan onnistu, koska se johtaisi pyöristysvirheisiin.
Pienellä päättelyllä voidaan osoittaa, että jakolaskun tekeminen ensin ei aiheuta ongelmia, jos murtoluvut supistetaan.
template <typename Int>
constexpr Int gcd(Int a, Int b) {
return b == 0 ? a : gcd<Int>(b, a % b);
}
template <typename Int>
Int binomial_with_loop_and_gcd(Int n, Int k) {
if (!normalize_k(n, k)) {
return 0;
}
Int result = 1;
for (Int i = 1, a = n - k + 1; i <= k; ++a, ++i) {
Int div = gcd(a, i);
result /= i / div;
result *= a / div;
}
return result;
}Tämä toteutus antaa vihdoin oikean tuloksen kaikille tapauksille, olettaen että se mahtuu käytettyyn tyyppiin. Koska kyseessä on puhdas funktio, se voidaan ilmaista suoraviivaisesti myös constexpr-funktiona:
template <typename Int>
constexpr Int binomial_constexpr_helper(Int n, Int k, Int result, Int a, Int i) {
return
i > k ? result :
binomial_constexpr_helper<Int>(n, k, result / (i / gcd(a, i)) * (a / gcd(a, i)), a + 1, i + 1);
}
template <typename Int>
constexpr Int binomial_constexpr(Int n, Int k) {
return
k > n ? 0 :
k > n / 2 ? binomial_constexpr<Int>(n, n - k) :
binomial_constexpr_helper<Int>(n, k, 1, n - k + 1, 1);
}Funktiota binomial_constexpr(n,k) voi käyttää käännösaikasessa ohjelmoinnissa. Ennen C++11-standardia tämä olisi vaatinut kikkailua, kuten Weggon koodivinkistä näkyy. C++14 laajentaa constexpr-funktioiden määritelmää lisää, ja itse asiassa binomial_with_loop_and_gcd kelpaa sellaisenaan.
Tässä lopuksi vielä pieni testiohjelma, joka havainnollistaa kertomatoteutuksen ja ensimmäisen silmukkatoteutuksen rajoja.
#include <stdint.h>
#include <functional>
#include <iomanip>
#include <iostream>
int main() {
auto binomial_with_factorial_safe = [](uint16_t n, uint16_t k) -> uint16_t {
// Avoid division by zero
if (factorial(k) * factorial(n - k) == 0) {
return 0;
}
return binomial_with_factorial(n, k);
};
std::vector<std::pair<char, std::function<uint16_t(uint16_t, uint16_t)>>> implementations{
{'F', binomial_with_factorial_safe},
{'L', binomial_with_loop<uint16_t>},
{'G', binomial_with_loop_and_gcd<uint16_t>},
{'C', binomial_constexpr<uint16_t>},
{'R', binomial_with_dp<uint16_t>},
};
size_t I = implementations.size();
std::cout.flags(std::ios::left);
uint16_t max_n = 20;
std::cout << std::setw(I) << "";
for (uint16_t k = 0; k <= max_n / 2; ++k) {
std::cout << "|k=";
std::cout << std::setw(I - 2) << k;
}
std::cout << '\n';
for (uint16_t n = 0; n <= max_n; ++n) {
std::cout << "n=" << std::setw(I - 2) << n;
for (uint16_t k = 0; k <= max_n / 2; ++k) {
std::cout << '|';
uint32_t correct = binomial_with_dp<uint32_t>(n, k);
if (correct > std::numeric_limits<uint16_t>::max()) {
std::cout << 'b' << std::string(I - 2, 'z') << 't';
}
else if (k > n / 2) {
std::cout << std::setw(I) << "";
}
else {
for (const auto& p : implementations) {
if (p.second(n, k) == correct) {
std::cout << p.first;
}
else {
std::cout << '.';
}
}
}
}
std::cout << '\n';
}
}Ohjelman tulostus on alla. Kertomakaava (F) toimii oikein vain tapauksille n <= 8. Silmukka (L) toimii oikein lähes aina, mutta antaa väärän tuloksen esimerkiksi kun (n,k) = (17,6). Muut funktiot (G,C,R) toimivat aina oikein, paitsi jos vastausta on mahdoton esittää annetulla tyypillä (bzzzt).
|k=0 |k=1 |k=2 |k=3 |k=4 |k=5 |k=6 |k=7 |k=8 |k=9 |k=10 n=0 |FLGCR| | | | | | | | | | n=1 |FLGCR| | | | | | | | | | n=2 |FLGCR|FLGCR| | | | | | | | | n=3 |FLGCR|FLGCR| | | | | | | | | n=4 |FLGCR|FLGCR|FLGCR| | | | | | | | n=5 |FLGCR|FLGCR|FLGCR| | | | | | | | n=6 |FLGCR|FLGCR|FLGCR|FLGCR| | | | | | | n=7 |FLGCR|FLGCR|FLGCR|FLGCR| | | | | | | n=8 |FLGCR|FLGCR|FLGCR|FLGCR|FLGCR| | | | | | n=9 |.LGCR|.LGCR|.LGCR|.LGCR|.LGCR| | | | | | n=10 |.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR| | | | | n=11 |.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR| | | | | n=12 |.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR| | | | n=13 |.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR| | | | n=14 |.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR| | | n=15 |.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR| | | n=16 |.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|..GCR|..GCR| | n=17 |.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|..GCR|..GCR|..GCR| | n=18 |.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|..GCR|..GCR|..GCR|..GCR| n=19 |.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|..GCR|..GCR|bzzzt|bzzzt|bzzzt n=20 |.LGCR|.LGCR|.LGCR|.LGCR|.LGCR|..GCR|..GCR|bzzzt|bzzzt|bzzzt|bzzzt
Hieno vinkki! Esittelee hienosti aiheen ja C++11:n uusia ominaisuuksia. Itselläni vaan pisti silmään tuon templaten nimi ("Int"), kun sen näyttää syntaksivirheeltä.
Tyyppiparametrien nimet aloitetaan yleensä isolla kirjaimella. Hyvin geneerisissä tapauksissa käytetään lyhyitä nimiä kuten T, U, X, T1, T2 yms., mutta kun parametrin oletetaan olevan tietynlainen tyyppi, nimeksi annetaan sitä kuvaava sana kuten vaikka Alloc, InputIterator tai Predicate. Ajattelin että Int olisi tämän konvention mukainen. (Kun konseptit vihdoin ja viimein lisätään C++:aan, tyyppiparametrille tehdyt oletukset voi ilmaista kätevästi.)
Toisaalta etsin C++14-standardista juuri esimerkkejä, ja siellä samantapaisissa tilanteissa nimenä on charT. Ehkä intT olisi tässä parempi kuin Int.