Calcolo dell'importo delle combinazioni
-
11-09-2019 - |
Domanda
Saluti,
So che puoi ottenere la quantità di combinazioni con la seguente formula (senza ripetizioni e l'ordine non è importante):
// Choose r from n n! / r!(n - r)!
Tuttavia, non so come implementarlo in C++, poiché ad esempio with
n = 52 n! = 8,0658175170943878571660636856404e+67
il numero diventa troppo grande anche per unsigned __int64
(O unsigned long long
).Esiste qualche soluzione alternativa per implementare la formula senza librerie "bigint" di terze parti?
Soluzione
Ecco un antico algoritmo che è esatto e non trabocchi a meno che il risultato è quello di grande per un long long
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;
}
Questo algoritmo è anche in Knuth "The Art of Computer Programming, 3rd Edition, Volume 2: Seminumerical Algoritmi" Penso
.UPDATE: C'è una piccola possibilità che l'algoritmo traboccherà sulla linea:
r *= n--;
molto di grandi dimensioni n. Un naive limite superiore è sqrt(std::numeric_limits<long long>::max())
che significa un n
inferiore rougly 4.000.000.000 di.
Altri suggerimenti
Ecco un antico algoritmo che è esatto e non trabocchi a meno che il risultato è quello di grande per un
long long
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; }
Questo algoritmo è anche in Knuth "The Art of Computer Programming, 3rd Edition, Volume 2: Seminumerical Algoritmi" Penso
.UPDATE: C'è una piccola possibilità che l'algoritmo traboccherà sulla linea:
r *= n--;
molto di grandi dimensioni n. Un naive limite superiore è
sqrt(std::numeric_limits<long long>::max())
che significa unn
inferiore rougly 4.000.000.000 di.
Consideriamo n == 67 e k == 33. È possibile che questo algoritmo trabocca di una lunga lunga 64 bit senza segno. Eppure la risposta corretta è rappresentabile in 64 bit: 14.226.520.737.620.288.370. E l'algoritmo di cui sopra tace circa la sua troppo pieno, scegliere (67, 33) restituisce:
8.829.174.638.479.413
Una risposta credibile ma non corretto.
Tuttavia l'algoritmo di cui sopra può essere leggermente modificato per non traboccare fino a quando la risposta finale è rappresentabile.
Il trucco sta nel riconoscere che ad ogni iterazione, la divisione r / d è esatta. Temporaneamente riscrittura:
r = r * n / d;
--n;
Per questo per essere esatti, significa che se si Estesa R, N e D nelle loro fattorizzazioni primi, allora si potrebbe facilmente annullare d, e lasciato con un valore modificato per n, lo chiamano t, e quindi il calcolo di R è semplicemente:
// compute t from r, n and d
r = r * t;
--n;
Un modo semplice e veloce per farlo è quello di trovare il massimo comun divisore di R & S, lo chiamano g:
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;
Ora possiamo fare la stessa cosa con d_temp e n (trovare il massimo comune divisore). Tuttavia poiché sappiamo a priori che R * n / d è esatto, allora sappiamo anche che MCD (d_temp, n) == d_temp, e quindi non abbiamo bisogno di calcolarlo. Così possiamo dividere n per d_temp:
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;
Pulizia:
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;
}
Ora è possibile calcolare scegliere (67, 33) senza troppopieno. E se si tenta scegliere (68, 33), si otterrà un'eccezione invece di una risposta sbagliata.
Il seguente procedura calcolerà la n-scegliere-k, usando la definizione ricorsiva e Memoizzazione. La routine è estremamente rapida e precisa:
inline unsigned long long n_choose_k(const unsigned long long& n,
const unsigned long long& k)
{
if (n < k) return 0;
if (0 == n) return 0;
if (0 == k) return 1;
if (n == k) return 1;
if (1 == k) return n;
typedef unsigned long long value_type;
value_type* table = new value_type[static_cast<std::size_t>(n * n)];
std::fill_n(table,n * n,0);
class n_choose_k_impl
{
public:
n_choose_k_impl(value_type* table,const value_type& dimension)
: table_(table),
dimension_(dimension)
{}
inline value_type& lookup(const value_type& n, const value_type& k)
{
return table_[dimension_ * n + k];
}
inline value_type compute(const value_type& n, const value_type& k)
{
if ((0 == k) || (k == n))
return 1;
value_type v1 = lookup(n - 1,k - 1);
if (0 == v1)
v1 = lookup(n - 1,k - 1) = compute(n - 1,k - 1);
value_type v2 = lookup(n - 1,k);
if (0 == v2)
v2 = lookup(n - 1,k) = compute(n - 1,k);
return v1 + v2;
}
value_type* table_;
value_type dimension_;
};
value_type result = n_choose_k_impl(table,n).compute(n,k);
delete [] table;
return result;
}
Ricordati che
n! / ( n - r )! = n * ( n - 1) * .. * (n - r + 1 )
quindi è molto più piccolo di n!.Quindi la soluzione è valutare n* ( n - 1 ) * ...* ( n - r + 1) invece di calcolare prima n!e poi dividendolo.
Ovviamente tutto dipende dalla grandezza relativa di n e r: se r è relativamente grande rispetto a n, allora non andrà comunque bene.
Bene, devo rispondere alla mia domanda. Stavo leggendo su triangolo di Pascal e per caso notato che siamo in grado di calcolare la quantità di combinazioni con esso:
#include <iostream>
#include <boost/cstdint.hpp>
boost::uint64_t Combinations(unsigned int n, unsigned int r)
{
if (r > n)
return 0;
/** We can use Pascal's triange to determine the amount
* of combinations. To calculate a single line:
*
* v(r) = (n - r) / r
*
* Since the triangle is symmetrical, we only need to calculate
* until r -column.
*/
boost::uint64_t v = n--;
for (unsigned int i = 2; i < r + 1; ++i, --n)
v = v * n / i;
return v;
}
int main()
{
std::cout << Combinations(52, 5) << std::endl;
}
Come la fattorizzazione prima del coefficiente binomiale è probabilmente il modo più efficiente per calcolarlo, soprattutto se la moltiplicazione è costoso. Questo è certamente vero per il relativo problema del calcolo fattoriale (vedi Clicca qui per esempio ).
Ecco un semplice algoritmo basato sul crivello di Eratostene che calcola la fattorizzazione prima. L'idea è sostanzialmente quella di passare attraverso i numeri primi, come li trovate usando il setaccio, ma poi anche per calcolare quanti di loro multipli rientrano negli intervalli [1, k] e [n-k + 1, n]. Il setaccio è essenzialmente un O (n \ log \ log n) algoritmo, ma non c'è la moltiplicazione fatto. Il numero effettivo di moltiplicazioni necessarie una volta che la fattorizzazione prima si trova è nella peggiore delle ipotesi O \ left (\ frac {n \ log \ log n} {\ log n} \ right), e probabilmente ci sono modi più veloce di quello.
prime_factors = []
n = 20
k = 10
composite = [True] * 2 + [False] * n
for p in xrange(n + 1):
if composite[p]:
continue
q = p
m = 1
total_prime_power = 0
prime_power = [0] * (n + 1)
while True:
prime_power[q] = prime_power[m] + 1
r = q
if q <= k:
total_prime_power -= prime_power[q]
if q > n - k:
total_prime_power += prime_power[q]
m += 1
q += p
if q > n:
break
composite[q] = True
prime_factors.append([p, total_prime_power])
print prime_factors
Semplificare la formula, però. Se non si vuole fare una divisione.
Una delle via più breve:
int nChoosek(int n, int k){
if (k > n) return 0;
if (k == 0) return 1;
return nChoosek(n - 1, k) + nChoosek(n - 1, k - 1);
}
Utilizzando un trucco sporco con un lungo doppio, è possibile ottenere la stessa accuratezza Howard Hinnant (e probabilmente più):
unsigned long long n_choose_k(int n, int k)
{
long double f = n;
for (int i = 1; i<k+1; i++)
f /= i;
for (int i=1; i<k; i++)
f *= n - i;
unsigned long long f_2 = std::round(f);
return f_2;
}
L'idea è quella di dividere la prima volta da k! e poi moltiplicare per n (n-1) ... (n-k + 1). L'approssimazione attraverso il doppio può essere evitato invertendo l'ordine del ciclo.
Se vuoi essere sicuro al 100% che non trabocchi si verificano fino a quando il risultato finale è entro il limite numerico, è possibile riassumere di Pascal Triangolo riga per riga:
for (int i=0; i<n; i++) {
for (int j=0; j<=i; j++) {
if (j == 0) current_row[j] = 1;
else current_row[j] = prev_row[j] + prev_row[j-1];
}
prev_row = current_row; // assume they are vectors
}
// result is now in current_row[r-1]
Tuttavia, questo algoritmo è molto più lento del moltiplicazione uno. Così forse si potrebbe usare la moltiplicazione per generare tutti i casi si sa che sono 'sicuro' e poi usare aggiunta da lì. (.. o si potrebbe utilizzare una libreria BigInt).