bessel_zero.c

多个常用的特殊函数的例子

 -5.5943e-13,
  1.4069e-13,
 -3.679e-14,
  1.119e-14,
 -4.99e-15,
  3.43e-15,
 -2.85e-15,
  2.3e-15,
 -1.7e-15,
  8.7e-16
};


/* Chebyshev expansion: j_{nu,14} = c_k T_k*(nu/28), nu <= 28 */
static const double coef_jnu14_a[] = {
  63.2256329577315566,
  19.5270342832914901,
 -0.42800190567884337,
  0.05859971627729398,
 -0.00987455163523582,
  0.00185641011402081,
 -0.00037352439419968,
  0.00007871886257265,
 -0.00001715728110045,
  3.83632624437e-6,
 -8.7518558668e-7,
  2.0291515353e-7,
 -4.767795233e-8,
  1.132844415e-8,
 -2.71734219e-9,
  6.5714886e-10,
 -1.6005342e-10,
  3.922557e-11,
 -9.66637e-12,
  2.39379e-12,
 -5.9541e-13,
  1.4868e-13,
 -3.726e-14,
  9.37e-15,
 -2.36e-15,
  6.0e-16
};


/* Chebyshev expansion: j_{nu,15} = c_k T_k*(nu/30), nu <= 30 */
static const double coef_jnu15_a[] = {
  67.7990939565631635,
  20.9237219226859859,
 -0.45840871823085836,
  0.06272149946755639,
 -0.01056229551143042,
  0.00198445078693100,
 -0.00039903958650729,
  0.00008404489865469,
 -0.00001830717574922,
  4.09103745566e-6,
 -9.3275533309e-7,
  2.1614056403e-7,
 -5.075725222e-8,
  1.205352081e-8,
 -2.88971837e-9,
  6.9846848e-10,
 -1.7002946e-10,
  4.164941e-11,
 -1.025859e-11,
  2.53921e-12,
 -6.3128e-13,
  1.5757e-13,
 -3.947e-14,
  9.92e-15,
 -2.50e-15,
  6.3e-16
};


/* Chebyshev expansion: j_{nu,16} = c_k T_k*(nu/32), nu <= 32 */
static const double coef_jnu16_a[] = {
  72.3725729616724770,
  22.32040402918608585,
 -0.48881449782358690,
  0.06684305681828766,
 -0.01124998690363398,
  0.00211247882775445,
 -0.00042455166484632,
  0.00008937015316346,
 -0.00001945687139551,
  4.34569739281e-6,
 -9.9031173548e-7,
  2.2936247195e-7,
 -5.383562595e-8,
  1.277835103e-8,
 -3.06202860e-9,
  7.3977037e-10,
 -1.8000071e-10,
  4.407196e-11,
 -1.085046e-11,
  2.68453e-12,
 -6.6712e-13,
  1.6644e-13,
 -4.168e-14,
  1.047e-14,
 -2.64e-15,
  6.7e-16
};


/* Chebyshev expansion: j_{nu,17} = c_k T_k*(nu/34), nu <= 34 */
static const double coef_jnu17_a[] = {
  76.9460667535209549,
  23.71708159112252670,
 -0.51921943142405352,
  0.07096442978067622,
 -0.01193763559341369,
  0.00224049662974902,
 -0.00045006122941781,
  0.00009469477941684,
 -0.00002060640777107,
  4.60031647195e-6,
 -1.04785755046e-6,
  2.4258161247e-7,
 -5.691327087e-8,
  1.350298805e-8,
 -3.23428733e-9,
  7.8105847e-10,
 -1.8996825e-10,
  4.649350e-11,
 -1.144205e-11,
  2.82979e-12,
 -7.0294e-13,
  1.7531e-13,
 -4.388e-14,
  1.102e-14,
 -2.78e-15,
  7.0e-16
};


