diff options
Diffstat (limited to 'zen/basic_math.h')
-rwxr-xr-x[-rw-r--r--] | zen/basic_math.h | 766 |
1 files changed, 383 insertions, 383 deletions
diff --git a/zen/basic_math.h b/zen/basic_math.h index eed23477..7836ae81 100644..100755 --- a/zen/basic_math.h +++ b/zen/basic_math.h @@ -1,383 +1,383 @@ -// ***************************************************************************** -// * This file is part of the FreeFileSync project. It is distributed under * -// * GNU General Public License: http://www.gnu.org/licenses/gpl-3.0 * -// * Copyright (C) Zenju (zenju AT freefilesync DOT org) - All Rights Reserved * -// ***************************************************************************** - -#ifndef BASIC_MATH_H_3472639843265675 -#define BASIC_MATH_H_3472639843265675 - -#include <algorithm> -#include <iterator> -#include <limits> -#include <cmath> -#include <functional> -#include <cassert> - - -namespace numeric -{ -template <class T> T abs(T value); -template <class T> auto dist(T a, T b); -template <class T> int sign(T value); //returns one of {-1, 0, 1} -template <class T> T min(T a, T b, T c); -template <class T> T max(T a, T b, T c); -template <class T> bool isNull(T value); - -template <class T> -void clamp(T& val, T minVal, T maxVal); //make sure minVal <= val && val <= maxVal -template <class T> -T clampCpy(T val, T minVal, T maxVal); - -template <class T, class InputIterator> //precondition: range must be sorted! -auto nearMatch(const T& val, InputIterator first, InputIterator last); - -int round(double d); //"little rounding function" - -template <class N, class D> -auto integerDivideRoundUp(N numerator, D denominator); - -template <size_t N, class T> -T power(T value); - -double radToDeg(double rad); //convert unit [rad] into [°] -double degToRad(double degree); //convert unit [°] into [rad] - -template <class InputIterator> -double arithmeticMean(InputIterator first, InputIterator last); - -template <class RandomAccessIterator> -double median(RandomAccessIterator first, RandomAccessIterator last); //note: invalidates input range! - -template <class InputIterator> -double stdDeviation(InputIterator first, InputIterator last, double* mean = nullptr); //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 <class RandomAccessIterator> -double mad(RandomAccessIterator first, RandomAccessIterator last); //note: invalidates input range! - -template <class InputIterator> -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 <class T> inline -T abs(T value) -{ - //static_assert(std::is_signed<T>::value, ""); - if (value < 0) - return -value; //operator "?:" caveat: may be different type than "value" - else - return value; -} - -template <class T> inline -auto dist(T a, T b) //return type might be different than T, e.g. std::chrono::duration instead of std::chrono::time_point -{ - return a > b ? a - b : b - a; -} - - -template <class T> inline -int sign(T value) //returns one of {-1, 0, 1} -{ - static_assert(std::is_signed<T>::value, ""); - return value < 0 ? -1 : (value > 0 ? 1 : 0); -} - - -template <class T> inline -T min(T a, T b, T c) //don't follow std::min's "const T&(const T&, const T&)" API -{ - if (a < b) - return a < c ? a : c; - else - return b < c ? b : c; - //return std::min(std::min(a, b), c); -} - - -template <class T> inline -T max(T a, T b, T c) -{ - if (a > b) - return a > c ? a : c; - else - return b > c ? b : c; - //return std::max(std::max(a, b), c); -} - - -template <class T> inline -T clampCpy(T val, T minVal, T maxVal) -{ - assert(minVal <= maxVal); - if (val < minVal) - return minVal; - else if (val > maxVal) - return maxVal; - return val; -} - -template <class T> inline -void clamp(T& val, T minVal, T maxVal) -{ - assert(minVal <= maxVal); - if (val < minVal) - val = minVal; - else if (val > maxVal) - val = maxVal; -} - - -/* -part of C++11 now! -template <class InputIterator, class Compare> inline -std::pair<InputIterator, InputIterator> minMaxElement(InputIterator first, InputIterator last, Compare compLess) -{ - //by factor 1.5 to 3 faster than boost::minmax_element (=two-step algorithm) for built-in types! - - InputIterator lowest = first; - InputIterator largest = first; - - if (first != last) - { - auto minVal = *lowest; //nice speedup on 64 bit! - auto maxVal = *largest; // - for (;;) - { - ++first; - if (first == last) - break; - const auto val = *first; - - if (compLess(maxVal, val)) - { - largest = first; - maxVal = val; - } - else if (compLess(val, minVal)) - { - lowest = first; - minVal = val; - } - } - } - return std::make_pair(lowest, largest); -} - - -template <class InputIterator> inline -std::pair<InputIterator, InputIterator> minMaxElement(InputIterator first, InputIterator last) -{ - return minMaxElement(first, last, std::less<typename std::iterator_traits<InputIterator>::value_type>()); -} -*/ - -template <class T, class InputIterator> inline -auto nearMatch(const T& val, InputIterator first, InputIterator last) -{ - if (first == last) - return static_cast<decltype(*first)>(0); - - assert(std::is_sorted(first, last)); - InputIterator it = std::lower_bound(first, last, val); - if (it == last) - return *--last; - if (it == first) - return *first; - - const auto nextVal = *it; - const auto prevVal = *--it; - return val - prevVal < nextVal - val ? prevVal : nextVal; -} - - -template <class T> inline -bool isNull(T value) -{ - return abs(value) <= std::numeric_limits<T>::epsilon(); //epsilon is 0 für integral types => less-equal -} - - -inline -int round(double d) -{ - assert(d - 0.5 >= std::numeric_limits<int>::min() && //if double is larger than what int can represent: - d + 0.5 <= std::numeric_limits<int>::max()); //=> undefined behavior! - return static_cast<int>(d < 0 ? d - 0.5 : d + 0.5); -} - - -template <class N, class D> inline -auto integerDivideRoundUp(N numerator, D denominator) -{ - static_assert(std::is_integral<N>::value && std::is_unsigned<N>::value, ""); - static_assert(std::is_integral<D>::value && std::is_unsigned<D>::value, ""); - assert(denominator > 0); - return (numerator + denominator - 1) / denominator; -} - - -namespace -{ -template <size_t N, class T> struct PowerImpl; -/* - template <size_t N, class T> -> let's use non-recursive specializations to help the compiler - struct PowerImpl { static T result(const T& value) { return PowerImpl<N - 1, T>::result(value) * value; } }; -*/ -template <class T> struct PowerImpl<2, T> { static T result(T value) { return value * value; } }; -template <class T> struct PowerImpl<3, T> { static T result(T value) { return value * value * value; } }; -} - -template <size_t n, class T> inline -T power(T value) -{ - return PowerImpl<n, T>::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 <class InputIterator> inline -double arithmeticMean(InputIterator first, InputIterator last) -{ - size_t n = 0; //avoid random-access requirement for iterator! - double sum_xi = 0; - - for (; first != last; ++first, ++n) - sum_xi += *first; - - return n == 0 ? 0 : sum_xi / n; -} - - -template <class RandomAccessIterator> 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; -} - - -template <class RandomAccessIterator> 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 which may have less than double precision! - - auto lessMedAbs = [m](double lhs, double rhs) { return abs(lhs - m) < abs(rhs - m); }; - - std::nth_element(first, first + n / 2, last, lessMedAbs); //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 - return 0.5 * (abs(*std::max_element(first, first + n / 2, lessMedAbs) - m) + midVal); - } - return 0; -} - - -template <class InputIterator> 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 degradation - - 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 <class InputIterator> 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_H_3472639843265675 +// *****************************************************************************
+// * This file is part of the FreeFileSync project. It is distributed under *
+// * GNU General Public License: http://www.gnu.org/licenses/gpl-3.0 *
+// * Copyright (C) Zenju (zenju AT freefilesync DOT org) - All Rights Reserved *
+// *****************************************************************************
+
+#ifndef BASIC_MATH_H_3472639843265675
+#define BASIC_MATH_H_3472639843265675
+
+#include <algorithm>
+#include <iterator>
+#include <limits>
+#include <cmath>
+#include <functional>
+#include <cassert>
+
+
+namespace numeric
+{
+template <class T> T abs(T value);
+template <class T> auto dist(T a, T b);
+template <class T> int sign(T value); //returns one of {-1, 0, 1}
+template <class T> T min(T a, T b, T c);
+template <class T> T max(T a, T b, T c);
+template <class T> bool isNull(T value);
+
+template <class T>
+void clamp(T& val, T minVal, T maxVal); //make sure minVal <= val && val <= maxVal
+template <class T>
+T clampCpy(T val, T minVal, T maxVal);
+
+template <class T, class InputIterator> //precondition: range must be sorted!
+auto nearMatch(const T& val, InputIterator first, InputIterator last);
+
+int round(double d); //"little rounding function"
+
+template <class N, class D>
+auto integerDivideRoundUp(N numerator, D denominator);
+
+template <size_t N, class T>
+T power(T value);
+
+double radToDeg(double rad); //convert unit [rad] into [°]
+double degToRad(double degree); //convert unit [°] into [rad]
+
+template <class InputIterator>
+double arithmeticMean(InputIterator first, InputIterator last);
+
+template <class RandomAccessIterator>
+double median(RandomAccessIterator first, RandomAccessIterator last); //note: invalidates input range!
+
+template <class InputIterator>
+double stdDeviation(InputIterator first, InputIterator last, double* mean = nullptr); //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 <class RandomAccessIterator>
+double mad(RandomAccessIterator first, RandomAccessIterator last); //note: invalidates input range!
+
+template <class InputIterator>
+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 <class T> inline
+T abs(T value)
+{
+ //static_assert(std::is_signed<T>::value, "");
+ if (value < 0)
+ return -value; //operator "?:" caveat: may be different type than "value"
+ else
+ return value;
+}
+
+template <class T> inline
+auto dist(T a, T b) //return type might be different than T, e.g. std::chrono::duration instead of std::chrono::time_point
+{
+ return a > b ? a - b : b - a;
+}
+
+
+template <class T> inline
+int sign(T value) //returns one of {-1, 0, 1}
+{
+ static_assert(std::is_signed<T>::value, "");
+ return value < 0 ? -1 : (value > 0 ? 1 : 0);
+}
+
+
+template <class T> inline
+T min(T a, T b, T c) //don't follow std::min's "const T&(const T&, const T&)" API
+{
+ if (a < b)
+ return a < c ? a : c;
+ else
+ return b < c ? b : c;
+ //return std::min(std::min(a, b), c);
+}
+
+
+template <class T> inline
+T max(T a, T b, T c)
+{
+ if (a > b)
+ return a > c ? a : c;
+ else
+ return b > c ? b : c;
+ //return std::max(std::max(a, b), c);
+}
+
+
+template <class T> inline
+T clampCpy(T val, T minVal, T maxVal)
+{
+ assert(minVal <= maxVal);
+ if (val < minVal)
+ return minVal;
+ else if (val > maxVal)
+ return maxVal;
+ return val;
+}
+
+template <class T> inline
+void clamp(T& val, T minVal, T maxVal)
+{
+ assert(minVal <= maxVal);
+ if (val < minVal)
+ val = minVal;
+ else if (val > maxVal)
+ val = maxVal;
+}
+
+
+/*
+part of C++11 now!
+template <class InputIterator, class Compare> inline
+std::pair<InputIterator, InputIterator> minMaxElement(InputIterator first, InputIterator last, Compare compLess)
+{
+ //by factor 1.5 to 3 faster than boost::minmax_element (=two-step algorithm) for built-in types!
+
+ InputIterator lowest = first;
+ InputIterator largest = first;
+
+ if (first != last)
+ {
+ auto minVal = *lowest; //nice speedup on 64 bit!
+ auto maxVal = *largest; //
+ for (;;)
+ {
+ ++first;
+ if (first == last)
+ break;
+ const auto val = *first;
+
+ if (compLess(maxVal, val))
+ {
+ largest = first;
+ maxVal = val;
+ }
+ else if (compLess(val, minVal))
+ {
+ lowest = first;
+ minVal = val;
+ }
+ }
+ }
+ return std::make_pair(lowest, largest);
+}
+
+
+template <class InputIterator> inline
+std::pair<InputIterator, InputIterator> minMaxElement(InputIterator first, InputIterator last)
+{
+ return minMaxElement(first, last, std::less<typename std::iterator_traits<InputIterator>::value_type>());
+}
+*/
+
+template <class T, class InputIterator> inline
+auto nearMatch(const T& val, InputIterator first, InputIterator last)
+{
+ if (first == last)
+ return static_cast<decltype(*first)>(0);
+
+ assert(std::is_sorted(first, last));
+ InputIterator it = std::lower_bound(first, last, val);
+ if (it == last)
+ return *--last;
+ if (it == first)
+ return *first;
+
+ const auto nextVal = *it;
+ const auto prevVal = *--it;
+ return val - prevVal < nextVal - val ? prevVal : nextVal;
+}
+
+
+template <class T> inline
+bool isNull(T value)
+{
+ return abs(value) <= std::numeric_limits<T>::epsilon(); //epsilon is 0 für integral types => less-equal
+}
+
+
+inline
+int round(double d)
+{
+ assert(d - 0.5 >= std::numeric_limits<int>::min() && //if double is larger than what int can represent:
+ d + 0.5 <= std::numeric_limits<int>::max()); //=> undefined behavior!
+ return static_cast<int>(d < 0 ? d - 0.5 : d + 0.5);
+}
+
+
+template <class N, class D> inline
+auto integerDivideRoundUp(N numerator, D denominator)
+{
+ static_assert(std::is_integral<N>::value && std::is_unsigned<N>::value, "");
+ static_assert(std::is_integral<D>::value && std::is_unsigned<D>::value, "");
+ assert(denominator > 0);
+ return (numerator + denominator - 1) / denominator;
+}
+
+
+namespace
+{
+template <size_t N, class T> struct PowerImpl;
+/*
+ template <size_t N, class T> -> let's use non-recursive specializations to help the compiler
+ struct PowerImpl { static T result(const T& value) { return PowerImpl<N - 1, T>::result(value) * value; } };
+*/
+template <class T> struct PowerImpl<2, T> { static T result(T value) { return value * value; } };
+template <class T> struct PowerImpl<3, T> { static T result(T value) { return value * value * value; } };
+}
+
+template <size_t n, class T> inline
+T power(T value)
+{
+ return PowerImpl<n, T>::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 <class InputIterator> inline
+double arithmeticMean(InputIterator first, InputIterator last)
+{
+ size_t n = 0; //avoid random-access requirement for iterator!
+ double sum_xi = 0;
+
+ for (; first != last; ++first, ++n)
+ sum_xi += *first;
+
+ return n == 0 ? 0 : sum_xi / n;
+}
+
+
+template <class RandomAccessIterator> 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;
+}
+
+
+template <class RandomAccessIterator> 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 which may have less than double precision!
+
+ auto lessMedAbs = [m](double lhs, double rhs) { return abs(lhs - m) < abs(rhs - m); };
+
+ std::nth_element(first, first + n / 2, last, lessMedAbs); //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
+ return 0.5 * (abs(*std::max_element(first, first + n / 2, lessMedAbs) - m) + midVal);
+ }
+ return 0;
+}
+
+
+template <class InputIterator> inline
+double stdDeviation(InputIterator first, InputIterator last, double* arithMean)
+{
+ //implementation minimizing rounding errors, see: http://en.wikipedia.org/wiki/Standard_deviation
+ //combined with technique avoiding overflow, see: http://www.netlib.org/blas/dnrm2.f -> only 10% performance degradation
+
+ 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 <class InputIterator> 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_H_3472639843265675
|