rmath.cc

Go to the documentation of this file.

#include "config.h"

#include "error.h"
#include "rmath.h"

/** A small set of numeric limits routines, since gcc prior to 3.x doesn't
 * have numeric_limits */

// ---------------------------------------------------------------------------

/** Find the maximum limit for a type, equivalent to
 * std::numeric_limits<T>.max() for systems that don't have the limits c++
 * header file. */
// template<typename T>
// static const T max_limit()
// {
//      T tmp = 0;
//      static T max = 0;
//      static bool found = false;
// 
//      if (found) {
//              return(max);
//      }
//      
//      tmp = 1;
//      max = 0;
//      while (tmp > max) {
//              max = tmp;
//              try {
//                      tmp = (tmp * T(2)) + T(1);
//              }
//              catch(...) { }
//      }
//      found = true;
//      return(max);
// }

/** Find the maximum limit for a type, equivalent to
 * std::numeric_limits<T>.min() for systems that don't have the limits c++
 * header file.
 */
// template<typename T>
// static const T min_limit()
// {
//      T tmp = 0;
//      static T min = 0;
//      static bool found = false;
// 
//      if (found) {
//              return(min);
//      }
//      
//      tmp = -1;
//      min = 0;
//      while (tmp < min) {
//              min = tmp;
//              try {
//                      tmp = (tmp * T(2)) - T(1);
//              }
//              catch(...) { }
//      }
//      
//      try {
//              tmp = - max_limit<T>() - T(1);
//      }
//      catch(...) { }
//      if (tmp < min) {
//              min = tmp;
//      }
// 
//      found = true;
//      return(min);
// }

/** Return the largest possible number that a float may hold. */
#ifdef OSX_WORKAROUND
template<>
const float max_limit<float>()
{
        return(FLT_MAX);
}
#endif

/** Return the smallest positive number that a float may hold.
 *
 * Caveat: This is in contrast to other types, where min_limit<T>() will return
 * either 0 or the largest possible negative number that the type may hold.  If
 * you are looking for the largest possible negative number for any given type,
 * use lowest_value<T>() instead.
 */
#ifdef OSX_WORKAROUND
template<>
const float min_limit<float>()
{
        return(FLT_MIN);
}
#endif

/** Return the largest possible number that a double may hold. */
#ifdef OSX_WORKAROUND
template<>
const double max_limit<double>()
{
        return(DBL_MAX);
}
#endif

/** Return the smallest positive number that a double may hold.
 *
 * Caveat: This is in contrast to other types, where min_limit<T>() will return
 * either 0 or the largest possible negative number that the type may hold.  If
 * you are looking for the largest possible negative number for any given type,
 * use lowest_value<T>() instead.
 */
#ifdef OSX_WORKAROUND
template<>
const double min_limit<double>()
{
        return(DBL_MIN);
}
#endif

// ---------------------------------------------------------------------------

/** Return the max_limit of a variable.
 *
 * This is handy to have because it means
 * that the author may change the type of a variable without having to track
 * down all uses of max_limit or min_limit to change the type they measure.
 */
// template<typename T>
// static const T max_limit(const T& a_arg)
// {
//      T value;
// 
//      value = max_limit<T>();
// 
//      return(value);
// }

/** Return the min_limit of a variable.
 *
 * This is handy to have because it means
 * that the author may change the type of a variable without having to track
 * down all uses of max_limit or min_limit to change the type they measure.
 */
// template<typename T>
// static const T min_limit(const T& a_arg)
// {
//      T value;
// 
//      value = min_limit<T>();
// 
//      return(value);
// }

// ---------------------------------------------------------------------------

/** Return the maximum possible value a type may hold.
 *
 * This is just a convenience function to match lowest_value<T>().  All it does
 * is return the value of max_limit<T>().
 */
// template<typename T>
// static const T highest_value(void)
// {
//      T value;
// 
//      value = max_limit<T>();
// 
//      return(value);
// }