/* Chebyshev expansion: j_{nu,18} = c_k T_k*(nu/36), nu <= 36 */
static const double coef_jnu18_a[] = {
  81.5195728368096659,
  25.11375537470259305,
 -0.54962366347317668,
  0.07508565026117689,
 -0.01262524908033818,
  0.00236850602019778,
 -0.00047556873651929,
  0.00010001889347161,
 -0.00002175581482429,
  4.85490251239e-6,
 -1.10539483940e-6,
  2.5579853343e-7,
 -5.999033352e-8,
  1.422747129e-8,
 -3.40650521e-9,
  8.2233565e-10,
 -1.9993286e-10,
  4.891426e-11,
 -1.203343e-11,
  2.97498e-12,
 -7.3875e-13,
  1.8418e-13,
 -4.608e-14,
  1.157e-14,
 -2.91e-15,
  7.4e-16
};


/* Chebyshev expansion: j_{nu,19} = c_k T_k*(nu/38), nu <= 38 */
static const double coef_jnu19_a[] = {
  86.0930892477047512,
  26.51042598308271729,
 -0.58002730731948358,
  0.07920674321589394,
 -0.01331283320930301,
  0.00249650841778073,
 -0.00050107453900793,
  0.00010534258471335,
 -0.00002290511552874,
  5.10946148897e-6,
 -1.16292517157e-6,
  2.6901365037e-7,
 -6.306692473e-8,
  1.495183048e-8,
 -3.57869025e-9,
  8.6360410e-10,
 -2.0989514e-10,
  5.133439e-11,
 -1.262465e-11,
  3.12013e-12,
 -7.7455e-13,
  1.9304e-13,
 -4.829e-14,
  1.212e-14,
 -3.05e-15,
  7.7e-16
};


/* Chebyshev expansion: j_{nu,20} = c_k T_k*(nu/40), nu <= 40 */
static const double coef_jnu20_a[] = {
  90.6666144195163770,
  27.9070938975436823,
 -0.61043045315390591,
  0.08332772844325554,
 -0.01400039260208282,
  0.00262450494035660,
 -0.00052657891389470,
  0.00011066592304919,
 -0.00002405432778364,
  5.36399803946e-6,
 -1.22044976064e-6,
  2.8222728362e-7,
 -6.614312964e-8,
  1.567608839e-8,
 -3.75084856e-9,
  9.0486546e-10,
 -2.1985553e-10,
  5.375401e-11,
 -1.321572e-11,
  3.26524e-12,
 -8.1033e-13,
  2.0190e-13,
 -5.049e-14,
  1.267e-14,
 -3.19e-15,
  8.0e-16,
 -2.0e-16
};


static const double * coef_jnu_a[] = {
  0,
  coef_jnu1_a,
  coef_jnu2_a,
  coef_jnu3_a,
  coef_jnu4_a,
  coef_jnu5_a,
  coef_jnu6_a,
  coef_jnu7_a,
  coef_jnu8_a,
  coef_jnu9_a,
  coef_jnu10_a,
  coef_jnu11_a,
  coef_jnu12_a,
  coef_jnu13_a,
  coef_jnu14_a,
  coef_jnu15_a,
  coef_jnu16_a,
  coef_jnu17_a,
  coef_jnu18_a,
  coef_jnu19_a,
  coef_jnu20_a
};

static const size_t size_jnu_a[] = {
  0,
  sizeof(coef_jnu1_a)/sizeof(double),
  sizeof(coef_jnu2_a)/sizeof(double),
  sizeof(coef_jnu3_a)/sizeof(double),
  sizeof(coef_jnu4_a)/sizeof(double),
  sizeof(coef_jnu5_a)/sizeof(double),
  sizeof(coef_jnu6_a)/sizeof(double),
  sizeof(coef_jnu7_a)/sizeof(double),
  sizeof(coef_jnu8_a)/sizeof(double),
  sizeof(coef_jnu9_a)/sizeof(double),
  sizeof(coef_jnu10_a)/sizeof(double),
  sizeof(coef_jnu11_a)/sizeof(double),
  sizeof(coef_jnu12_a)/sizeof(double),
  sizeof(coef_jnu13_a)/sizeof(double),
  sizeof(coef_jnu14_a)/sizeof(double),
  sizeof(coef_jnu15_a)/sizeof(double),
  sizeof(coef_jnu16_a)/sizeof(double),
  sizeof(coef_jnu17_a)/sizeof(double),
  sizeof(coef_jnu18_a)/sizeof(double),
  sizeof(coef_jnu19_a)/sizeof(double),
  sizeof(coef_jnu20_a)/sizeof(double)
};


