jv.cpp

强大的C++库

/*!
 * \file 
 * \brief Implementation of Bessel functions of noninteager order
 * \author Tony Ottosson
 *
 * $Date: 2006-08-07 10:38:19 +0200 (Mon, 07 Aug 2006) $
 * $Revision: 582 $
 *
 * -------------------------------------------------------------------------
 *
 * IT++ - C++ library of mathematical, signal processing, speech processing,
 *        and communications classes and functions
 *
 * Copyright (C) 1995-2006  (see AUTHORS file for a list of contributors)
 *
 * This program 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 program 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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 *
 * -------------------------------------------------------------------------
 *
 * This is slightly modified routine from the Cephes library:
 * http://www.netlib.org/cephes/ 
 */

#include <itpp/base/bessel/bessel_internal.h>
#include <itpp/base/scalfunc.h>
#include <cmath>

using namespace itpp;

/*
 * Bessel function of noninteger order
 *
 * double v, x, y, jv();
 *
 * y = jv( v, x );
 *
 * DESCRIPTION:
 *
 * Returns Bessel function of order v of the argument,
 * where v is real.  Negative x is allowed if v is an integer.
 *
 * Several expansions are included: the ascending power
 * series, the Hankel expansion, and two transitional
 * expansions for large v.  If v is not too large, it
 * is reduced by recurrence to a region of best accuracy.
 * The transitional expansions give 12D accuracy for v > 500.
 *
 * ACCURACY:
 * Results for integer v are indicated by *, where x and v
 * both vary from -125 to +125.  Otherwise,
 * x ranges from 0 to 125, v ranges as indicated by "domain."
 * Error criterion is absolute, except relative when |jv()| > 1.
 *
 * arithmetic  v domain  x domain    # trials      peak       rms
 *    IEEE      0,125     0,125      100000      4.6e-15    2.2e-16
 *    IEEE   -125,0       0,125       40000      5.4e-11    3.7e-13
 *    IEEE      0,500     0,500       20000      4.4e-15    4.0e-16
 * Integer v:
 *    IEEE   -125,125   -125,125      50000      3.5e-15*   1.9e-16*
 */

/*
  Cephes Math Library Release 2.8:  June, 2000
  Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
*/

static double recur(double *, double, double *, int);
static double jvs(double, double);
static double hankel(double, double);
static double jnx(double, double);
static double jnt(double, double);

#define MAXGAM 171.624376956302725
#define MAXNUM 1.79769313486231570815E308    /* 2**1024*(1-MACHEP) */
#define MACHEP 1.11022302462515654042E-16   /* 2**-53 */
#define MAXLOG 7.08396418532264106224E2     /* log 2**1022 */
#define MINLOG -7.08396418532264106224E2    /* log 2**-1022 */
#define PI 3.14159265358979323846       /* pi */

#define BIG  1.44115188075855872E+17



// ---------------------------- jv() -------------------------------------------------------
double jv(double n, double x)
{
  double k, q, t, y, an;
  int i, sign, nint;

  nint = 0;	/* Flag for integer n */
  sign = 1;	/* Flag for sign inversion */
  an = fabs( n );
  y = floor( an );
  if( y == an )
    {
      nint = 1;
      i = int(an - 16384.0 * floor( an/16384.0 ));
      if( n < 0.0 )
	{
	  if( i & 1 )
	    sign = -sign;
	  n = an;
	}
      if( x < 0.0 )
	{
	  if( i & 1 )
	    sign = -sign;
	  x = -x;
	}
      if( n == 0.0 ) // use 0th order bessel function
#ifdef _MSC_VER
	return _j0(x);
#else
	return j0(x);
#endif
      if( n == 1.0 ) // use 1th order bessel function
#ifdef _MSC_VER
	return (sign * _j1(x));
#else
	return (sign * j1(x));
#endif
    }

  it_error_if( (x < 0.0) && (y != an), "besselj:: negative values only allowed for integer orders.");

  y = fabs(x);

  if( y < MACHEP )
    goto underf;

  k = 3.6 * sqrt(y);
  t = 3.6 * sqrt(an);
  if( (y < t) && (an > 21.0) )
    return( sign * jvs(n,x) );
  if( (an < k) && (y > 21.0) )
    return( sign * hankel(n,x) );

  if( an < 500.0 )
    {
      /* Note: if x is too large, the continued
       * fraction will fail; but then the
       * Hankel expansion can be used.
       */
      if( nint != 0 )
	{
	  k = 0.0;
	  q = recur( &n, x, &k, 1 );
	  if( k == 0.0 )
	    {
#ifdef _MSC_VER
	      y = _j0(x)/q;
#else
	      y = j0(x)/q;
#endif
	      goto done;
	    }
	  if( k == 1.0 )
	    {
#ifdef _MSC_VER
	      y = _j1(x)/q;
#else
	      y = j1(x)/q;
#endif
	      goto done;
	    }
	}

      if( an > 2.0 * y )
	goto rlarger;

      if( (n >= 0.0) && (n < 20.0)
	  && (y > 6.0) && (y < 20.0) )
	{
	  /* Recur backwards from a larger value of n
	   */
	rlarger:
	  k = n;

	  y = y + an + 1.0;
	  if( y < 30.0 )
	    y = 30.0;
	  y = n + floor(y-n);
	  q = recur( &y, x, &k, 0 );
	  y = jvs(y,x) * q;
	  goto done;
	}

      if( k <= 30.0 )
	{
	  k = 2.0;
	}
      else if( k < 90.0 )
	{
	  k = (3*k)/4;
	}
      if( an > (k + 3.0) )
	{
	  if( n < 0.0 )
	    k = -k;
	  q = n - floor(n);
	  k = floor(k) + q;
	  if( n > 0.0 )
	    q = recur( &n, x, &k, 1 );
	  else
	    {
	      t = k;
	      k = n;
	      q = recur( &t, x, &k, 1 );
	      k = t;
	    }
	  if( q == 0.0 )
	    {
	    underf:
	      y = 0.0;
	      goto done;
	    }
	}
      else
	{
	  k = n;
	  q = 1.0;
	}

      /* boundary between convergence of
       * power series and Hankel expansion
       */
      y = fabs(k);
      if( y < 26.0 )
	t = (0.0083*y + 0.09)*y + 12.9;
      else
	t = 0.9 * y;

      if( x > t )
	y = hankel(k,x);
      else
	y = jvs(k,x);

      if( n > 0.0 )
	y /= q;
      else
	y *= q;
    }

  else
    {
      /* For large n, use the uniform expansion
       * or the transitional expansion.
       * But if x is of the order of n**2,
       * these may blow up, whereas the
       * Hankel expansion will then work.
       */

      if( n < 0.0 )
	{
	  it_warning("besselj:: partial loss of precision");
	  y = 0.0;
	  goto done;
	}
      t = x/n;
      t /= n;
      if( t > 0.3 )
	y = hankel(n,x);
      else
	y = jnx(n,x);
    }

 done:	return( sign * y);
}