/** Return the maximum possible value of a variable. */
// template<typename T>
// static const T highest_value(const T& a_arg)
// {
//      T value;
// 
//      value = highest_value<T>();
// 
//      return(value);
// }

/** Return 0 for unsigned types, or the maximum negative value that the type
 * may hold.
 */
// template<typename T>
// static const T lowest_value(void)
// {
//      T value;
// 
//      value = -max_limit<T>();
//      if (value > min_limit<T>())
//              value = min_limit<T>();
//      
//      return(value);
// }

/** Return 0 for unsigned types, or the maximum negative value that a variable
 * may hold.
 */
// template<typename T>
// static const T lowest_value(const T& a_arg)
// {
//      T value;
// 
//      value = lowest_value<T>();
// 
//      return(value);
// }

// ---------------------------------------------------------------------------

/** Return the absolute value of a numeric type 
 *
 * Caveat: For some types, the maximum negative value is one larger than the
 * maximum positive value: specifically char.  Depending on the type and value,
 * it may be impossible to return the absolute value.  For such types under
 * such circumstances an exception is thrown.
 */
// template<typename T>
// T absolute(const T& a_num)
// {
//      T num;
//      std::string es;
// 
// // std::cerr << "absolute(" << static_cast<long long>(a_num) << ") : BEGIN" << std::endl;
//      num = a_num;
//      if ((num < 0) && (num == min_limit<T>())) {
//              TRY_nomem(es = "Absolute value is not containable by this type: ");
//              if (is_char(num)) {
//                      TRY_nomem(es += estring(static_cast<int>(num)));
//              }
//              else {
//                      TRY_nomem(es += estring(num));
//              }
//              throw(ERROR(0,es));
//      }
//      if (num < 0) {
//              num = -num;
//              if (num < 0) {
//                      TRY_nomem(es = "Absolute value is still negative: ");
//                      if (is_char(num)) {
//                              TRY_nomem(es += estring(static_cast<int>(num)));
//                      }
//                      else {
//                              TRY_nomem(es += estring(num));
//                      }
//                      throw(INTERNAL_ERROR(0,es));
//              }
// // else std::cerr << "absolute(" << static_cast<long long>(a_num) << ") >= 0" << std::endl;
//      }
// 
// // std::cerr << "absolute(" << static_cast<long long>(a_num) << ") = " << static_cast<long long>(num) << std::endl;
// // std::cerr << "absolute(" << static_cast<long long>(a_num) << ") : END" << std::endl;
//      return(num);
// }

// ---------------------------------------------------------------------------

