NunoSempere 2023-07-23 12:53:46 +02:00
parent 9e1d4ee6d4
commit 32033b5c86
2 changed files with 5 additions and 4 deletions

Binary file not shown.

View File

@ -105,9 +105,9 @@ float sample_gamma(float alpha, uint64_t* seed)
v = 1.0 + c * x; v = 1.0 + c * x;
} while (v <= 0.0); } while (v <= 0.0);
v = pow(v, 3); v = v * v * v;
u = sample_unit_uniform(seed); u = sample_unit_uniform(seed);
if (u < 1.0 - 0.0331 * pow(x, 4)) { // Condition 1 if (u < 1.0 - 0.0331 * (x * x * x * x)) { // Condition 1
// the 0.0331 doesn't inspire much confidence // the 0.0331 doesn't inspire much confidence
// however, this isn't the whole story // however, this isn't the whole story
// by knowing that Condition 1 implies condition 2 // by knowing that Condition 1 implies condition 2
@ -115,7 +115,7 @@ float sample_gamma(float alpha, uint64_t* seed)
// i.e., of not using the logarithms // i.e., of not using the logarithms
return d * v; return d * v;
} }
if (log(u) < 0.5 * pow(x, 2) + d * (1.0 - v + log(v))) { // Condition 2 if (log(u) < 0.5 * (x * x) + d * (1.0 - v + log(v))) { // Condition 2
return d * v; return d * v;
} }
} }
@ -161,7 +161,8 @@ float array_std(float* array, int length)
float mean = array_mean(array, length); float mean = array_mean(array, length);
float std = 0.0; float std = 0.0;
for (int i = 0; i < length; i++) { for (int i = 0; i < length; i++) {
std += pow(array[i] - mean, 2.0); std += (array[i] - mean);
std *= std;
} }
std = sqrt(std / length); std = sqrt(std / length);
return std; return std;