From d6beb76209faadc4b44f39b7e954d6301434f0b9 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (l)" <j.wuttke@fz-juelich.de>
Date: Sun, 21 Aug 2016 15:07:13 +0200
Subject: [PATCH] replace double_epsilon by direct use of <limits>
---
Core/Aggregate/FTDecayFunctions.cpp | 2 +-
Core/Aggregate/FTDistributions1D.cpp | 4 ++--
Core/Aggregate/FTDistributions2D.cpp | 4 ++--
Core/Aggregate/InterferenceFunction2DParaCrystal.cpp | 3 +--
.../InterferenceFunctionRadialParaCrystal.cpp | 3 +--
Core/HardParticle/FormFactorCone.cpp | 4 ++--
Core/HardParticle/FormFactorFullSpheroid.cpp | 10 ++++------
Core/HardParticle/FormFactorHemiEllipsoid.cpp | 8 +++-----
Core/HardParticle/FormFactorTruncatedSphere.cpp | 6 +++---
Core/HardParticle/FormFactorTruncatedSpheroid.cpp | 10 ++++------
Core/Mask/Line.cpp | 3 ++-
Core/SoftParticle/FormFactorSphereUniformRadius.cpp | 8 ++++----
Core/Tools/MathFunctions.cpp | 6 +++---
Core/Tools/Numeric.cpp | 11 ++++++-----
Core/Tools/Numeric.h | 4 ----
15 files changed, 38 insertions(+), 48 deletions(-)
diff --git a/Core/Aggregate/FTDecayFunctions.cpp b/Core/Aggregate/FTDecayFunctions.cpp
index b1b7b4b34ae..ef1562d199d 100644
--- a/Core/Aggregate/FTDecayFunctions.cpp
+++ b/Core/Aggregate/FTDecayFunctions.cpp
@@ -116,7 +116,7 @@ FTDecayFunction1DCosine::FTDecayFunction1DCosine(double omega)
double FTDecayFunction1DCosine::evaluate(double q) const
{
double qw = std::abs(q*m_omega);
- if (std::abs(qw/Pi::PI-1.0) < Numeric::double_epsilon) {
+ if (std::abs(qw/Pi::PI-1.0) < std::numeric_limits<double>::epsilon()) {
return m_omega/2.0;
}
else {
diff --git a/Core/Aggregate/FTDistributions1D.cpp b/Core/Aggregate/FTDistributions1D.cpp
index d7e32203f70..186d69662df 100644
--- a/Core/Aggregate/FTDistributions1D.cpp
+++ b/Core/Aggregate/FTDistributions1D.cpp
@@ -17,10 +17,10 @@
#include "BornAgainNamespace.h"
#include "IntegratorReal.h"
#include "MathFunctions.h"
-#include "Numeric.h"
#include "ParameterPool.h"
#include "Pi.h"
#include "RealParameter.h"
+#include <limits>
void IFTDistribution1D::print(std::ostream& ostr) const
{
@@ -94,7 +94,7 @@ FTDistribution1DCosine::FTDistribution1DCosine(double omega)
double FTDistribution1DCosine::evaluate(double q) const
{
double qw = std::abs(q*m_omega);
- if (std::abs(qw/Pi::PI-1.0) < Numeric::double_epsilon)
+ if (std::abs(qw/Pi::PI-1.0) < std::numeric_limits<double>::epsilon())
return 0.5;
return MathFunctions::sinc(qw)/(1.0-qw*qw/Pi::PI/Pi::PI);
}
diff --git a/Core/Aggregate/FTDistributions2D.cpp b/Core/Aggregate/FTDistributions2D.cpp
index d9666deb903..ae72d771eb2 100644
--- a/Core/Aggregate/FTDistributions2D.cpp
+++ b/Core/Aggregate/FTDistributions2D.cpp
@@ -17,10 +17,10 @@
#include "BornAgainNamespace.h"
#include "IntegratorReal.h"
#include "MathFunctions.h"
-#include "Numeric.h"
#include "ParameterPool.h"
#include "Pi.h"
#include "RealParameter.h"
+#include <limits>
IFTDistribution2D::IFTDistribution2D(
double coherence_length_x, double coherence_length_y, double gamma, double delta)
@@ -100,7 +100,7 @@ FTDistribution2DCone::FTDistribution2DCone(
double FTDistribution2DCone::evaluate(double qx, double qy) const
{
double scaled_q = std::sqrt(sumsq(qx,qy));
- if (scaled_q<Numeric::double_epsilon)
+ if (scaled_q<std::numeric_limits<double>::epsilon())
return 1.0 - 3.0*scaled_q*scaled_q/40.0;
auto integrator = make_integrator_real(this, &FTDistribution2DCone::coneIntegrand2);
double integral = integrator->integrate(0.0, scaled_q);
diff --git a/Core/Aggregate/InterferenceFunction2DParaCrystal.cpp b/Core/Aggregate/InterferenceFunction2DParaCrystal.cpp
index 3d64448d67e..31c87aff2d5 100644
--- a/Core/Aggregate/InterferenceFunction2DParaCrystal.cpp
+++ b/Core/Aggregate/InterferenceFunction2DParaCrystal.cpp
@@ -18,7 +18,6 @@
#include "Exceptions.h"
#include "ISampleVisitor.h"
#include "IntegratorReal.h"
-#include "Numeric.h"
#include "ParameterPool.h"
#include "Pi.h"
#include "RealParameter.h"
@@ -202,7 +201,7 @@ double InterferenceFunction2DParaCrystal::interference1D(
if (n<1) {
result = ((1.0 + fp)/(1.0 - fp)).real();
} else {
- if (std::norm(1.0-fp) < Numeric::double_epsilon )
+ if (std::norm(1.0-fp) < std::numeric_limits<double>::epsilon() )
result = nd;
// for (1-fp)*nd small, take the series expansion to second order in nd*(1-fp)
else if (std::abs(1.0-fp)*nd < 2e-4) {
diff --git a/Core/Aggregate/InterferenceFunctionRadialParaCrystal.cpp b/Core/Aggregate/InterferenceFunctionRadialParaCrystal.cpp
index 20be390db46..aeaae8e8d91 100644
--- a/Core/Aggregate/InterferenceFunctionRadialParaCrystal.cpp
+++ b/Core/Aggregate/InterferenceFunctionRadialParaCrystal.cpp
@@ -17,7 +17,6 @@
#include "BornAgainNamespace.h"
#include "Exceptions.h"
#include "ISampleVisitor.h"
-#include "Numeric.h"
#include "ParameterPool.h"
#include "RealParameter.h"
#include <limits>
@@ -83,7 +82,7 @@ double InterferenceFunctionRadialParaCrystal::evaluate(const kvector_t q) const
if (n<1) {
result = ((1.0 + fp)/(1.0 - fp)).real();
} else {
- if (std::norm(1.0-fp) < Numeric::double_epsilon ) {
+ if (std::norm(1.0-fp) < std::numeric_limits<double>::epsilon() ) {
result = nd;
}
// for (1-fp)*nd small, take the series expansion to second order in nd*(1-fp)
diff --git a/Core/HardParticle/FormFactorCone.cpp b/Core/HardParticle/FormFactorCone.cpp
index 22c0fcdafa4..27c012667b2 100644
--- a/Core/HardParticle/FormFactorCone.cpp
+++ b/Core/HardParticle/FormFactorCone.cpp
@@ -17,9 +17,9 @@
#include "BornAgainNamespace.h"
#include "Exceptions.h"
#include "MathFunctions.h"
-#include "Numeric.h"
#include "Pi.h"
#include "RealParameter.h"
+#include <limits>
//! @param radius of circular base
//! @param height of frustum
@@ -58,7 +58,7 @@ complex_t FormFactorCone::Integrand(double Z) const
complex_t FormFactorCone::evaluate_for_q(const cvector_t q) const
{
m_q = q;
- if ( std::abs(m_q.mag()) < Numeric::double_epsilon) {
+ if ( std::abs(m_q.mag()) < std::numeric_limits<double>::epsilon()) {
double R = m_radius;
double H = m_height;
double tga = std::tan(m_alpha);
diff --git a/Core/HardParticle/FormFactorFullSpheroid.cpp b/Core/HardParticle/FormFactorFullSpheroid.cpp
index 8e253543fd2..9cdfd98d94c 100644
--- a/Core/HardParticle/FormFactorFullSpheroid.cpp
+++ b/Core/HardParticle/FormFactorFullSpheroid.cpp
@@ -16,9 +16,9 @@
#include "FormFactorFullSpheroid.h"
#include "BornAgainNamespace.h"
#include "MathFunctions.h"
-#include "Numeric.h"
#include "Pi.h"
#include "RealParameter.h"
+#include <limits>
//! @param radius of the two equal axes
//! @param height total height of the spheroid, i.e. twice the radius of the third axis
@@ -51,10 +51,8 @@ complex_t FormFactorFullSpheroid::evaluate_for_q(const cvector_t q) const
double R = m_radius;
m_q = q;
- if (std::abs(m_q.mag()) <= Numeric::double_epsilon) {
+ if (std::abs(m_q.mag()) <= std::numeric_limits<double>::epsilon())
return Pi::PI2*R*R*H/3.;
- } else {
- complex_t qzH_half = H/2*q.z();
- return 4 * Pi::PI * mP_integrator->integrate(0.0, H/2.0) * exp_I(qzH_half);
- }
+ complex_t qzH_half = H/2*q.z();
+ return 4 * Pi::PI * mP_integrator->integrate(0.0, H/2.0) * exp_I(qzH_half);
}
diff --git a/Core/HardParticle/FormFactorHemiEllipsoid.cpp b/Core/HardParticle/FormFactorHemiEllipsoid.cpp
index 7baa217a803..51da0225f0c 100644
--- a/Core/HardParticle/FormFactorHemiEllipsoid.cpp
+++ b/Core/HardParticle/FormFactorHemiEllipsoid.cpp
@@ -16,9 +16,9 @@
#include "FormFactorHemiEllipsoid.h"
#include "BornAgainNamespace.h"
#include "MathFunctions.h"
-#include "Numeric.h"
#include "Pi.h"
#include "RealParameter.h"
+#include <limits>
//! @param radius_x half length of one horizontal main axes
//! @param radius_y half length of the other horizontal main axes
@@ -64,9 +64,7 @@ complex_t FormFactorHemiEllipsoid::evaluate_for_q(const cvector_t q) const
double W = m_radius_y;
double H = m_height;
- if (std::abs(m_q.mag()) <= Numeric::double_epsilon) {
+ if (std::abs(m_q.mag()) <= std::numeric_limits<double>::epsilon())
return Pi::PI2*R*W*H/3.;
- } else {
- return Pi::PI2*mP_integrator->integrate(0.,H );
- }
+ return Pi::PI2*mP_integrator->integrate(0.,H );
}
diff --git a/Core/HardParticle/FormFactorTruncatedSphere.cpp b/Core/HardParticle/FormFactorTruncatedSphere.cpp
index 41ee234b706..4e831a63d66 100644
--- a/Core/HardParticle/FormFactorTruncatedSphere.cpp
+++ b/Core/HardParticle/FormFactorTruncatedSphere.cpp
@@ -18,11 +18,11 @@
#include "Exceptions.h"
#include "Limits.h"
#include "MathFunctions.h"
-#include "Numeric.h"
#include "Pi.h"
#include "RealParameter.h"
+#include <limits>
-using namespace BornAgain;
+using namespace BornAgain;
FormFactorTruncatedSphere::FormFactorTruncatedSphere(double radius, double height)
: m_radius(radius), m_height(height)
@@ -63,7 +63,7 @@ complex_t FormFactorTruncatedSphere::Integrand(double Z) const
complex_t FormFactorTruncatedSphere::evaluate_for_q(const cvector_t q) const
{
m_q = q;
- if ( std::abs(q.mag()) < Numeric::double_epsilon) {
+ if ( std::abs(q.mag()) < std::numeric_limits<double>::epsilon()) {
double HdivR = m_height/m_radius;
return Pi::PI/3.*m_radius*m_radius*m_radius
*(3.*HdivR -1. - (HdivR - 1.)*(HdivR - 1.)*(HdivR - 1.));
diff --git a/Core/HardParticle/FormFactorTruncatedSpheroid.cpp b/Core/HardParticle/FormFactorTruncatedSpheroid.cpp
index af08c277cd9..dbda19c0d21 100644
--- a/Core/HardParticle/FormFactorTruncatedSpheroid.cpp
+++ b/Core/HardParticle/FormFactorTruncatedSpheroid.cpp
@@ -17,9 +17,9 @@
#include "BornAgainNamespace.h"
#include "Exceptions.h"
#include "MathFunctions.h"
-#include "Numeric.h"
#include "Pi.h"
#include "RealParameter.h"
+#include <limits>
FormFactorTruncatedSpheroid::FormFactorTruncatedSpheroid(
double radius, double height, double height_flattening)
@@ -68,10 +68,8 @@ complex_t FormFactorTruncatedSpheroid::evaluate_for_q(const cvector_t q) const
double fp = m_height_flattening;
m_q = q;
- if (std::abs(m_q.mag()) <= Numeric::double_epsilon) {
+ if (std::abs(m_q.mag()) <= std::numeric_limits<double>::epsilon())
return Pi::PI*R*H*H/fp*(1.-H/(3.*fp*R));
- } else {
- complex_t z_part = std::exp(complex_t(0.0, 1.0)*m_q.z()*(H-fp*R));
- return Pi::PI2 * z_part *mP_integrator->integrate(fp*R-H,fp*R );
- }
+ complex_t z_part = std::exp(complex_t(0.0, 1.0)*m_q.z()*(H-fp*R));
+ return Pi::PI2 * z_part *mP_integrator->integrate(fp*R-H,fp*R );
}
diff --git a/Core/Mask/Line.cpp b/Core/Mask/Line.cpp
index 3b6191d83c4..5ca4ae9a9f1 100644
--- a/Core/Mask/Line.cpp
+++ b/Core/Mask/Line.cpp
@@ -17,6 +17,7 @@
#include "Bin.h"
#include "Macros.h"
#include "Numeric.h"
+#include <limits>
GCC_DIAG_OFF(unused-parameter)
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
@@ -42,7 +43,7 @@ bool Line::contains(double x, double y) const
double d = distance(p, line);
- return d<Numeric::double_epsilon;
+ return d<std::numeric_limits<double>::epsilon();
}
// Calculates if line crosses the box made out of our bins.
diff --git a/Core/SoftParticle/FormFactorSphereUniformRadius.cpp b/Core/SoftParticle/FormFactorSphereUniformRadius.cpp
index 0479eb536d7..19d6fd888b1 100644
--- a/Core/SoftParticle/FormFactorSphereUniformRadius.cpp
+++ b/Core/SoftParticle/FormFactorSphereUniformRadius.cpp
@@ -16,9 +16,9 @@
#include "FormFactorSphereUniformRadius.h"
#include "BornAgainNamespace.h"
#include "Exceptions.h"
-#include "Numeric.h"
#include "Pi.h"
#include "RealParameter.h"
+#include <limits>
FormFactorSphereUniformRadius::FormFactorSphereUniformRadius(double mean,
double full_width)
@@ -40,13 +40,13 @@ complex_t FormFactorSphereUniformRadius::evaluate_for_q(const cvector_t q) const
double W = m_full_width;
double q2 = std::norm(q.x()) + std::norm(q.y()) + std::norm(q.z());
double q_r = std::sqrt(q2);
- if (q_r*R < Numeric::double_epsilon)
+ if (q_r*R < std::numeric_limits<double>::epsilon())
return (4.0*Pi::PI*R*R*R + Pi::PI*R*W*W)/3.0;
double qR = q_r*R;
double qW = q_r*W;
double nominator = 4*Pi::PI*( 4*std::sin(qR)*std::sin(qW/2.0)
- - qW*std::cos(qW/2.0)*std::sin(qR)
- - 2.0*qR*std::cos(qR)*std::sin(qW/2.0) );
+ - qW*std::cos(qW/2.0)*std::sin(qR)
+ - 2.0*qR*std::cos(qR)*std::sin(qW/2.0) );
return nominator/(q2*q2*W);
}
diff --git a/Core/Tools/MathFunctions.cpp b/Core/Tools/MathFunctions.cpp
index a429139923f..71a433b3bb5 100644
--- a/Core/Tools/MathFunctions.cpp
+++ b/Core/Tools/MathFunctions.cpp
@@ -14,7 +14,6 @@
// ************************************************************************** //
#include "MathFunctions.h"
-#include "Numeric.h"
#include "Pi.h"
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_erf.h>
@@ -23,6 +22,7 @@
#include <fftw3.h>
#include <chrono>
#include <cstring>
+#include <limits>
#include <random>
#include <stdexcept> // need overlooked by g++ 5.4
@@ -75,7 +75,7 @@ complex_t MathFunctions::sinc(const complex_t z) // Sin(x)/x
complex_t MathFunctions::tanhc(const complex_t z) // tanh(x)/x
{
- if(std::abs(z)<Numeric::double_epsilon)
+ if(std::abs(z)<std::numeric_limits<double>::epsilon())
return 1.0;
return std::tanh(z)/z;
}
@@ -84,7 +84,7 @@ complex_t MathFunctions::Laue(const complex_t z, size_t N) // Exp(iNx/2)*Sin((N+
{
if (N==0)
return 1.0;
- if(std::abs(z)<Numeric::double_epsilon)
+ if(std::abs(z)<std::numeric_limits<double>::epsilon())
return N+1.0;
return exp_I(N/2.0*z)*std::sin(z*(N+1.0)/2.0)/std::sin(z/2.0);
}
diff --git a/Core/Tools/Numeric.cpp b/Core/Tools/Numeric.cpp
index ec87d422646..ac1db42e1b1 100644
--- a/Core/Tools/Numeric.cpp
+++ b/Core/Tools/Numeric.cpp
@@ -22,22 +22,23 @@
namespace Numeric {
-//! Returns true if two doubles agree within epsilon*tolerance
+//! Returns true if two doubles agree within epsilon*tolerance.
bool areAlmostEqual(double a, double b, double tolerance)
{
constexpr double eps = std::numeric_limits<double>::epsilon();
return std::abs(a-b) <= eps * std::max( tolerance*eps, std::max(1., tolerance)*std::abs(b) );
}
-//! calculates safe relative difference |(a-b)/b|
+//! Returns the safe relative difference, which is |(a-b)/b| except in special cases.
double get_relative_difference(double a, double b)
{
+ constexpr double eps = std::numeric_limits<double>::epsilon();
// return 0.0 if relative error smaller than epsilon
- if (std::abs(a-b) <= double_epsilon*std::abs(b))
+ if (std::abs(a-b) <= eps*std::abs(b))
return 0.0;
// for small numbers, divide by epsilon (to avoid catastrophic cancellation)
- if (std::abs(b) <= double_epsilon)
- return std::abs((a-b)/double_epsilon);
+ if (std::abs(b) <= eps)
+ return std::abs((a-b)/eps);
return std::abs((a-b)/b);
}
diff --git a/Core/Tools/Numeric.h b/Core/Tools/Numeric.h
index f1580ca7cbf..b1f073b03d6 100644
--- a/Core/Tools/Numeric.h
+++ b/Core/Tools/Numeric.h
@@ -23,12 +23,8 @@
namespace Numeric {
-static const double double_epsilon = std::numeric_limits<double>::epsilon();
-
-//! compare two doubles
bool BA_CORE_API_ areAlmostEqual(double a, double b, double tolerance_factor=1.0);
-//! calculates safe relative difference |(a-b)/b|
double BA_CORE_API_ get_relative_difference(double a, double b);
} // Numeric namespace
--
GitLab