1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
|
// ************************************************************************************************
//
// BornAgain: simulate and fit reflection and scattering
//
//! @file Base/Math/Numeric.cpp
//! @brief Implements "almost equal" in namespace Numeric.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
// ************************************************************************************************
#include "Base/Math/Numeric.h"
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <limits>
//! For tests only.
void check_scalar(double a, double b, int ulp)
{
if (!Numeric::almostEqual(a, b, ulp))
exit(1);
}
double Numeric::relativeDifference(double a, double b)
{
return 2 * std::abs(a - b) / (std::abs(a) + std::abs(b));
}
bool Numeric::almostEqual(double a, double b, int ulp)
{
constexpr double eps = std::numeric_limits<double>::epsilon();
return std::abs(a - b) <= eps * ulp * (std::abs(a) + std::abs(b)) / 2;
}
bool Numeric::almostEqual(complex_t a, complex_t b, int ulp)
{
constexpr double eps = std::numeric_limits<double>::epsilon();
return std::abs(a - b) <= eps * ulp * (std::abs(a) + std::abs(b)) / 2;
}
bool Numeric::almostEqual(const R3& a, const R3& b, int ulp)
{
return almostEqual(a.x(), b.x(), ulp) && almostEqual(a.y(), b.y(), ulp)
&& almostEqual(a.z(), b.z(), ulp);
}
double Numeric::ignoreDenormalized(double value)
{
return (std::fpclassify(value) == FP_SUBNORMAL) ? 0.0 : value;
}
double Numeric::round_decimal(double val, double digits) // copied from Frida/trivia
{
if (val == 0 || digits <= 0)
return 0;
double v = fabs(val);
int n = (int)floor(digits - log10(v));
// cout << " n=" << n;
double p = pow(10., n);
double r = round(v * p) / p;
return (val < 0 ? -1 : +1) * r;
}
int Numeric::orderOfMagnitude(double x)
{
if (x == 0)
return 0;
return std::floor(std::log10(std::abs(x)));
}
|
