#ifndef   math_utility_H__
#define   math_utility_H__

#include <cstdlib>
#include <vector>
#include <algorithm>
#include <cmath>

namespace meow {

//! 圓周率...
static const double PI = 3.14159265358979323846264338327950288;

/*!
 * @brief 將角度調整於0~2PI
 */
template<class T>
inline T circle(T x) {
  while (x < 0) x += 2.0 * PI;
  while (2.0 * PI <= x) x -= 2.0 * PI;
  return x;
}

/*!
 * @brief 如果abs(輸入的數值) < eps, 則回傳0, 否則回傳輸入的數值
 */
template<class T>
inline T noEPS(T value, T eps = 1e-9) {
  T epsp((eps < T(0)) ? -eps : eps);
  return ((value < -epsp || value > epsp) ? value : T(0));
}

/*!
 * @brief \c (value-lower)/(upper-lower)
 */
template<class T>
inline T normalize(T lower, T upper, T value) {
  return (value - lower) / (upper - lower);
}

/*!
 * @brief \c (lower+_ratio*(upper-lower))
 */
template<class T>
inline T denormalize(T lower, T upper, T _ratio) {
  return lower + _ratio * (upper - lower);
}

/*!
 * @brief \c denormalize(l2,u2,normalize(l1,u1,m1))
 */
template<class T>
inline T ratioMapping(T l1, T u1, T m1, T l2, T u2) {
  return denormalize(l2, u2, normalize(l1, u1, m1));
}

/*!
 * @brief \c std::min(mx,std::max(mn,v))
 */
template<class T>
inline T inRange(T const& mn, T const& mx, T const& v) {
  return std::min(mx, std::max(mn, v));
}

/*!
 * @brief (mn <= x && x <= mx)
 */
template<class T>
inline T isInRange(T const& mn, T const& mx, T const& x) {
  return (mn <= x && x <= mx);
}

/*!
 * @brief \c x*x
 */
template<class T>
inline T squ(T const& x) {
  return x * x;
}

/*!
 * @brief \c x*x*x
 */
template<class T>
inline T cub(T const& x) {
  return x * x * x;
}

/*!
 * @brief 只將 \c sigs 個標準差以內的數據拿來取平均
 */
template<class T>
inline double average(T const& beg, T const& end, double sigs) {
  int N = 0;
  double av = 0;
  for (T it = beg; it != end; ++it, ++N) {
    av += *it;
  }
  av /= N;
  double sig = 0;
  for (T it = beg; it != end; ++it) {
    sig += (*it - av) * (*it - av);
  }
  sig = sqrt(sig / N);
  double lower = av - sig * sigs, upper = av + sig * sigs;
  double ret = 0, retn = 0;
  for (T it = beg; it != end; ++it) {
    if (lower <= *it && *it <= upper) {
      ret += *it;
      retn++;
    }
  }
  return ret / retn;
}

/*!
 * @brief 只將 \c sigs 個標準差以內的數據拿來取平均, 不過這次用 \c p 來加權平均
 */
template<class T>
inline double average(T const& beg, T const& end, T const& p, double sigs) {
  int N = 0;
  double ps = 0;
  for (T it = beg, ip = p; it != end; ++it, ++N, ++ip) {
    ps += *ip;
  }
  double av = 0;
  for (T it = beg, ip = p; it != end; ++it, ++ip) {
    av += *it * *ip / ps;
  }
  double sig = 0;
  for (T it = beg, ip = p; it != end; ++it, ++ip) {
    sig += *ip / ps * (*it - av) * (*it - av);
  }
  sig = sqrt(sig);
  double lower = av - sig * sigs, upper = av + sig * sigs;
  double ret = 0, retn = 0;
  for (T it = beg, ip = p; it != end; ++it, ++ip) {
    if (lower <= *it && *it <= upper) {
      ret  += *it * *ip;
      retn += *ip;
    }
  }
  if (retn <= 1e-10) return av;
  return ret / retn;
}

/*!
 * @brief 就只是個取絕對值
 */
template<class T>
inline T tAbs(T const& t) {
  return (t < 0 ? -t : t);
}

} // meow

#endif // math_utility_H__