---

more/math/math.h

Go to the documentation of this file.

//  more/math.h -- various mathematical functions
//  Copyright (C) 1998--2002  Petter Urkedal
//
//  This file is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.
//
//  This file is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
//  As a special exception, you may use this file as part of a free
//  software library without restriction.  Specifically, if other files
//  instantiate templates or use macros or inline functions from this
//  file, or you compile this file and link it with other files to
//  produce an executable, this file does not by itself cause the
//  resulting executable to be covered by the GNU General Public
//  License.  This exception does not however invalidate any other
//  reasons why the executable file might be covered by the GNU General
//  Public License.
//
//  $Id: math.h,v 1.2 2002/08/24 13:53:19 petter_urkedal Exp $

//  This file is primarily meant for signatures of float, double and
//  long double.  But some generic templates may be defined here.

#ifndef MORE_MATH_MATH_H
#define MORE_MATH_MATH_H

#include <stdexcept>
#include <string>
#include <more/bits/config.h>
#include <more/bits/mathbits.h>
#include <more/bits/infnan.h>
#include <cmath>
#include <cstdlib>
#include <cfloat>

namespace more {
namespace math {

  extern "C" {
    double besjr(int, double);
    double besIr(double, double);
    double besir(int, double);
    double lagLc(int, int);
  }

  using std::ceil;  using std::floor;
  using std::sin;   using std::cos;   using std::tan;
  using std::asin;  using std::acos;  using std::atan;
  using std::sinh;  using std::cosh;  using std::tanh;
  using std::exp;   using std::log;   using std::log10;
  using std::sqrt;
  using std::fabs;  using std::abs;
#ifdef CXX_HAVE_STD_DRAND48
  using std::drand48; using std::srand48;
#elif CXX_HAVE_DRAND48
  using ::drand48; using ::srand48;
#endif

  using std::atan2;
  using std::fmod;
  using std::ldexp;
  using std::frexp;
  using std::pow;


  /** Numerical constants.
   **
   ** The names of the constants are created in the following way.
   **
   ** Operators and numbers are separated by spaces.  Operators
   ** applies to the following item(s).  Binary operators are
   ** <ul>
   **     <li>\c f — fraction
   **     <li>\c p — product
   **     <li>\c x — exponentiation; \c x_a_b means <tt>pow(a, b)</tt>
   ** </ul>
   ** To express two binary operators applied to three arguments we
   ** use \c pp_x_y_z instead of \c p_p_x_y_z or \c p_x_p_y_z, etc.
   **
   ** If unary operators, like \c log is followed by a number (before
   ** the underscore) it belongs to the operator, and in the case of
   ** \c log, specifies the base.  Unary operators are \c sqrt (square
   ** root), \c log, \c log2, \c log10 (logarithms)
   **/
  namespace numbers {
    const double e =            2.7182818284590452354;
    const double log2_e =       1.4426950408889634074;  // log_2(e)
    const double log2_10 =      3.32192809488736234789; // log2(10)
    const double log_2 =        0.69314718055994530942; // log_e(2)
    const double log_10 =       2.30258509299404568402; // log_e(10)
    const double log10_2 =      0.30102999566398119521; // log10(2)
    const double log10_e =      0.43429448190325182765; // log_10(e)
    const double log_p_2_pi =   1.837877066409345;
    const double pi =           3.14159265358979323846; // pi
    const double p_4_pi =       12.56637061435917;      // 4*pi
    const double p_2_pi =       6.28318530717958647692; // 2*pi
    const double sqrt_pi =      1.77245385090551602729; // sqrt(pi)
    const double f_pi_2 =       1.57079632679489661923; // pi/2
    const double f_pi_4 =       0.78539816339744830962; // pi/4
    const double f_1_pi =       0.31830988618379067154; // 1/pi
    const double f_2_pi =       0.63661977236758134308; // 2/pi
    const double f_2_sqrt_pi =  1.12837916709551257390; // 2/sqrt(pi)
    const double sqrt_2 =       1.41421356237309504880; // sqrt(2)
    const double f_1_sqrt_2 =   0.70710678118654752440; // 1/sqrt(2)
    const double f_1_6 =        0.16666666666666666667;
    const double f_1_3 =        0.33333333333333333333;
    const double eulC =         0.57721566490153286061;  // Euler's constant
    const double catG =         0.915965594;  // Catalan's constant
  }