/* Reduce the order by backward recurrence.
 * AMS55 #9.1.27 and 9.1.73.
 */

static double recur(double *n, double x, double *newn, int cancel)
{
  double pkm2, pkm1, pk, qkm2, qkm1;
  /* double pkp1; */
  double k, ans, qk, xk, yk, r, t, kf;
  static double big = BIG;
  int nflag, ctr;

  /* continued fraction for Jn(x)/Jn-1(x)  */
  if( *n < 0.0 )
    nflag = 1;
  else
    nflag = 0;

 fstart:
  pkm2 = 0.0;
  qkm2 = 1.0;
  pkm1 = x;
  qkm1 = *n + *n;
  xk = -x * x;
  yk = qkm1;
  ans = 1.0;
  ctr = 0;
  do
    {
      yk += 2.0;
      pk = pkm1 * yk +  pkm2 * xk;
      qk = qkm1 * yk +  qkm2 * xk;
      pkm2 = pkm1;
      pkm1 = pk;
      qkm2 = qkm1;
      qkm1 = qk;
      if( qk != 0 )
	r = pk/qk;
      else
	r = 0.0;
      if( r != 0 )
	{
	  t = fabs( (ans - r)/r );
	  ans = r;
	}
      else
	t = 1.0;

      if( ++ctr > 1000 )
	{
	  it_warning("besselj:: Underflow");
	  //mtherr( "jv", UNDERFLOW );
	  goto done;
	}

      if( t < MACHEP )
	goto done;

      if( fabs(pk) > big )
	{
	  pkm2 /= big;
	  pkm1 /= big;
	  qkm2 /= big;
	  qkm1 /= big;
	}
    }
  while( t > MACHEP );

 done:

  /* Change n to n-1 if n < 0 and the continued fraction is small
   */
  if( nflag > 0 )
    {
      if( fabs(ans) < 0.125 )
	{
	  nflag = -1;
	  *n = *n - 1.0;
	  goto fstart;
	}
    }


  kf = *newn;

  /* backward recurrence
   *              2k
   *  J   (x)  =  --- J (x)  -  J   (x)
   *   k-1         x   k         k+1
   */

  pk = 1.0;
  pkm1 = 1.0/ans;
  k = *n - 1.0;
  r = 2 * k;
  do
    {
      pkm2 = (pkm1 * r  -  pk * x) / x;
      /*	pkp1 = pk; */
      pk = pkm1;
      pkm1 = pkm2;
      r -= 2.0;
      /*
	t = fabs(pkp1) + fabs(pk);
	if( (k > (kf + 2.5)) && (fabs(pkm1) < 0.25*t) )
	{
	k -= 1.0;
	t = x*x;
	pkm2 = ( (r*(r+2.0)-t)*pk - r*x*pkp1 )/t;
	pkp1 = pk;
	pk = pkm1;
	pkm1 = pkm2;
	r -= 2.0;
	}
      */
      k -= 1.0;
    }
  while( k > (kf + 0.5) );

  /* Take the larger of the last two iterates
   * on the theory that it may have less cancellation error.
   */

  if( cancel )
    {
      if( (kf >= 0.0) && (fabs(pk) > fabs(pkm1)) )
	{
	  k += 1.0;
	  pkm2 = pk;
	}
    }
  *newn = k;

  return( pkm2 );
}



/* Ascending power series for Jv(x).
 * AMS55 #9.1.10.
 */