bessel_olver.c

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

static double olver_A1(double z, double abs_zeta, double * err)
{
  if(z < 0.98) {
    double t = 1.0/sqrt(1.0-z*z);
    double rz = sqrt(abs_zeta);
    double t2 = t*t;
    double term1 =  t2*(81.0 - 462.0*t2 + 385.0*t2*t2)/1152.0;
    double term2 = -455.0/(4608.0*abs_zeta*abs_zeta*abs_zeta);
    double term3 =  7.0*t*(-3.0 + 5.0*t2)/(1152.0*rz*rz*rz);
    *err = 2.0 * GSL_DBL_EPSILON * (fabs(term1) + fabs(term2) + fabs(term3));
    return term1 + term2 + term3;
  }
  else if(z < 1.02) {
    const double a = 1.0-z;
    const double c0 = -0.00444444444444444444444444444444;
    const double c1 = -0.00184415584415584415584415584416;
    const double c2 =  0.00056812076812076812076812076812;
    const double c3 =  0.00168137865661675185484709294233;
    const double c4 =  0.00186744042139000122193399504324;
    const double c5 =  0.00161330105833747826430066790326;
    const double c6 =  0.00123177312220625816558607537838;
    const double c7 =  0.00087334711007377573881689318421;
    const double c8 =  0.00059004942455353250141217015410;
    const double sum = c0+a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*c8)))))));
    *err = 2.0 * GSL_DBL_EPSILON * fabs(sum);
    return sum;
  }
  else {
    const double t = 1.0/(z*sqrt(1.0 - 1.0/(z*z)));
    const double rz = sqrt(abs_zeta);
    const double t2 = t*t;
    const double term1 = -t2*(81.0 + 462.0*t2 + 385.0*t2*t2)/1152.0;
    const double term2 =  455.0/(4608.0*abs_zeta*abs_zeta*abs_zeta);
    const double term3 = -7.0*t*(3.0 + 5.0*t2)/(1152.0*rz*rz*rz);
    *err = 2.0 * GSL_DBL_EPSILON * (fabs(term1) + fabs(term2) + fabs(term3));
    return term1 + term2 + term3;
  }
}


static double olver_A2(double z, double abs_zeta)
{
  if(z < 0.88) {
    double t  = 1.0/sqrt(1.0-z*z);
    double t2 = t*t;
    double t4 = t2*t2;
    double t6 = t4*t2;
    double t8 = t4*t4;
    double rz = sqrt(abs_zeta);
    double z3 = abs_zeta*abs_zeta*abs_zeta;
    double z32 = rz*rz*rz;
    double z92 = z3*z32;
    double term1 = t4*(4465125.0 - 94121676.0*t2 + 349922430.0*t4 - 446185740.0*t6  + 185910725.0*t8)/39813120.0;
    double term2 = -40415375.0/(127401984.0*z3*z3);
    double term3 = -95095.0/15925248.0*t*(3.0-5.0*t2)/z92;
    double term4 = -455.0/5308416.0 *t2*(81.0 - 462.0*t2 + 385.0*t4)/z3;
    double term5 = -7.0/19906560.0*t*t2*(30375.0 - 369603.0*t2  + 765765.0*t4  - 425425.0*t6)/z32;
    return term1 + term2 + term3 + term4 + term5;
  }
  else if(z < 1.12) {
    double a = 1.0-z;
    const double c0  =  0.000693735541354588973636592684210;
    const double c1  =  0.000464483490365843307019777608010;
    const double c2  = -0.000289036254605598132482570468291;
    const double c3  = -0.000874764943953712638574497548110;
    const double c4  = -0.001029716376139865629968584679350;
    const double c5  = -0.000836857329713810600584714031650;
    const double c6  = -0.000488910893527218954998270124540;
    const double c7  = -0.000144236747940817220502256810151;
    const double c8  =  0.000114363800986163478038576460325;
    const double c9  =  0.000266806881492777536223944807117;
    const double c10 = -0.011975517576151069627471048587000;
    return c0+a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*(c9+a*c10)))))))));
  }
  else {
    const double t  = 1.0/(z*sqrt(1.0 - 1.0/(z*z)));
    const double t2 = t*t;
    const double t4 = t2*t2;
    const double t6 = t4*t2;
    const double t8 = t4*t4;
    const double rz = sqrt(abs_zeta);
    const double z3 = abs_zeta*abs_zeta*abs_zeta;
    const double z32 = rz*rz*rz;
    const double z92 = z3*z32;
    const double term1 = t4*(4465125.0 + 94121676.0*t2 + 349922430.0*t4 + 446185740.0*t6  + 185910725.0*t8)/39813120.0;
    const double term2 = -40415375.0/(127401984.0*z3*z3);
    const double term3 =  95095.0/15925248.0*t*(3.0+5.0*t2)/z92;
    const double term4 = -455.0/5308416.0 *t2*(81.0 + 462.0*t2 + 385.0*t4)/z3;
    const double term5 =  7.0/19906560.0*t*t2*(30375.0 + 369603.0*t2  + 765765.0*t4  + 425425.0*t6)/z32;
    return term1 + term2 + term3 + term4 + term5;
  }
}