/** Safely manipulate numbers without worryiung about over/underflow error */
// template<typename T>
// class safe_num
// {
// public:
//      /** C'tor */
//      safe_num()
//      {
//              clear();
//      }
// 
//      /** C'tor */
//      safe_num(const T a_num)
//      {
//              clear();
//              m_num = a_num;
//      }
// 
//      /** C'tor */
//      safe_num(const safe_num& a_class)
//      {
//              clear();
//              m_num = a_class.value();
//      }
// 
//      /** Clear the value */
//      void clear(void)
//      {
//              m_num = static_cast<T>(0);
//      }
// 
//      /** Return the value */
//      const T value(void) const
//      {
//              return(m_num);
//      }
// 
//      /** Assign a value */
//      void assign(const T& a_arg)
//      {
//              m_num = a_arg;
//      }
// 
//      /** Add a value */
//      void add(const T& a_arg) 
//      {
//              bool overflow = false;
//              T num;
// 
//              if (a_arg < static_cast<T>(0)) {
//                      num = absolute(a_arg);
//                      subtract(num);
//                      return;
//              }
// 
//              if (highest_value<T>() - a_arg < m_num)
//                      overflow = true;
// 
//              if (overflow) {
//                      estring es;
// 
//                      es = "Addition overflow error detected: ";
//                      if (is_char(m_num))
//                              es += estring(static_cast<int>(m_num));
//                      else
//                              es += estring(m_num);
//                      es += " + ";
//                      if (is_char(a_arg))
//                              es += estring(static_cast<int>(a_arg));
//                      else
//                              es += estring(a_arg);
//                      throw(INTERNAL_ERROR(0,es));
//              }
//              m_num += a_arg;
//      }
// 
//      /** Subtract a value */
//      void subtract(const T& a_arg) 
//      {
//              bool underflow = false;
//              T num;
// 
//              if (a_arg < static_cast<T>(0)) {
//                      num = absolute(a_arg);
//                      add(num);
//                      return;
//              }
// 
//              if (lowest_value<T>() < 0) {
//                      if (m_num < lowest_value<T>() + a_arg)
//                              underflow = true;
//              }
//              else {
//                      if (m_num - lowest_value<T>() < a_arg)
//                              underflow = true;
//              }
// 
//              if (underflow) {
//                      estring es;
// 
//                      es = "Subtraction underflow error detected: ";
//                      if (is_char(m_num))
//                              es += estring(static_cast<int>(m_num));
//                      else
//                              es += estring(m_num);
//                      es += " - ";
//                      if (is_char(a_arg))
//                              es += estring(static_cast<int>(a_arg));
//                      else
//                              es += estring(a_arg);
//                      throw(INTERNAL_ERROR(0,es));
//              }
//              m_num -= a_arg;
//      }
// 
//      /** Multiply by a value */
//      void multiply(const T& a_arg) 
//      {
//              bool overflow = false;
// 
//              if ((a_arg == 0) || (m_num == 0) || (a_arg == 1) || (m_num == 1)) {
//                      m_num *= a_arg;
//                      return;
//              }
// 
//              if ((lowest_value<T>() < 0) && (m_num < 0) && (a_arg < 0)) {
//                      if (-highest_value<T>() > a_arg)
//                              overflow = true;
//                      if (-(highest_value<T>() / a_arg) < -m_num)
//                              overflow = true;
//              }
//              else
//              if ((lowest_value<T>() < 0) && (m_num < 0) && (a_arg >= 0)) {
//                      if (lowest_value<T>() / a_arg > m_num)
//                              overflow = true;
//              }
//              else
//              if ((lowest_value<T>() < 0) && (m_num >= 0) && (a_arg < 0)) {
//                      if (lowest_value<T>() / m_num > a_arg)
//                              overflow = true;
//              }
//              else
//              if ((lowest_value<T>() < 0) && (m_num >= 0) && (a_arg >= 0)) {
//                      if (highest_value<T>() / a_arg < m_num)
//                              overflow = true;
//              }
//              else
//              if ((lowest_value<T>() >= 0) && (m_num >= 0) && (a_arg >= 0)) {
//                      if (highest_value<T>() / a_arg < m_num)
//                              overflow = true;
//              }
//              else {
//                      // This should never happen
//                      ASSERT(0);
//              }
// 
//              if (overflow) {
//                      estring es;
// 
//                      es = "Multiplication overflow error detected: ";
//                      if (is_char(m_num))
//                              es += estring(static_cast<int>(m_num));
//                      else
//                              es += estring(m_num);
//                      es += " * ";
//                      if (is_char(a_arg))
//                              es += estring(static_cast<int>(a_arg));
//                      else
//                              es += estring(a_arg);
//                      throw(INTERNAL_ERROR(0,es));
//              }
//              m_num *= a_arg;
//      }
// 
//      /** Divide by a value */
//      void divide(const T& a_arg) 
//      {
//              T num;
// 
//              if (a_arg == static_cast<T>(0)) {
//                      estring es;
// 
//                      es = "Division by zero error detected: ";
//                      if (is_char(m_num))
//                              es += estring(static_cast<int>(m_num));
//                      else
//                              es += estring(m_num);
//                      es += " / ";
//                      if (is_char(a_arg))
//                              es += estring(static_cast<int>(a_arg));
//                      else
//                              es += estring(a_arg);
//                      throw(INTERNAL_ERROR(0,es));
//              }
// 
//              num = m_num;
//              num /= a_arg;
//              if ((m_num < 0) && (a_arg < 0) && (num < 0)) {
//                      estring es;
// 
//                      es = "Division result has incorrect sign: ";
//                      if (is_char(m_num))
//                              es += estring(static_cast<int>(m_num));
//                      else
//                              es += estring(m_num);
//                      es += " / ";
//                      if (is_char(a_arg))
//                              es += estring(static_cast<int>(a_arg));
//                      else
//                              es += estring(a_arg);
//                      es += " != ";
//                      if (is_char(num))
//                              es += estring(static_cast<int>(num));
//                      else
//                              es += estring(num);
//                      throw(INTERNAL_ERROR(0,es));
//              }
// 
//              m_num = num;
//              return;
//      }
// 
//      /** Assign a safe_num */
//      void assign(const safe_num<T>& a_class) 
//      {
//              assign(a_class.value());
//      }
// 
//      /** Add a safe_num */
//      void add(const safe_num<T>& a_class) 
//      {
//              add(a_class.value());
//      }
// 
//      /** Subtract a safe_num */
//      void subtract(const safe_num<T>& a_class)
//      {
//              subtract(a_class.value());
//      }
// 
//      /** Multiply by a safe_num */
//      void multiply(const safe_num<T>& a_class) 
//      {
//              multiply(a_class.value());
//      }
// 
//      /** Divide by a safe_num */
//      void divide(const safe_num<T>& a_class)
//      {
//              divide(a_class.value());
//      }
// 
//      /** Boolean operator */
//      const bool operator==(const T& a_arg) const 
//      {
//              bool value;
// 
//              value = (m_num == a_arg);
// 
//              return(value);
//      }
// 
//      /** Boolean operator */
//      const bool operator==(const safe_num<T>& a_class) const 
//      {
//              bool value;
// 
//              value = (m_num == a_class.value());
// 
//              return(value);
//      }
// 
//      /** Boolean operator */
//      const bool operator!=(const T& a_arg) const
//      {
//              bool value;
// 
//              value = (m_num != a_arg);
// 
//              return(value);
//      }
// 
//      /** Boolean operator */
//      const bool operator!=(const safe_num<T>& a_class) const
//      {
//              bool value;
// 
//              value = (m_num != a_class.value());
// 
//              return(value);
//      }
// 
//      /** Boolean operator */
//      const bool operator<(const T& a_arg) const 
//      {
//              bool value;
// 
//              value = (m_num < a_arg);
// 
//              return(value);
//      }
// 
//      /** Boolean operator */
//      const bool operator<(const safe_num<T>& a_class) const
//      {
//              bool value;
// 
//              value = (m_num < a_class.value());
// 
//              return(value);
//      }
// 
//      /** Boolean operator */
//      const bool operator>(const T& a_arg) const
//      {
//              bool value;
// 
//              value = (m_num > a_arg);
// 
//              return(value);
//      }
// 
//      /** Boolean operator */
//      const bool operator>(const safe_num<T>& a_class) const 
//      {
//              bool value;
// 
//              value = (m_num > a_class.value());
// 
//              return(value);
//      }
// 
//      /** Boolean operator */
//      const bool operator<=(const T& a_arg) const
//      {
//              bool value;
// 
//              value = (m_num <= a_arg);
// 
//              return(value);
//      }
// 
//      /** Boolean operator */
//      const bool operator<=(const safe_num<T>& a_class) const
//      {
//              bool value;
// 
//              value = (m_num <= a_class.value());
// 
//              return(value);
//      }
// 
//      /** Boolean operator */
//      const bool operator>=(const T& a_arg) const
//      {
//              bool value;
// 
//              value = (m_num >= a_arg);
// 
//              return(value);
//      }
// 
//      /** Boolean operator */
//      const bool operator>=(const safe_num<T>& a_class) const
//      {
//              bool value;
// 
//              value = (m_num >= a_class.value());
// 
//              return(value);
//      }
// 
//      /** Arithmetic operator */
//      safe_num<T>& operator+=(safe_num<T> a_class)
//      {
//              add(a_class);
// 
//              return(*this);
//      }
// 
//      /** Arithmetic operator */
//      safe_num<T>& operator-=(safe_num<T> a_class)
//      {
//              subtract(a_class);
// 
//              return(*this);
//      }
// 
//      /** Arithmetic operator */
//      safe_num<T>& operator*=(safe_num<T> a_class)
//      {
//              multiply(a_class);
// 
//              return(*this);
//      }
// 
//      /** Arithmetic operator */
//      safe_num<T>& operator/=(safe_num<T> a_class)
//      {
//              divide(a_class);
// 
//              return(*this);
//      }
// 
//      /** Arithmetic operator */
//      safe_num<T>& operator%=(safe_num<T> a_class)
//      {
//              m_num %= a_class.value();
// 
//              return(*this);
//      }
// 
// private:
//      T m_num;
// };