static const double * coef_jnu_b[] = {
  0,
  coef_jnu1_b,
  coef_jnu2_b,
  coef_jnu3_b,
  coef_jnu4_b,
  coef_jnu5_b,
  coef_jnu6_b,
  coef_jnu7_b,
  coef_jnu8_b,
  coef_jnu9_b,
  coef_jnu10_b
};

static const size_t size_jnu_b[] = {
  0,
  sizeof(coef_jnu1_b)/sizeof(double),
  sizeof(coef_jnu2_b)/sizeof(double),
  sizeof(coef_jnu3_b)/sizeof(double),
  sizeof(coef_jnu4_b)/sizeof(double),
  sizeof(coef_jnu5_b)/sizeof(double),
  sizeof(coef_jnu6_b)/sizeof(double),
  sizeof(coef_jnu7_b)/sizeof(double),
  sizeof(coef_jnu8_b)/sizeof(double),
  sizeof(coef_jnu9_b)/sizeof(double),
  sizeof(coef_jnu10_b)/sizeof(double)
};



/* Evaluate Clenshaw recurrence for
 * a T* Chebyshev series.
 * sizeof(c) = N+1
 */
static double
clenshaw(const double * c, int N, double u)
{
  double B_np1 = 0.0;
  double B_n   = c[N];
  double B_nm1;
  int n;
  for(n=N; n>0; n--) {
    B_nm1 = 2.0*(2.0*u-1.0) * B_n - B_np1 + c[n-1];
    B_np1 = B_n;
    B_n   = B_nm1;
  }
  return B_n - (2.0*u-1.0)*B_np1;
}



/* correction terms to leading McMahon expansion
 * [Abramowitz+Stegun 9.5.12]
 * [Olver, Royal Society Math. Tables, v. 7]
 * We factor out a beta, so that this is a multiplicative
 * correction:
 *   j_{nu,s} = beta(s,nu) * mcmahon_correction(nu, beta(s,nu))
 *   macmahon_correction --> 1 as s --> Inf
 */
static double
mcmahon_correction(const double mu, const double beta)
{
  const double eb   = 8.0*beta;
  const double ebsq = eb*eb;

  if(mu < GSL_DBL_EPSILON) {
    /* Prevent division by zero below. */
    const double term1 =  1.0/ebsq;
    const double term2 = -4.0*31.0/(3*ebsq*ebsq);
    const double term3 =  32.0*3779.0/(15.0*ebsq*ebsq*ebsq);
    const double term4 = -64.0*6277237.0/(105.0*ebsq*ebsq*ebsq*ebsq);
    const double term5 =  512.0*2092163573.0/(315.0*ebsq*ebsq*ebsq*ebsq*ebsq);
    return 1.0 + 8.0*(term1 + term2 + term3 + term4 + term5);
  }
  else {
    /* Here we do things in terms of 1/mu, which
     * is purely to prevent overflow in the very
     * unlikely case that mu is really big.
     */
    const double mi   = 1.0/mu;
    const double r  = mu/ebsq;
    const double n2 = 4.0/3.0    * (7.0 - 31.0*mi);
    const double n3 = 32.0/15.0  * (83.0 + (-982.0 + 3779.0*mi)*mi);
    const double n4 = 64.0/105.0 * (6949.0 + (-153855.0 + (1585743.0 - 6277237.0*mi)*mi)*mi);
    const double n5 = 512.0/315.0 * (70197.0 + (-2479316.0 + (48010494.0 + (-512062548.0 + 2092163573.0*mi)*mi)*mi)*mi);
    const double n6 = 2048.0/3465.0 * (5592657.0 + (-287149133.0 + (8903961290.0 + (-179289628602.0 + (1982611456181.0 - 8249725736393.0*mi)*mi)*mi)*mi)*mi);
    const double term1 = (1.0 - mi) * r;
    const double term2 = term1 * n2 * r;
    const double term3 = term1 * n3 * r*r;
    const double term4 = term1 * n4 * r*r*r;
    const double term5 = term1 * n5 * r*r*r*r;
    const double term6 = term1 * n6 * r*r*r*r*r;
    return 1.0 - 8.0*(term1 + term2 + term3 + term4 + term5 + term6);
  }
}


