legendre_con.c

开放gsl矩阵运算

        result->val  = pre * (E.val - K.val);
	result->err  = pre * (E.err + K.err);
	result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
        return stat_K;
      }
    }
  }
  else if(   (x <= 0.0 && lambda < 1000.0)
          || (x <  0.1 && lambda < 17.0)
	  || (x <  0.2 && lambda < 5.0 )
    ) {
    return conicalP_xlt1_hyperg_A(1.0, lambda, x, result);
  }
  else if(   (x <= 0.2 && lambda < 17.0)
          || (x <  1.5 && lambda < 20.0)
    ) {
    const double arg = fabs(x*x - 1.0);
    const double sgn = GSL_SIGN(1.0 - x);
    const double pre = 0.5*(lambda*lambda + 0.25) * sgn * sqrt(arg);
    gsl_sf_result F;
    int stat_F = gsl_sf_hyperg_2F1_conj_e(1.5, lambda, 2.0, (1.0-x)/2, &F);
    result->val  = pre * F.val;
    result->err  = fabs(pre) * F.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_F;
  }
  else if(1.5 <= x && lambda < GSL_MAX(x,20.0)) {
    gsl_sf_result P;
    double lm;
    int stat_P = gsl_sf_conicalP_large_x_e(1.0, lambda, x,
                                              &P, &lm
                                              );
    int stat_e = gsl_sf_exp_mult_err_e(lm, 2.0 * GSL_DBL_EPSILON * fabs(lm),
    					  P.val, P.err,
    					  result);
    return GSL_ERROR_SELECT_2(stat_e, stat_P);
  }
  else {
    double V0, V1;
    if(x < 1.0) {
      const double sqrt_1mx = sqrt(1.0 - x);
      const double sqrt_1px = sqrt(1.0 + x);
      const double th  = acos(x);
      const double sth = sqrt_1mx * sqrt_1px;  /* sin(th) */
      gsl_sf_result I0, I1;
      int stat_I0 = gsl_sf_bessel_I0_scaled_e(th * lambda, &I0);
      int stat_I1 = gsl_sf_bessel_I1_scaled_e(th * lambda, &I1);
      int stat_I  = GSL_ERROR_SELECT_2(stat_I0, stat_I1);
      int stat_V  = conicalP_1_V(th, x/sth, lambda, -1.0, &V0, &V1);
      double bessterm = V0 * I0.val + V1 * I1.val;
      double besserr  =  fabs(V0) * I0.err + fabs(V1) * I1.err
                       + 2.0 * GSL_DBL_EPSILON * fabs(V0 * I0.val)
		       + 2.0 * GSL_DBL_EPSILON * fabs(V1 * I1.val);
      double arg1 = th * lambda;
      double sqts = sqrt(th/sth);
      int stat_e = gsl_sf_exp_mult_err_e(arg1, 2.0 * GSL_DBL_EPSILON * fabs(arg1),
                                            sqts * bessterm, sqts * besserr,
					    result);
      result->err *= 1.0/sqrt_1mx;
      return GSL_ERROR_SELECT_3(stat_e, stat_V, stat_I);
    }
    else {
      const double sqrt_xm1 = sqrt(x - 1.0);
      const double sqrt_xp1 = sqrt(x + 1.0);
      const double sh = sqrt_xm1 * sqrt_xp1;  /* sinh(xi)      */
      const double xi = log(x + sh);          /* xi = acosh(x) */
      const double xi_lam = xi * lambda;
      gsl_sf_result J0, J1;
      const int stat_J0 = gsl_sf_bessel_J0_e(xi_lam, &J0);
      const int stat_J1 = gsl_sf_bessel_J1_e(xi_lam, &J1);
      const int stat_J  = GSL_ERROR_SELECT_2(stat_J0, stat_J1);
      const int stat_V  = conicalP_1_V(xi, x/sh, lambda, 1.0, &V0, &V1);
      const double bessterm = V0 * J0.val + V1 * J1.val;
      const double besserr  = fabs(V0) * J0.err + fabs(V1) * J1.err
                       + 512.0 * 2.0 * GSL_DBL_EPSILON * fabs(V0 * J0.val)
		       + 512.0 * 2.0 * GSL_DBL_EPSILON * fabs(V1 * J1.val)
                       + GSL_DBL_EPSILON * fabs(xi_lam * V0 * J1.val)
                       + GSL_DBL_EPSILON * fabs(xi_lam * V1 * J0.val);
      const double pre = sqrt(xi/sh);
      result->val  = pre * bessterm;
      result->err  = pre * besserr * sqrt_xp1 / sqrt_xm1;
      result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
      return GSL_ERROR_SELECT_2(stat_V, stat_J);
    }
  }
}


