---

more/math/special_functions.h

Go to the documentation of this file.

//  more/special_functions.h -- special mathematical functions
//  Copyright (C) 1999--2001  Petter Urkedal (petter.urkedal@matfys.lth.se)

//  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: special_functions.h,v 1.1 2002/05/30 18:01:38 petter_urkedal Exp $


#ifndef MORE_SPECIAL_FUNCTIONS_H
#define MORE_SPECIAL_FUNCTIONS_H

#include <more/math/qtypes.h>
#include <more/math/math.h>
#include <more/math/complex.h>


//  OBS:  The hypergeometric and legP with non-integer indices does not
//  work.

namespace more {
namespace math {

  // until we have numeric_limits
  double mathfunc_precision(double*) { return 1e-10; }
  float mathfunc_precision(float*) { return 1e-6; }

  // Chose the most efficient of abs() and norm().
  template<typename T> inline T abs_or_norm(T const& x) { return abs(x); }
  template<typename T> inline T abs_or_norm(complex<T> const& x)
    { return norm(x); }

  extern "C" {
      double f1c(int);
      double f2c(int);
      double fhc(int);
      double f1h(int);
  }

  //  --- hypergeometric ---

  // T and U should be float, double, complex<float> or complex<double>
  // T may also be int or halfint.
  template<typename T, typename U> inline U
    nonterminating_hypergeometric(T alpha, T beta, T gamma, U z) {
        typedef typename norm_type_<T>::eval RT;
        typedef typename norm_type_<U>::eval RU;
        RT test = real(alpha+beta-gamma);
        if(abs(z) > RU(1) || test >= RT(1) || test >= RT(0) && z == U(1))
            throw runtime_error(
                "arguments to hypergeometric is out of range");

        U sum = one;
        U term = one;
        T fact = one;
        do {
            term *= alpha*beta*z;
            term /= gamma*fact;
            sum += term;
            alpha += T(1);
            beta += T(1);
            gamma += T(1);
            fact += T(1);
        } while(abs(term) >= abs(sum)*mathfunc_precision((RU*)0)
                && fact < T(500));
        assert(fact < T(500));
        return sum;
    }

  template<typename T, typename U> inline U
    terminating_hypergeometric(int alpha, T beta, T gamma, U z) {
        assert(alpha <= 0);
        U sum = one;
        U term = one;
        T fact = one;
        while(alpha) {
            term *= alpha*beta*z;
            term /= gamma*fact;
            sum += term;
            ++alpha;
            beta += T(1);
            gamma += T(1);
            fact += T(1);
        }
        return sum;
    }

  //  2-specialized
  template<typename T, typename U> inline U
    hypergeometric(halfint alpha, halfint beta, T gamma, U z) {
        if(mod1(alpha) == 0 && alpha <= 0)
            return terminating_hypergeometric(itrunc(alpha), beta, gamma, z);
        else if(mod1(beta) == 0 && beta <= 0)
            return terminating_hypergeometric(itrunc(beta), alpha, gamma, z);
        else
            return nontermating_hypergeometric(alpha, beta, gamma, z);
    }
  template<typename T, typename U> inline U
    hypergeometric(int alpha, int beta, T gamma, U z) {
        if(alpha <= 0)
            return terminating_hypergeometric(alpha, beta, gamma, z);
        if(beta <= 0)
            return terminating_hypergeometric(beta, alpha, gamma, z);
        else
            return nonterminating_hypergeometric(alpha, beta, gamma, z);
    }
  //  1-specialized
  template<typename T1, typename T2, typename U> inline U
    hypergeometric(halfint alpha, T1 beta, T2 gamma, U z) {
        if(mod1(alpha) == 0 && alpha <= 0)
            return terminating_hypergeometric(itrunc(alpha), beta, gamma, z);
        else
            return nontermating_hypergeometric(alpha, beta, gamma, z);
    }
  template<typename T1, typename T2, typename U> inline U
    hypergeometric(int alpha, T1 beta, T2 gamma, U z) {
        if(alpha <= 0)
            return terminating_hypergeometric(alpha, beta, gamma, z);
        else
            return nonterminating_hypergeometric(alpha, beta, gamma, z);
    }
  template<typename T1, typename T2, typename U> inline U
    hypergeometric(T1 alpha, int beta, T2 gamma, U z) {
        return hypergeometric(beta, alpha, gamma, z);
    }
  template<typename T1, typename T2, typename U> inline U
    hypergeometric(T1 alpha, halfint beta, T2 gamma, U z) {
        return hypergeometric(beta, alpha, gamma, z);
    }
  //  generic
  template<typename T1, typename T2, typename T3, typename U> inline U
    hypergeometric(T1 alpha, T2 beta, T3 gamma, U z) {
        return nonterminating_hypergeometric(alpha, beta, gamma, z);
    }


  //  --- gamma ---

  //  To extend lgamma (if present) to the complex plane,
  template<typename T> inline T lgamma(T z) {
      if(real(z) < 10.0) {
          T q = 1.0;
          do {
              q *= z;
              z += 1.0;
          } while(real(z) < 10.0);
          return lgamma(z) - log(q);
      }

      // Stirling's series for z! = Gamma(z+1)

      z -= 1.0;
      T sum = (z+.5)*log(z) + .5*numbers::log_2pi - z;
      const T coeff[] = {
          .0833333333333333333333333333333333333333, //  B_2/2
          -.0027777777777777777777777777777777777777,// -B_4/(4*3)
          .0007936507936507936507936507936507936507, //  B_6/(6*5)
          -.0005952380952380952380952380952380952380,// -B_8/(8*7)
          .0008417508417508417508417508417508417508, //  B_10/(10*9)
          -.0019175269175269175269175269175269175269 // -B_12/(11*10)
      };
      T zinvn = 1.0/z;
      T zfact = zinvn*zinvn;
      sum += zinvn*coeff[0];
      for(int i = 1; i < 6; ++i) {
          zinvn *= zfact;
          sum += zinvn*coeff[i];
      }
      return sum;
  }
  template<typename T> inline T gamma(T z) { return exp(lgamma(z)); }
  inline int gamma(int i) { return factorial(i-1); }
  inline double lgamma(int i) { return log((double)factorial(i-1)); }