  // --- some additional functions and signatures ---

# ifdef CXX_HAVE_STD_FINITE
  using std::finite;
# elif defined(CXX_HAVE_FINITE)
  using ::finite;
# else
#   if defined (CXX_HAVE_STD_ISNAN) && defined(CXX_HAVE_STD_ISINF)
  using std::isnan; using std::isinf;
  inline int finite(double x) { return !isinf(x) && !isnan(x); }
#   elif defined(CXX_HAVE_ISNAN_ISINF)
  inline int finite(double x) { return !isinf(x) && !isnan(x); }
#   else
  inline int finite(double x) { return 1; }
#   endif
# endif

  template<typename T>
    inline T max(T const& x, T const& y) { return x < y? y : x; }
  template<typename T>
    inline T min(T const& x, T const& y) { return x < y? x : y; }
  template<typename T>
    inline T abs(T const& x) { return x >= T(0)? x : -x; }

  // ISO/IEC 14882:1998(E) 4.9.1.  An rvalue of a floating point type
  // can be converted to an rvalue of an integer type. The conversion
  // truncates; that is, the fractional part is discarded. The
  // behavior is undefined if the truncated value cannot be
  // represented in the destination type. [Note: If the destination
  // type is bool, see 4.12. ]
//   inline int ifloor(double x)
//      { return x >= 0.0? (int)x : (int)(x-x*DBL_EPSILON)-1; }
//   inline int iceil(double x)
//      { return x > 0.0? (int)(x-x*DBL_EPSILON)+1 : (int)x; }
//   inline int itrunc(double x) { return (int)x; }
//   inline int iround(double x) { return ifloor(x+0.5); }

#if HAVE_STD_LROUND
  using std::lround(double);
#elif HAVE_LROUND
  using ::lround(double);
  inline int iround(double x) { return lround(x); }
#else
  inline long lround(double x) { return long(x >= 0.0? x + 0.5 : x - 0.5); }
  inline int iround(double x) { return int(x >= 0.0? x + 0.5 : x - 0.5); }
#endif
  inline long lfloor(double x) { return long(floor(x)); }
  inline long lceil(double x)  { return long(ceil(x)); }
  inline long ltrunc(double x) { return long(x); }

  inline int ifloor(double x) { return int(floor(x)); }
  inline int iceil(double x) { return int(ceil(x)); }
  inline int itrunc(double x) { return int(x); }

  inline float real(float x) { return x; }
  inline double real(double x) { return x; }
  inline long double real(long double x) { return x; }
  inline float imag(float x) { return 0.0; }
  inline double imag(double x) { return 0.0; }
  inline long double imag(long double x) { return 0.0L; }

  template<typename T>
    inline T pow2(const T& x) { return x*x; }
  template<typename T>
    inline T pow3(const T& x) { return x*x*x; }
  template<typename T>
    inline T pow4(const T& x) { return pow2(pow2(x)); }
  template<typename T>
    inline T pow5(const T& x) { return x*pow4(x); }
  template<typename T>
    inline T pow6(const T& x) { return pow2(pow3(x)); }
  template<typename T>
    inline T pow7(const T& x) { return x*pow6(x); }
  template<typename T>
    inline T pow8(const T& x) { return pow2(pow4(x)); }
  template<typename T>
    inline T pow9(const T& x) { return pow8(x)*x; }
  template<typename T>
    inline T pow10(const T& x) { return pow2(pow5(x)); }

  inline float conj(float x) { return x; }
  inline double conj(double x) { return x; }
  inline long double conj(long double x) { return x; }

  inline unsigned int norm(unsigned int x) { return x*x; }
  inline int norm(int x) { return x*x; }
  inline float norm(float x) { return x*x; }
  inline double norm(double x) { return x*x; }
  inline long double norm(long double x) { return x*x; }

  inline int ldexp(int x, int n) { return x << n; }
  inline int ldexp(long x, int n) { return x << n; }
  inline unsigned int ldexp(unsigned int x, int n) { return x << n; }
  inline unsigned long ldexp(unsigned long x, int n) { return x << n; }

