TerrainGen coverage

terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
42 / 42
100.0%
Branches:
30 / 33
90.9%
Functions:
3 / 3
100.0%
util/special_functions.cpp
Line Branch Exec Source
1 #include "terraingen/util/special_functions.h"
2
3 #include <cmath>
4 #include <limits>
5 #include <numbers>
6 #include <stdexcept>
7
8 namespace terraingen {
9
10 namespace {
11
12 constexpr int kMaxIterations = 200;
13 constexpr double kEpsilon = std::numeric_limits<double>::epsilon();
14
15 /// E_n(x) for x > 1 by its continued fraction, evaluated with Lentz's method.
16 2012149 double ExponentialIntegralFraction(int n, double x) {
17 2012149 constexpr double kTiny = std::numeric_limits<double>::min() / kEpsilon;
18 2012149 const int nm1 = n - 1;
19 2012149 double b = x + n;
20 2012149 double c = 1.0 / kTiny;
21 2012149 double d = 1.0 / b;
22 2012149 double h = d;
23
1/2
✓ Branch 0 taken 102196326 times.
✗ Branch 1 not taken.
 102196326 for (int i = 1; i <= kMaxIterations; ++i) {
24 102196326 const double a = -static_cast<double>(i) * (nm1 + i);
25 102196326 b += 2.0;
26 102196326 d = 1.0 / (a * d + b);
27 102196326 c = b + a / c;
28 102196326 const double delta = c * d;
29 102196326 h *= delta;
30
2/2
✓ Branch 0 taken 2012149 times.
✓ Branch 1 taken 100184177 times.
 102196326 if (std::fabs(delta - 1.0) <= kEpsilon) {
31 2012149 break;
32 }
33 }
34 2012149 return h * std::exp(-x);
35 }
36
37 /// E_n(x) for 0 < x <= 1 by its power series.
38 1903802 double ExponentialIntegralSeries(int n, double x) {
39 1903802 const int nm1 = n - 1;
40
2/2
✓ Branch 0 taken 1903799 times.
✓ Branch 1 taken 3 times.
 1903802 double sum = nm1 != 0 ? 1.0 / nm1 : -std::log(x) - std::numbers::egamma;
41 1903802 double factor = 1.0;
42
1/2
✓ Branch 0 taken 29410797 times.
✗ Branch 1 not taken.
 29410797 for (int i = 1; i <= kMaxIterations; ++i) {
43 29410797 factor *= -x / i;
44 29410797 double delta = 0.0;
45
2/2
✓ Branch 0 taken 27506998 times.
✓ Branch 1 taken 1903799 times.
 29410797 if (i != nm1) {
46 27506998 delta = -factor / (i - nm1);
47 } else {
48 1903799 double psi = -std::numbers::egamma;
49
2/2
✓ Branch 0 taken 1903800 times.
✓ Branch 1 taken 1903799 times.
 3807599 for (int k = 1; k <= nm1; ++k) {
50 1903800 psi += 1.0 / k;
51 }
52 1903799 delta = factor * (-std::log(x) + psi);
53 }
54 29410797 sum += delta;
55
2/2
✓ Branch 0 taken 1903802 times.
✓ Branch 1 taken 27506995 times.
 29410797 if (std::fabs(delta) <= std::fabs(sum) * kEpsilon) {
56 1903802 break;
57 }
58 }
59 1903802 return sum;
60 }
61
62 } // namespace
63
64 3915956 double ExponentialIntegral(int n, double x) {
65
11/12
✓ Branch 0 taken 3915955 times.
✓ Branch 1 taken 1 time.
✓ Branch 3 taken 3915955 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3915954 times.
✓ Branch 6 taken 1 time.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 3915952 times.
✓ Branch 9 taken 1 time.
✓ Branch 10 taken 1 time.
✓ Branch 11 taken 3 times.
✓ Branch 12 taken 3915953 times.
 3915956 if (n < 0 || !(std::isfinite(x)) || x < 0.0 || (x == 0.0 && n <= 1)) {
66
1/1
✓ Branch 2 taken 3 times.
 3 throw std::domain_error("The exponential integral E_n(x) needs n >= 0 and x > 0");
67 }
68
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3915952 times.
 3915953 if (n == 0) {
69 1 return std::exp(-x) / x;
70 }
71
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3915951 times.
 3915952 if (x == 0.0) {
72 1 return 1.0 / (n - 1);
73 }
74
2/2
✓ Branch 0 taken 2012149 times.
✓ Branch 1 taken 1903802 times.
 3915951 return x > 1.0 ? ExponentialIntegralFraction(n, x) : ExponentialIntegralSeries(n, x);
75 }
76
77 } // namespace terraingen
78