/* P^{1/2}_{-1/2 + I lambda} (x)
 * [Abramowitz+Stegun 8.6.8, 8.6.12]
 * checked OK [GJ] Fri May  8 12:24:36 MDT 1998 
 */
int gsl_sf_conicalP_half_e(const double lambda, const double x,
                              gsl_sf_result * result
                              )
{
  /* CHECK_POINTER(result) */

  if(x <= -1.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 1.0) {
    double err_amp = 1.0 + 1.0/(GSL_DBL_EPSILON + fabs(1.0-fabs(x)));
    double ac  = acos(x);
    double den = sqrt(sqrt(1.0-x)*sqrt(1.0+x));
    result->val  = Root_2OverPi_ / den * cosh(ac * lambda);
    result->err  = err_amp * 3.0 * GSL_DBL_EPSILON * fabs(result->val);
    result->err *= fabs(ac * lambda) + 1.0;
    return GSL_SUCCESS;
  }
  else if(x == 1.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    /* x > 1 */
    double err_amp = 1.0 + 1.0/(GSL_DBL_EPSILON + fabs(1.0-fabs(x)));
    double sq_term = sqrt(x-1.0)*sqrt(x+1.0);
    double ln_term = log(x + sq_term);
    double den = sqrt(sq_term);
    double carg_val = lambda * ln_term;
    double carg_err = 2.0 * GSL_DBL_EPSILON * fabs(carg_val);
    gsl_sf_result cos_result;
    int stat_cos = gsl_sf_cos_err_e(carg_val, carg_err, &cos_result);
    result->val  = Root_2OverPi_ / den * cos_result.val;
    result->err  = err_amp * Root_2OverPi_ / den * cos_result.err;
    result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_cos;
  }
}


/* P^{-1/2}_{-1/2 + I lambda} (x)
 * [Abramowitz+Stegun 8.6.9, 8.6.14]
 * checked OK [GJ] Fri May  8 12:24:43 MDT 1998 
 */
int gsl_sf_conicalP_mhalf_e(const double lambda, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= -1.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 1.0) {
    double ac  = acos(x);
    double den = sqrt(sqrt(1.0-x)*sqrt(1.0+x));
    double arg = ac * lambda;
    double err_amp = 1.0 + 1.0/(GSL_DBL_EPSILON + fabs(1.0-fabs(x)));
    if(fabs(arg) < GSL_SQRT_DBL_EPSILON) {
      result->val  = Root_2OverPi_ / den * ac;
      result->err  = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      result->err *= err_amp;
    }
    else {
      result->val  = Root_2OverPi_ / (den*lambda) * sinh(arg);
      result->err  = GSL_DBL_EPSILON * (fabs(arg)+1.0) * fabs(result->val);
      result->err *= err_amp;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    }
    return GSL_SUCCESS;
  }
  else if(x == 1.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    /* x > 1 */
    double sq_term = sqrt(x-1.0)*sqrt(x+1.0);
    double ln_term = log(x + sq_term);
    double den = sqrt(sq_term);
    double arg_val = lambda * ln_term;
    double arg_err = 2.0 * GSL_DBL_EPSILON * fabs(arg_val);
    if(arg_val < GSL_SQRT_DBL_EPSILON) {
      result->val = Root_2OverPi_ / den * ln_term;
      result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return GSL_SUCCESS;
    }
    else {
      gsl_sf_result sin_result;
      int stat_sin = gsl_sf_sin_err_e(arg_val, arg_err, &sin_result);
      result->val  = Root_2OverPi_ / (den*lambda) * sin_result.val;
      result->err  = Root_2OverPi_ / fabs(den*lambda) * sin_result.err;
      result->err += 3.0 * GSL_DBL_EPSILON * fabs(result->val);
      return stat_sin;
    }
  }
}