/** Arithmetic operator */
// template<typename T>
// safe_num<T> operator+(safe_num<T> a_class1, safe_num<T> a_class2)
// {
//      safe_num<T> result;
// 
//      result.assign(a_class1);
//      result.add(a_class2);
// 
//      return(result);
// }

/** Arithmetic operator */
// template<typename T>
// safe_num<T> operator-(safe_num<T> a_class1, safe_num<T> a_class2)
// {
//      safe_num<T> result;
// 
//      result.assign(a_class1);
//      result.subtract(a_class2);
//      
//      return(result);
// }

/** Arithmetic operator */
// template<typename T>
// safe_num<T> operator*(safe_num<T> a_class1, safe_num<T> a_class2)
// {
//      safe_num<T> result;
// 
//      result.assign(a_class1);
//      result.multiply(a_class2);
// 
//      return(result);
// }

/** Arithmetic operator */
// template<typename T>
// safe_num<T> operator/(safe_num<T> a_class1, safe_num<T> a_class2)
// {
//      safe_num<T> result;
// 
//      result.assign(a_class1);
//      result.divide(a_class2);
// 
//      return(result);
// }

/** Arithmetic operator */
// template<typename T>
// safe_num<T> operator%(safe_num<T> a_class1, safe_num<T> a_class2)
// {
//      safe_num<T> result;
// 
//      result.assign(a_class1.value() % a_class2.value());
// 
//      return(result);
// }