static double olver_A3(double z, double abs_zeta)
{
  if(z < 0.9) {
    const double x = 20.0*z/9.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&A3_lt1_cs, x, &c);
    return c.val;
  }
  else if(z < 1.1) {
    double a = 1.0-z;
    const double c0 = -0.000354211971457743840771125759200;
    const double c1 = -0.000312322527890318832782774881353;
    const double c2 =  0.000277947465383133980329617631915;
    const double c3 =  0.000919803044747966977054155192400;
    const double c4 =  0.001147600388275977640983696906320;
    const double c5 =  0.000869239326123625742931772044544;
    const double c6 =  0.000287392257282507334785281718027;
    return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*c6)))));
  }
  else {
    const double x   = 11.0/(5.0*z) - 1.0;
    const double zi2 = 1.0/(z*z);
    gsl_sf_result c;
    cheb_eval_e(&A3_gt1_cs, x, &c);
    return  c.val * zi2*zi2*zi2;
  }
}


static double olver_A4(double z, double abs_zeta)
{
  if(z < 0.8) {
    const double x = 5.0*z/2.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&A4_lt1_cs, x, &c);
    return c.val;
  }
  else if(z < 1.2) {
    double a = 1.0-z;
    const double c0 =  0.00037819419920177291402661228437;
    const double c1 =  0.00040494390552363233477213857527;
    const double c2 = -0.00045764735528936113047289344569;
    const double c3 = -0.00165361044229650225813161341879;
    const double c4 = -0.00217527517983360049717137015539;
    const double c5 = -0.00152003287866490735107772795537;
    return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*c5))));
  }
  else {
    const double x   = 12.0/(5.0*z) - 1.0;
    const double zi2 = 1.0/(z*z);
    gsl_sf_result c;
    cheb_eval_e(&A4_gt1_cs, x, &c);
    return c.val * zi2*zi2*zi2*zi2;
  }
}

inline
static double olver_Asum(double nu, double z, double abs_zeta, double * err)
{
  double nu2 = nu*nu;
  double A1_err;
  double A1 = olver_A1(z, abs_zeta, &A1_err);
  double A2 = olver_A2(z, abs_zeta);
  double A3 = olver_A3(z, abs_zeta);
  double A4 = olver_A4(z, abs_zeta);
  *err = A1_err/nu2 + GSL_DBL_EPSILON;
  return 1.0 + A1/nu2 + A2/(nu2*nu2) + A3/(nu2*nu2*nu2) + A4/(nu2*nu2*nu2*nu2);
}

inline
static double olver_Bsum(double nu, double z, double abs_zeta)
{
  double nu2 = nu*nu;
  double B0 = olver_B0(z, abs_zeta);
  double B1 = olver_B1(z, abs_zeta);
  double B2 = olver_B2(z, abs_zeta);
  double B3 = olver_B3(z, abs_zeta);
  return B0 + B1/nu2 + B2/(nu2*nu2) + B3/(nu2*nu2*nu2*nu2);
}


/* uniform asymptotic, nu -> Inf, [Abramowitz+Stegun, 9.3.35]
 *
 * error:
 *    nu =  2: uniformly good to >  6D
 *    nu =  5: uniformly good to >  8D
 *    nu = 10: uniformly good to > 10D
 *    nu = 20: uniformly good to > 13D
 *
 */
