You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 

72 lines
2.2 KiB

#ifndef ALGOIM_BINOMIAL_HPP
#define ALGOIM_BINOMIAL_HPP
// algoim::Binomial implements a few simple methods to calculate binomial
// coefficients. For especially large coefficients, there are better methods
// (faster, more stable) than the ones used here, but these are sufficient
// for their application in the Algoim library
#include <array>
#include <vector>
#include <unordered_map>
#include "real.hpp"
namespace algoim
{
struct Binomial
{
// Compute a row of binomial coefficients binom(n,i), i = 0,...,n;
// it is assumed that 'out' has length at least n+1
static void compute_row(int n, real* out)
{
out[0] = 1;
if (n == 0)
return;
out[1] = n;
for (int k = 2; k <= n/2; ++k)
out[k] = (out[k - 1] * (n + 1 - k)) / k;
for (int k = 0; k <= n/2; ++k)
out[n-k] = out[k];
}
// Compute (and cache) all binomial coefficients binom(n,i), i = 0,...,n;
// for n small enough, the row is obtained directly from a dense table;
// for n large, the row is computed and cached in a sparse hash table
static const real* row(int n)
{
static constexpr int m = 31;
static const auto precomputed = []()
{
std::array<real,((m+1)*(m+2))/2> d;
int i = 0;
for (int n = 0; n <= m; ++n)
{
compute_row(n, &d[i]);
i += n + 1;
}
return d;
}();
if (n <= m)
return &precomputed[ (n*(n+1))/2 ];
// All other n are stored in a thread_local hash table
static thread_local std::unordered_map<int,std::vector<real>> coeff;
auto& b = coeff[n];
if (b.empty())
{
b.resize(n + 1);
compute_row(n, b.data());
}
return b.data();
};
// returns the binomial coefficient (n \\ i)
static real c(int n, int i)
{
return row(n)[i];
}
};
} // namespace algoim
#endif