  //  --- associated Legendre functions ---

  template<typename T1, typename T2, typename U>
    inline U legP(T1 nu, T2 mu, U z) {
        return pow((z+U(1))/(z-U(1)), half*mu);
        return hypergeometric(-nu, nu+1, 1-mu, half*(U(1)-z))/(T2(1)-mu);
    }

//   template<typename Real>
//     inline Real legP(int m, int n, Real x) {
//      assert(m >= 0);
//      return (m%2? -1 : 1)*factorial(n+m, n-m)/factorial(m)
//          *pow(sqrt(1-x*x)/2, m);
//     }

  template<typename T>
    inline T legP(int n, int m, T x) {
        assert(n-m >= 0); // or: if(n-m < 0) return 0;
        int i0 = (n-m)%2;
        int n_plus_m_plus_i_minus_1 = n+m+i0-1;
        int n_minus_m = n-m;
        int n_minus_m_minus_i_plus_2 = n_minus_m-i0+2;
        T sqr_x = x*x;
        T term = n_plus_m_plus_i_minus_1 < 2? 1
            : /*double_factorial*/f2c(n_plus_m_plus_i_minus_1);
        term /= n_minus_m-i0 < 2? 1 : /*double_factorial*/f2c(n_minus_m-i0);
        term *= pow(1-sqr_x, half*m);
        if(i0) term *= x;
        if((n+m-i0)%4) term = -term;
        T sum = term;
        for(int i = i0+2; i <= n_minus_m; i += 2) {
            n_plus_m_plus_i_minus_1 += 2;
            n_minus_m_minus_i_plus_2 -= 2;
            term *= sqr_x;
            term *= -n_plus_m_plus_i_minus_1*n_minus_m_minus_i_plus_2;
            term /= i*(i-1);
            sum += term;
        }
        return sum;
    }
  //  These were used for testing.  Maybe we don't need them.
#ifdef MORE_ORTHOPOLY_SPECIAL_CASES
  template<typename T>
    inline T legP_0(T x) { return 1; }
  template<typename T>
    inline T legP_1(T x) { return x; }
  template<typename T>
    inline T legP_2(T x) { return T(1.5)*x*x-T(0.5); }
  template<typename T>
    inline T legP_3(T x) { return (T(2.5)*x*x-T(1.5))*x; }
  template<typename T>
    inline T legP_1_1(T x) { return -sqrt(T(1)-x*x); }
  template<typename T>
    inline T legP_2_1(T x) { return -T(3)*sqrt(T(1)-x*x)*x; }
  template<typename T>
    inline T legP_2_2(T x) { return T(3)*(T(1)-x*x); }
  template<typename T>
    inline T legP_3_1(T x) {
        T sqr_x = x*x;
        return -T(1.5)*sqrt(T(1)-sqr_x)*(T(5)*sqr_x-T(1));
    }
  template<typename T>
    inline T legP_3_2(T x) { return T(15)*(T(1)-x*x)*x; }
  template<typename T>
    inline T legP_3_3(T x) { return -T(15)*pow(T(1)-x*x, 3*half); }
#endif

  //  This could be specialized to avoid the extra computation of
  //  the factorials in legP.
  template<typename T>
    inline complex<typename norm_type_<T>::eval>
    sphY(int l, int m, T theta, T phi) {
        int abs_m = abs(m);
        T c = (2*l+1)/(numbers::p_4_pi*/*factorial(l+abs_m, l-abs_m)*/f1c(l+abs_m)/f1c(l-abs_m));
        c = sqrt(c);
        if(m%2) c = -c;
        return c*exp(onei*((T)m*phi))*legP(l, abs_m, cos(theta));
    }


  //  Hermite polynomials
  //  fixme: more efficient implementation
  template<typename T>
    inline T herH(int n, T const& x) {
#ifdef MORE_THROW_LOGIC_ERROR
        if(n < 0) throw std::domain_error
                      ("First index of herH must be non-negative.");
#endif
        if(n == 0) return T(1);
        else if(n == 1) return T(2)*x;
        else return T(2)*(x*herH(n-1, x) - T(n-1)*herH(n-2, x));
    }

  //  Laguerre polynomials; 1st index non-negative int, 2nd arbitrary
  //  fixme: recursive definitions are indeed easy...
  template<typename T, typename U>
    inline T lagL(int n, U alpha, T x) {
#ifdef MORE_THROW_LOGIC_ERROR
        if(n < 0) throw std::domain_error
                      ("First index of lagL must be non-negative.");
#endif
        if(n == 0) return T(1);
        else if(n == 1) return T(alpha+1)-x;
        else return ((T(2*n-1+alpha)-x)*lagL(n-1, alpha, x)
                     -T(n-1+alpha)*lagL(n-2, alpha, x))/T(n);
    }

  template<typename T> T besj(int, T);
}
  using namespace math;
}

#include "special_functions.tcc"
#endif

---