/* Assumes z >= 1.0 */
static double
olver_b0(double z, double minus_zeta)
{
  if(z < 1.02) {
    const double a = 1.0-z;
    const double c0 =  0.0179988721413553309252458658183;
    const double c1 =  0.0111992982212877614645974276203;
    const double c2 =  0.0059404069786014304317781160605;
    const double c3 =  0.0028676724516390040844556450173;
    const double c4 =  0.0012339189052567271708525111185;
    const double c5 =  0.0004169250674535178764734660248;
    const double c6 =  0.0000330173385085949806952777365;
    const double c7 = -0.0001318076238578203009990106425;
    const double c8 = -0.0001906870370050847239813945647;
    return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*(c6 + a*(c7 + a*c8)))))));
  }
  else {
    const double abs_zeta = minus_zeta;
    const double t = 1.0/(z*sqrt(1.0 - 1.0/(z*z)));
    return -5.0/(48.0*abs_zeta*abs_zeta) + t*(3.0 + 5.0*t*t)/(24.0*sqrt(abs_zeta));
  }
}


inline
static double
olver_f1(double z, double minus_zeta)
{
  const double b0 = olver_b0(z, minus_zeta);
  const double h2 = sqrt(4.0*minus_zeta/(z*z-1.0)); /* FIXME */
  return 0.5 * z * h2 * b0;
}


int
gsl_sf_bessel_zero_J0_e(unsigned int s, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(s == 0){
    result->val = 0.0;
    result->err = 0.0;
    GSL_ERROR ("error", GSL_EINVAL);
  }
  else {
    /* See [F. Lether, J. Comp. Appl .Math. 67, 167 (1996)]. */

    static const double P[] = { 1567450796.0/12539606369.0,
                                8903660.0/2365861.0,
                                10747040.0/536751.0,
                                17590991.0/1696654.0
                              };
    static const double Q[] = { 1.0,
                                29354255.0/954518.0,
                                76900001.0/431847.0,
                                67237052.0/442411.0
                              };

    const double beta = (s - 0.25) * M_PI;
    const double bi2  = 1.0/(beta*beta);
    const double R33num = P[0] + bi2 * (P[1] + bi2 * (P[2] + P[3] * bi2));
    const double R33den = Q[0] + bi2 * (Q[1] + bi2 * (Q[2] + Q[3] * bi2));
    const double R33 = R33num/R33den;
    result->val = beta + R33/beta;
    result->err = fabs(3.0e-15 * result->val);
    return GSL_SUCCESS;
  }
}


int
gsl_sf_bessel_zero_J1_e(unsigned int s, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(s == 0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    /* See [M. Branders et al., J. Comp. Phys. 42, 403 (1981)]. */

    static const double a[] = { -0.362804405737084,
                                 0.120341279038597,
				 0.439454547101171e-01,
				 0.159340088474713e-02
                              };
    static const double b[] = {  1.0,
                                -0.325641790801361,
				-0.117453445968927,
				-0.424906902601794e-02
                              };

    const double beta = (s + 0.25) * M_PI;
    const double bi2  = 1.0/(beta*beta);
    const double Rnum = a[3] + bi2 * (a[2] + bi2 * (a[1] + bi2 * a[0]));
    const double Rden = b[3] + bi2 * (b[2] + bi2 * (b[1] + bi2 * b[0]));
    const double R = Rnum/Rden;
    result->val = beta * (1.0 + R*bi2);
    result->err = fabs(2.0e-14 * result->val);
    return GSL_SUCCESS;
  }
}


