Commit a51084cd authored by Mikaël Salson's avatar Mikaël Salson

tools: nChoosek stores some computations to avoid recomputing them again and again.

We have a 500x500 array (half of it is useless…) which takes about 2MB
but avoid lots of computations when computing the e-value.
parent 4484f619
......@@ -218,16 +218,25 @@ string reverse(const string &text) {
return string(text.rbegin(), text.rend());
}
double nChoosek_stored[NB_N_CHOOSE_K_STORED][NB_N_CHOOSE_K_STORED] = {};
double nChoosek(unsigned n, unsigned k)
{
if (k > n) return 0;
if (k * 2 > n) k = n-k;
if (k == 0) return 1;
double result = n;
for( unsigned i = 2; i <= k; ++i ) {
result *= (n-i+1);
result /= i;
if (n >= NB_N_CHOOSE_K_STORED || nChoosek_stored[n][k] == 0) {
double result = 1;
unsigned i;
for (i = 0; i < k && ((n-i) >= NB_N_CHOOSE_K_STORED || nChoosek_stored[n-i][k-i] == 0); i++ ) {
result *= (n - i)*1./(k - i);
}
if (i < k) {
result *= nChoosek_stored[n-i][k-i];
}
if (n < NB_N_CHOOSE_K_STORED)
nChoosek_stored[n][k] = result;
return result;
}
return result;
return nChoosek_stored[n][k];
}
......@@ -22,6 +22,8 @@ using namespace std;
#define SEED_S12 "######-######"
#define SEED_S13 "#######-######"
#define NB_N_CHOOSE_K_STORED 500
string seed_contiguous(int k);
int seed_weight(const string &seed);
......@@ -149,6 +151,7 @@ string reverse(const string &text);
*/
Sequence create_sequence(string label_full, string label, string sequence, string quality);
extern double nChoosek_stored[NB_N_CHOOSE_K_STORED][NB_N_CHOOSE_K_STORED];
/**
* @return the combinatorial of k among n
* @see http://stackoverflow.com/a/9331125/1192742
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment