From bd6336c629841c6db3a6ca53a936d629d34db53b Mon Sep 17 00:00:00 2001 From: Daniel Wilhelm Date: Fri, 18 Apr 2014 17:15:16 +0200 Subject: 4.1 --- zen/basic_math.h | 356 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 356 insertions(+) create mode 100644 zen/basic_math.h (limited to 'zen/basic_math.h') diff --git a/zen/basic_math.h b/zen/basic_math.h new file mode 100644 index 00000000..24bcf27a --- /dev/null +++ b/zen/basic_math.h @@ -0,0 +1,356 @@ +// ************************************************************************** +// * This file is part of the FreeFileSync project. It is distributed under * +// * GNU General Public License: http://www.gnu.org/licenses/gpl.html * +// * Copyright (C) 2008-2011 ZenJu (zhnmju123 AT gmx.de) * +// ************************************************************************** + +#ifndef BASIC_MATH_HEADER_34726398432 +#define BASIC_MATH_HEADER_34726398432 + +#include +#include +#include +#include + + +namespace numeric +{ +template +T abs(T value); + +template +T dist(T a, T b); + +template +int sign(T value); //returns -1/0/1 + +template +const T& min(const T& a, const T& b, const T& c); + +template +const T& max(const T& a, const T& b, const T& c); + +template +void confine(T& val, const T& minVal, const T& maxVal); //make sure minVal <= val && val <= maxVal + +template +std::pair minMaxElement(InputIterator first, InputIterator last); +template +std::pair minMaxElement(InputIterator first, InputIterator last, Compare comp); + +template +bool isNull(T value); + +int round(double d); //little rounding function + +template +T power(const T& value); + +double radToDeg(double rad); //convert unit [rad] into [°] +double degToRad(double degree); //convert unit [°] into [rad] + +template +double arithmeticMean(InputIterator first, InputIterator last); + +template +double median(RandomAccessIterator first, RandomAccessIterator last); //note: invalidates input range! + +template +double stdDeviation(InputIterator first, InputIterator last, double* mean = NULL); //estimate standard deviation (and thereby arithmetic mean) + +//median absolute deviation: "mad / 0.6745" is a robust measure for standard deviation of a normal distribution +template +double mad(RandomAccessIterator first, RandomAccessIterator last); //note: invalidates input range! + + +template +double norm2(InputIterator first, InputIterator last); + +//constants +const double pi = 3.14159265358979323846; +const double e = 2.71828182845904523536; +const double sqrt2 = 1.41421356237309504880; +const double ln2 = 0.693147180559945309417; +//---------------------------------------------------------------------------------- + + + + + + + + + + + + + + + + + + + + + + +//################# inline implementation ######################### +template inline +T abs(T value) +{ + return value < 0 ? -1 * value : value; +} + +template inline +T dist(T a, T b) +{ + return a > b ? a - b : b - a; +} + + +template inline +int sign(T value) //returns -1/0/1 +{ + return value < 0 ? -1 : (value > 0 ? 1 : 0); +} + + +template inline +const T& min(const T& a, const T& b, const T& c) +{ + return std::min(std::min(a, b), c); +} + + +template inline +const T& max(const T& a, const T& b, const T& c) +{ + return std::max(std::max(a, b), c); +} + + +template inline +void confine(T& val, const T& minVal, const T& maxVal) +{ + assert(minVal <= maxVal); + if (val < minVal) + val = minVal; + else if (val > maxVal) + val = maxVal; +} + + +template inline +std::pair minMaxElement(InputIterator first, InputIterator last, Compare comp) +{ + InputIterator lowest = first; + InputIterator largest = first; + + if (first != last) + while (++first != last) + { + if (comp(*largest, *first)) // or: if (comp(*largest,*lowest)) for the comp version + largest = first; + else if (comp(*first, *lowest)) + lowest = first; + } + return std::make_pair(lowest, largest); +} + + +template inline +std::pair minMaxElement(InputIterator first, InputIterator last) +{ + return minMaxElement(first, last, std::less::value_type>()); +} + + +template inline +bool isNull(T value) +{ + return abs(value) <= std::numeric_limits::epsilon(); //epsilon is 0 für integral types => less-equal +} + + +inline +int round(double d) +{ + return static_cast(d < 0 ? d - 0.5 : d + 0.5); +} + + +namespace +{ +template +struct PowerImpl +{ + static T result(const T& value) + { + return PowerImpl::result(value) * value; + } +}; + +template +struct PowerImpl<2, T> +{ + static T result(const T& value) + { + return value * value; + } +}; + +template +struct PowerImpl<0, T>; //not defined: invalidates power<0> and power<1> + +template +struct PowerImpl<10, T>; //not defined: invalidates power for N >= 10 +} + +template inline +T power(const T& value) +{ + return PowerImpl::result(value); +} + + +inline +double radToDeg(double rad) +{ + return rad * 180.0 / numeric::pi; +} + + +inline +double degToRad(double degree) +{ + return degree * numeric::pi / 180.0; +} + + +template inline +double arithmeticMean(InputIterator first, InputIterator last) +{ + //low level implementation to avoid random-access requirement on iterator + size_t n = 0; + double sum_xi = 0; + + for (; first != last; ++first, ++n) + sum_xi += *first; + + return n == 0 ? 0 : sum_xi / n; +} + + +template inline +double median(RandomAccessIterator first, RandomAccessIterator last) //note: invalidates input range! +{ + const size_t n = last - first; + if (n > 0) + { + std::nth_element(first, first + n / 2, last); //complexity: O(n) + const double midVal = *(first + n / 2); + + if (n % 2 != 0) + return midVal; + else //n is even and >= 2 in this context: return mean of two middle values + return 0.5 * (*std::max_element(first, first + n / 2) + midVal); //this operation is the reason why median() CANNOT support a comparison predicate!!! + } + return 0; +} + + +class LessMinusMedAbs : public std::binary_function +{ +public: + LessMinusMedAbs(double median) : median_(median) {} + bool operator()(double lhs, double rhs) const + { + return abs(lhs - median_) < abs(rhs - median_); + } +private: + double median_; +}; + + +template inline +double mad(RandomAccessIterator first, RandomAccessIterator last) //note: invalidates input range! +{ + //http://en.wikipedia.org/wiki/Median_absolute_deviation + + const size_t n = last - first; + if (n > 0) + { + const double m = median(first, last); + + //the second median needs to operate on absolute residuals => avoid transforming input range as it may decrease precision! + + std::nth_element(first, first + n / 2, last, LessMinusMedAbs(m)); //complexity: O(n) + const double midVal = abs(*(first + n / 2) - m); + + if (n % 2 != 0) + return midVal; + else //n is even and >= 2 in this context: return mean of two middle values + { + const double midVal2 = abs(*std::max_element(first, first + n / 2, LessMinusMedAbs(m)) - m); + return 0.5 * (midVal2 + midVal); + } + } + return 0; +} + + +template inline +double stdDeviation(InputIterator first, InputIterator last, double* arithMean) +{ + //implementation minimizing rounding errors, see: http://en.wikipedia.org/wiki/Standard_deviation + //combined with techinque avoiding overflow, see: http://www.netlib.org/blas/dnrm2.f -> only 10% performance degregation + + size_t n = 0; + double mean = 0; + double q = 0; + double scale = 1; + + for (; first != last; ++first) + { + ++n; + const double val = *first - mean; + + if (abs(val) > scale) + { + q = (n - 1.0) / n + q * power<2>(scale / val); + scale = abs(val); + } + else + q += (n - 1.0) * power<2>(val / scale) / n; + + mean += val / n; + } + + if (arithMean) + *arithMean = mean; + + return n <= 1 ? 0 : std::sqrt(q / (n - 1)) * scale; +} + + +template inline +double norm2(InputIterator first, InputIterator last) +{ + double result = 0; + double scale = 1; + for (; first != last; ++first) + { + const double tmp = abs(*first); + if (tmp > scale) + { + result = 1 + result * power<2>(scale / tmp); + scale = tmp; + } + else + result += power<2>(tmp / scale); + } + return std::sqrt(result) * scale; +} +} + +#endif //BASIC_MATH_HEADER_34726398432 -- cgit