int
gsl_sf_bessel_zero_Jnu_e(double nu, unsigned int s, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(nu <= -1.0) {
    DOMAIN_ERROR(result);
  }
  else if(s == 0) {
    result->val = 0.0;
    result->err = 0.0;
    if (nu == 0.0) {
      GSL_ERROR ("no zero-th root for nu = 0.0", GSL_EINVAL);
    }
    return GSL_SUCCESS;
  }
  else if(nu < 0.0) {
    /* This can be done, I'm just lazy now. */
    result->val = 0.0;
    result->err = 0.0;
    GSL_ERROR("unimplemented", GSL_EUNIMPL);
  }
  else if(s == 1) {
    /* Chebyshev fits for the first positive zero.
     * For some reason Nemeth made this different from the others.
     */
    if(nu < 2.0) {
      const double * c = coef_jnu_a[s];
      const size_t   L = size_jnu_a[s];
      const double arg = nu/2.0;
      const double chb = clenshaw(c, L-1, arg);
      result->val = chb;
      result->err = 2.0e-15 * result->val;
    }
    else {
      const double * c = coef_jnu_b[s];
      const size_t   L = size_jnu_b[s];
      const double arg = pow(2.0/nu, 2.0/3.0);
      const double chb = clenshaw(c, L-1, arg);
      result->val = nu * chb;
      result->err = 2.0e-15 * result->val;
    }
    return GSL_SUCCESS;
  }
  else if(s <= 10) {
    /* Chebyshev fits for the first 10 positive zeros. */
    if(nu < s) {
      const double * c = coef_jnu_a[s];
      const size_t   L = size_jnu_a[s];
      const double arg = nu/s;
      const double chb = clenshaw(c, L-1, arg);
      result->val = chb;
      result->err = 2.0e-15 * result->val;
    }
    else {
      const double * c = coef_jnu_b[s];
      const size_t   L = size_jnu_b[s];
      const double arg = pow(s/nu, 2.0/3.0);
      const double chb = clenshaw(c, L-1, arg);
      result->val = nu * chb;
      result->err = 2.0e-15 * result->val;

      /* FIXME: truth in advertising for the screwed up
       * s = 5 fit. Need to fix that.
       */
      if(s == 5) {
        result->err *= 5.0e+06;
      }
    }
    return GSL_SUCCESS;
  }
  else if(s > 0.5*nu && s <= 20) {
    /* Chebyshev fits for 10 < s <= 20. */
    const double * c = coef_jnu_a[s];
    const size_t   L = size_jnu_a[s];
    const double arg = nu/(2.0*s);
    const double chb = clenshaw(c, L-1, arg);
    result->val = chb;
    result->err = 4.0e-15 * chb;
    return GSL_SUCCESS;
  }
  else if(s > 2.0 * nu) {
    /* McMahon expansion if s is large compared to nu. */
    const double beta = (s + 0.5*nu - 0.25) * M_PI;
    const double mc   = mcmahon_correction(4.0*nu*nu, beta);
    gsl_sf_result rat12;
    gsl_sf_pow_int_e(nu/beta, 14, &rat12);
    result->val  = beta * mc;
    result->err  = 4.0 * fabs(beta) * rat12.val;
    result->err += 4.0 * fabs(GSL_DBL_EPSILON * result->val);
    return GSL_SUCCESS;
  }
  else {
    /* Olver uniform asymptotic. */
    gsl_sf_result as;
    const int stat_as = gsl_sf_airy_zero_Ai_e(s, &as);
    const double minus_zeta = -pow(nu,-2.0/3.0) * as.val;
    const double z  = gsl_sf_bessel_Olver_zofmzeta(minus_zeta);
    const double f1 = olver_f1(z, minus_zeta);
    result->val  = nu * (z + f1/(nu*nu));
    result->err  = 0.001/(nu*nu*nu);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_as;
  }
}


/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_zero_J0(unsigned int s)
{
  EVAL_RESULT(gsl_sf_bessel_zero_J0_e(s, &result));
}

double gsl_sf_bessel_zero_J1(unsigned int s)
{
  EVAL_RESULT(gsl_sf_bessel_zero_J1_e(s, &result));
}

double gsl_sf_bessel_zero_Jnu(double nu, unsigned int s)
{
  EVAL_RESULT(gsl_sf_bessel_zero_Jnu_e(nu, s, &result));
}