  inline int rank(float) { return 0; }
  inline int rank(double) { return 0; }
  inline int rank(long double) { return 0; }
  inline int size0(float) { return 1; }
  inline int size0(double) { return 1; }
  inline int size0(long double) { return 1; }
  inline int size1(float) { return 1; }
  inline int size1(double) { return 1; }
  inline int size1(long double) { return 1; }
  inline int size2(float) { return 1; }
  inline int size2(double) { return 1; }
  inline int size2(long double) { return 1; }
  template<typename T> inline int rank(T const&x) { return x.rank(); }
  template<typename T> inline int size0(T const&x) { return x.size0(); }
  template<typename T> inline int size1(T const&x) { return x.size1(); }
  template<typename T> inline int size2(T const&x) { return x.size2(); }
  template<typename T> inline T trace(T const& x) { return x; }

  inline float dot(float x, float y) { return x*y; }
  inline double dot(double x, double y) { return x*y; }
  inline long double dot(long double x, long double y) { return x*y; }

  template<typename T>
    inline T
    subscript(T const& x, int i, int j=0, int k=0)
    {
#ifdef MORE_CHECK_BOUNDS
        if (i != 0 || j != 0 || k != 0)
            throw std::logic_error("more::math::subscript: "
                                   "Index out of bounds in subscript().");
#endif
        return x;
    }

  template<typename T>
    inline T&
    subscript(T& x, int i, int j=0, int k=0)
    {
#ifdef MORE_CHECK_BOUNDS
        if (i != 0 || j != 0 || k != 0)
            throw std::logic_error("more::math::subscript: "
                                   "Index out of bounds in subscript().");
#endif
        return x;
    }


  // --- combinatorics ---

  //  XXX optimize
  inline int
  factorial(int n, int m = 0)
  {
#ifdef MORE_THROW_LOGIC_ERROR
      if (n < m)
          throw std::domain_error("more::math::factorial(int n, int m): "
                                  "Undefined for n < m.");
#endif
      if(n == m) return 1;
      int k = n;
      while(--n > m) k *= n;
      return k;
  }

  //  XXX optimize
  inline int
  double_factorial(int n, int m = 0)
  {
#ifdef MORE_THROW_LOGIC_ERROR
      if(n < m)
          throw std::domain_error("more::math::double_factorial"
                                  "(int n, int m): "
                                  "Undefined for n < m.");
      else if (n%2 != m%2)
          throw std::domain_error("more::math::double_factorial"
                                  "(int n, int m): "
                                  "Only defined for n, m both even"
                                  " or both odd.");
#endif
      int k = 1;
      while(n > m) { k *= n; n -= 2; }
      return k;
  }

  //  XXX optimize
  inline int
  binominal(int n, int i)
  {
#ifdef MORE_THROW_LOGIC_ERROR
      if (n < 0)
          throw std::domain_error("more::math::binominal(int n, int): "
                                  "Undefined for n < 0.");
#endif
      if (i >= n)
          return i == n? 1 : 0;
      if (i <= 0)
          return i == 0? 1 : 0;
      int j = n - i;
      if (i > j) {
          j = i;
          i = n - j;
      }
      return factorial(n, j)/factorial(i);
  }


  // --- miscellaneous ---

//   template<typename T>
//     inline T partial_exp(int n, T x) {
//      T res = 1.0;
//      T term = 1.0;
//      for(int i = 1; i <= n; ++i) {
//          term *= x/i;
//          res += term;
//      }
//      return res;
//     }

  // --- a fractal function ---

  unsigned int bitpow(unsigned int, int);
  unsigned long bitpow(unsigned long, int);
#ifdef CXX_HAVE_LONG_LONG
  unsigned long long bitpow(unsigned long long, int);
#endif
  float bitpow(float, int);
  double bitpow(double, int);
  long double bitpow(long double, int);

  int upper_bit(unsigned int);
  int lower_bit(unsigned int);
//   int lowbit(unsigned long);
// #ifdef CXX_HAVE_LONG_LONG
//   int lowbit(unsigned long long);
// #endif

}
#if defined(MORE_BACKWARD) && MORE_BACKWARD < 20010615
using namespace math;
#endif
}

#endif

---