int gsl_sf_conicalP_sph_reg_e(const int l, const double lambda,
                                 const double x,
                                 gsl_sf_result * result
                                 )
{
  /* CHECK_POINTER(result) */

  if(x <= -1.0 || l < -1) {
    DOMAIN_ERROR(result);
  }
  else if(l == -1) {
    return gsl_sf_conicalP_half_e(lambda, x, result);
  }
  else if(l == 0) {
    return gsl_sf_conicalP_mhalf_e(lambda, x, result);
  }
  else if(x == 1.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x < 0.0) {
    double c = 1.0/sqrt(1.0-x*x);
    gsl_sf_result r_Pellm1;
    gsl_sf_result r_Pell;
    int stat_0 = gsl_sf_conicalP_half_e(lambda, x, &r_Pellm1);  /* P^( 1/2) */
    int stat_1 = gsl_sf_conicalP_mhalf_e(lambda, x, &r_Pell);   /* P^(-1/2) */
    int stat_P = GSL_ERROR_SELECT_2(stat_0, stat_1);
    double Pellm1 = r_Pellm1.val;
    double Pell   = r_Pell.val;
    double Pellp1;
    int ell;

    for(ell=0; ell<l; ell++) {
      double d = (ell+1.0)*(ell+1.0) + lambda*lambda;
      Pellp1 = (Pellm1 - (2.0*ell+1.0)*c*x * Pell) / d;
      Pellm1 = Pell;
      Pell   = Pellp1;
    }

    result->val  = Pell;
    result->err  = (0.5*l + 1.0) * GSL_DBL_EPSILON * fabs(Pell);
    result->err += GSL_DBL_EPSILON * l * fabs(result->val);
    return stat_P;
  }
  else if(x < 1.0) {
    const double xi = x/(sqrt(1.0-x)*sqrt(1.0+x));
    gsl_sf_result rat;
    gsl_sf_result Phf;
    int stat_CF1 = conicalP_negmu_xlt1_CF1(0.5, l, lambda, x, &rat);
    int stat_Phf = gsl_sf_conicalP_half_e(lambda, x, &Phf);
    double Pellp1 = rat.val * GSL_SQRT_DBL_MIN;
    double Pell   = GSL_SQRT_DBL_MIN;
    double Pellm1;
    int ell;

    for(ell=l; ell>=0; ell--) {
      double d = (ell+1.0)*(ell+1.0) + lambda*lambda;
      Pellm1 = (2.0*ell+1.0)*xi * Pell + d * Pellp1;
      Pellp1 = Pell;
      Pell   = Pellm1;
    }

    result->val  = GSL_SQRT_DBL_MIN * Phf.val / Pell;
    result->err  = GSL_SQRT_DBL_MIN * Phf.err / fabs(Pell);
    result->err += fabs(rat.err/rat.val) * (l + 1.0) * fabs(result->val);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_ERROR_SELECT_2(stat_Phf, stat_CF1);
  }
  else if(x == 1.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    /* x > 1.0 */

    const double xi = x/sqrt((x-1.0)*(x+1.0));
    gsl_sf_result rat;
    int stat_CF1 = conicalP_negmu_xgt1_CF1(0.5, l, lambda, x, &rat);
    int stat_P;
    double Pellp1 = rat.val * GSL_SQRT_DBL_MIN;
    double Pell   = GSL_SQRT_DBL_MIN;
    double Pellm1;
    int ell;

    for(ell=l; ell>=0; ell--) {
      double d = (ell+1.0)*(ell+1.0) + lambda*lambda;
      Pellm1 = (2.0*ell+1.0)*xi * Pell - d * Pellp1;
      Pellp1 = Pell;
      Pell   = Pellm1;
    }

    if(fabs(Pell) > fabs(Pellp1)){
      gsl_sf_result Phf;
      stat_P = gsl_sf_conicalP_half_e(lambda, x, &Phf);
      result->val  =       GSL_SQRT_DBL_MIN * Phf.val / Pell;
      result->err  = 2.0 * GSL_SQRT_DBL_MIN * Phf.err / fabs(Pell);
      result->err += 2.0 * fabs(rat.err/rat.val) * (l + 1.0) * fabs(result->val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    }
    else {
      gsl_sf_result Pmhf;
      stat_P = gsl_sf_conicalP_mhalf_e(lambda, x, &Pmhf);
      result->val  =       GSL_SQRT_DBL_MIN * Pmhf.val / Pellp1;
      result->err  = 2.0 * GSL_SQRT_DBL_MIN * Pmhf.err / fabs(Pellp1);
      result->err += 2.0 * fabs(rat.err/rat.val) * (l + 1.0) * fabs(result->val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    }

    return GSL_ERROR_SELECT_2(stat_P, stat_CF1);
  }
}


int gsl_sf_conicalP_cyl_reg_e(const int m, const double lambda,
                                 const double x,
                                 gsl_sf_result * result
                                 )
{
  /* CHECK_POINTER(result) */

  if(x <= -1.0 || m < -1) {
    DOMAIN_ERROR(result);
  }
  else if(m == -1) {
    return gsl_sf_conicalP_1_e(lambda, x, result);
  }
  else if(m == 0) {
    return gsl_sf_conicalP_0_e(lambda, x, result);
  }
  else if(x == 1.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x < 0.0) {
    double c = 1.0/sqrt(1.0-x*x);
    gsl_sf_result r_Pkm1;
    gsl_sf_result r_Pk;
    int stat_0 = gsl_sf_conicalP_1_e(lambda, x, &r_Pkm1);  /* P^1 */
    int stat_1 = gsl_sf_conicalP_0_e(lambda, x, &r_Pk);    /* P^0 */
    int stat_P = GSL_ERROR_SELECT_2(stat_0, stat_1);
    double Pkm1 = r_Pkm1.val;
    double Pk   = r_Pk.val;
    double Pkp1;
    int k;

    for(k=0; k<m; k++) {
      double d = (k+0.5)*(k+0.5) + lambda*lambda;
      Pkp1 = (Pkm1 - 2.0*k*c*x * Pk) / d;
      Pkm1 = Pk;
      Pk   = Pkp1;
    }

    result->val  = Pk;
    result->err  = (m + 2.0) * GSL_DBL_EPSILON * fabs(Pk);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

    return stat_P;
  }
  else if(x < 1.0) {
    const double xi = x/(sqrt(1.0-x)*sqrt(1.0+x));
    gsl_sf_result rat;
    gsl_sf_result P0;
    int stat_CF1 = conicalP_negmu_xlt1_CF1(0.0, m, lambda, x, &rat);
    int stat_P0  = gsl_sf_conicalP_0_e(lambda, x, &P0);
    double Pkp1 = rat.val * GSL_SQRT_DBL_MIN;
    double Pk   = GSL_SQRT_DBL_MIN;
    double Pkm1;
    int k;

    for(k=m; k>0; k--) {
      double d = (k+0.5)*(k+0.5) + lambda*lambda;
      Pkm1 = 2.0*k*xi * Pk + d * Pkp1;
      Pkp1 = Pk;
      Pk   = Pkm1;
    }

    result->val  = GSL_SQRT_DBL_MIN * P0.val / Pk;
    result->err  = 2.0 * GSL_SQRT_DBL_MIN * P0.err / fabs(Pk);
    result->err += 2.0 * fabs(rat.err/rat.val) * (m + 1.0) * fabs(result->val);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_ERROR_SELECT_2(stat_P0, stat_CF1);
  }
  else if(x == 1.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    /* x > 1.0 */

    const double xi = x/sqrt((x-1.0)*(x+1.0));
    gsl_sf_result rat;
    int stat_CF1 = conicalP_negmu_xgt1_CF1(0.0, m, lambda, x, &rat);
    int stat_P;
    double Pkp1 = rat.val * GSL_SQRT_DBL_MIN;
    double Pk   = GSL_SQRT_DBL_MIN;
    double Pkm1;
    int k;

    for(k=m; k>-1; k--) {
      double d = (k+0.5)*(k+0.5) + lambda*lambda;
      Pkm1 = 2.0*k*xi * Pk - d * Pkp1;
      Pkp1 = Pk;
      Pk   = Pkm1;
    }

    if(fabs(Pk) > fabs(Pkp1)){
      gsl_sf_result P1;
      stat_P = gsl_sf_conicalP_1_e(lambda, x, &P1);
      result->val  = GSL_SQRT_DBL_MIN * P1.val / Pk;
      result->err  = 2.0 * GSL_SQRT_DBL_MIN * P1.err / fabs(Pk);
      result->err += 2.0 * fabs(rat.err/rat.val) * (m+2.0) * fabs(result->val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    }
    else {
      gsl_sf_result P0;
      stat_P = gsl_sf_conicalP_0_e(lambda, x, &P0);
      result->val  = GSL_SQRT_DBL_MIN * P0.val / Pkp1;
      result->err  = 2.0 * GSL_SQRT_DBL_MIN * P0.err / fabs(Pkp1);
      result->err += 2.0 * fabs(rat.err/rat.val) * (m+2.0) * fabs(result->val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    }

    return GSL_ERROR_SELECT_2(stat_P, stat_CF1);
  }
}


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

#include "eval.h"

double gsl_sf_conicalP_0(const double lambda, const double x)
{
  EVAL_RESULT(gsl_sf_conicalP_0_e(lambda, x, &result));
}

double gsl_sf_conicalP_1(const double lambda, const double x)
{
  EVAL_RESULT(gsl_sf_conicalP_1_e(lambda, x, &result));
}

double gsl_sf_conicalP_half(const double lambda, const double x)
{
  EVAL_RESULT(gsl_sf_conicalP_half_e(lambda, x, &result));
}

double gsl_sf_conicalP_mhalf(const double lambda, const double x)
{
  EVAL_RESULT(gsl_sf_conicalP_mhalf_e(lambda, x, &result));
}

double gsl_sf_conicalP_sph_reg(const int l, const double lambda, const double x)
{
  EVAL_RESULT(gsl_sf_conicalP_sph_reg_e(l, lambda, x, &result));
}

double gsl_sf_conicalP_cyl_reg(const int m, const double lambda, const double x)
{
  EVAL_RESULT(gsl_sf_conicalP_cyl_reg_e(m, lambda, x, &result));
}