/** Arithmetic operator */
// template<typename T>
// safe_num<T>& operator++(safe_num<T>& a_class)
// {
//      a_class.add(static_cast<T>(1));
// 
//      return(a_class);
// }

/** Arithmetic operator */
// template<typename T>
// safe_num<T> operator++(safe_num<T>& a_class, int)
// {
//      safe_num<T> result;
// 
//      result.assign(a_class);
//      a_class.add(static_cast<T>(1));
// 
//      return(result);
// }

/** Arithmetic operator */
// template<typename T>
// safe_num<T>& operator--(safe_num<T>& a_class)
// {
//      a_class.subtract(static_cast<T>(1));
// 
//      return(a_class);
// }

/** Arithmetic operator */
// template<typename T>
// safe_num<T> operator--(safe_num<T>& a_class, int)
// {
//      safe_num<T> result;
// 
//      result.assign(a_class);
//      a_class.subtract(static_cast<T>(1));
// 
//      return(result);
// }

/** Arithmetic operator */
// template<typename T>
// std::ostream& operator<<(std::ostream& a_out, safe_num<T> a_class)
// {
//      a_out << a_class.value();
// 
//      return(a_out);
// }

/** Arithmetic operator */
// template<typename T>
// std::istream& operator>>(std::istream& a_in, safe_num<T>& a_class)
// {
//      T num;
// 
//      a_in >> num;
//      a_class = num;
// 
//      return(a_in);
// }

// ---------------------------------------------------------------------------