int gsl_sf_bessel_Jnu_asymp_Olver_e(double nu, double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0 || nu <= 0.0) {
    DOMAIN_ERROR(result);
  }  
  else {
    double zeta, abs_zeta;
    double arg;
    double pre;
    double asum, bsum, asum_err;
    gsl_sf_result ai;
    gsl_sf_result aip;
    double z = x/nu;
    double crnu = pow(nu, 1.0/3.0);
    double nu3  = nu*nu*nu;
    double nu11 = nu3*nu3*nu3*nu*nu;
    int stat_a, stat_ap;

    if(fabs(1.0-z) < 0.02) {
      const double a = 1.0-z;
      const double c0 = 1.25992104989487316476721060728;
      const double c1 = 0.37797631496846194943016318218;
      const double c2 = 0.230385563409348235843147082474;
      const double c3 = 0.165909603649648694839821892031;
      const double c4 = 0.12931387086451008907;
      const double c5 = 0.10568046188858133991;
      const double c6 = 0.08916997952268186978;
      const double c7 = 0.07700014900618802456;
      pre = c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*(c6 + a*c7))))));
      zeta = a * pre;
      pre  = sqrt(2.0*sqrt(pre/(1.0+z)));
      abs_zeta = fabs(zeta);
    }
    else if(z < 1.0) {
      double rt   = sqrt(1.0 - z*z);
      abs_zeta = pow(1.5*(log((1.0+rt)/z) - rt), 2.0/3.0);
      zeta = abs_zeta;
      pre  = sqrt(2.0*sqrt(abs_zeta/(rt*rt)));
    }
    else {
      /* z > 1 */
      double rt = z * sqrt(1.0 - 1.0/(z*z));
      abs_zeta = pow(1.5*(rt - acos(1.0/z)), 2.0/3.0);
      zeta = -abs_zeta;
      pre  = sqrt(2.0*sqrt(abs_zeta/(rt*rt)));
    }

    asum = olver_Asum(nu, z, abs_zeta, &asum_err);
    bsum = olver_Bsum(nu, z, abs_zeta);

    arg  = crnu*crnu * zeta;
    stat_a  = gsl_sf_airy_Ai_e(arg, GSL_MODE_DEFAULT, &ai);
    stat_ap = gsl_sf_airy_Ai_deriv_e(arg, GSL_MODE_DEFAULT, &aip);

    result->val  = pre * (ai.val*asum/crnu + aip.val*bsum/(nu*crnu*crnu));
    result->err  = pre * (ai.err * fabs(asum/crnu));
    result->err += pre * fabs(ai.val) * asum_err / crnu;
    result->err += pre * fabs(ai.val * asum) / (crnu*nu11);
    result->err += 8.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_ERROR_SELECT_2(stat_a, stat_ap);
  }
}


/* uniform asymptotic, nu -> Inf,  [Abramowitz+Stegun, 9.3.36]
 *
 * error:
 *    nu =  2: uniformly good to >  6D
 *    nu =  5: uniformly good to >  8D
 *    nu = 10: uniformly good to > 10D
 *    nu = 20: uniformly good to > 13D
 */
int gsl_sf_bessel_Ynu_asymp_Olver_e(double nu, double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0 || nu <= 0.0) {
    DOMAIN_ERROR(result);
  }  
  else {
    double zeta, abs_zeta;
    double arg;
    double pre;
    double asum, bsum, asum_err;
    gsl_sf_result bi;
    gsl_sf_result bip;
    double z = x/nu;
    double crnu = pow(nu, 1.0/3.0);
    double nu3  = nu*nu*nu;
    double nu11 = nu3*nu3*nu3*nu*nu;
    int stat_b, stat_d;

    if(fabs(1.0-z) < 0.02) {
      const double a = 1.0-z;
      const double c0 = 1.25992104989487316476721060728;
      const double c1 = 0.37797631496846194943016318218;
      const double c2 = 0.230385563409348235843147082474;
      const double c3 = 0.165909603649648694839821892031;
      const double c4 = 0.12931387086451008907;
      const double c5 = 0.10568046188858133991;
      const double c6 = 0.08916997952268186978;
      const double c7 = 0.07700014900618802456;
      pre = c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*(c6 + a*c7))))));
      zeta = a * pre;
      pre  = sqrt(2.0*sqrt(pre/(1.0+z)));
      abs_zeta = fabs(zeta);
    }
    else if(z < 1.0) {
      double rt   = sqrt(1.0 - z*z);
      abs_zeta = pow(1.5*(log((1.0+rt)/z) - rt), 2.0/3.0);
      zeta = abs_zeta;
      pre  = sqrt(2.0*sqrt(abs_zeta/(rt*rt)));
    }
    else {
      /* z > 1 */
      double rt = z * sqrt(1.0 - 1.0/(z*z));
      double ac = acos(1.0/z);
      abs_zeta = pow(1.5*(rt - ac), 2.0/3.0);
      zeta = -abs_zeta;
      pre  = sqrt(2.0*sqrt(abs_zeta)/rt);
    }

    asum = olver_Asum(nu, z, abs_zeta, &asum_err);
    bsum = olver_Bsum(nu, z, abs_zeta);

    arg  = crnu*crnu * zeta;
    stat_b = gsl_sf_airy_Bi_e(arg, GSL_MODE_DEFAULT, &bi);
    stat_d = gsl_sf_airy_Bi_deriv_e(arg, GSL_MODE_DEFAULT, &bip);

    result->val  = -pre * (bi.val*asum/crnu + bip.val*bsum/(nu*crnu*crnu));
    result->err  =  pre * (bi.err * fabs(asum/crnu));
    result->err +=  pre * fabs(bi.val) * asum_err / crnu;
    result->err +=  pre * fabs(bi.val*asum) / (crnu*nu11);
    result->err +=  8.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_ERROR_SELECT_2(stat_b, stat_d);
  }
}