airy_der.c

开放gsl矩阵运算

   -0.0000001628202185026,
   -0.0000000090991442687,
   -0.0000000007438647126,
   -0.0000000000795494752,
   -0.0000000000104050944,
   -0.0000000000015932426,
   -0.0000000000002770648,
   -0.0000000000000535343,
   -0.0000000000000113062,
   -0.0000000000000025772,
   -0.0000000000000006278,
   -0.0000000000000001621,
   -0.0000000000000000441
};
static cheb_series an20_cs = {
  an20_data,
  15,
  -1, 1,
  8
};

/*
 series for aph2 on the interval -1.00000e+00 to -1.25000e-01
					with weighted error   2.94e-17
					 log weighted error  16.53
			       significant figures required  15.58
				    decimal places required  17.28
*/
static double aph2_data[32] = {
   -0.2057088719781465107,
    0.0422196961357771922,
    0.0020482560511207275,
    0.0002607800735165006,
    0.0000474824268004729,
    0.0000105102756431612,
    0.0000026353534014668,
    0.0000007208824863499,
    0.0000002103236664473,
    0.0000000644975634555,
    0.0000000205802377264,
    0.0000000067836273921,
    0.0000000022974015284,
    0.0000000007961306765,
    0.0000000002813860610,
    0.0000000001011749057,
    0.0000000000369306738,
    0.0000000000136615066,
    0.0000000000051142751,
    0.0000000000019351689,
    0.0000000000007393607,
    0.0000000000002849792,
    0.0000000000001107281,
    0.0000000000000433412,
    0.0000000000000170801,
    0.0000000000000067733,
    0.0000000000000027017,
    0.0000000000000010835,
    0.0000000000000004367,
    0.0000000000000001769,
    0.0000000000000000719,
    0.0000000000000000294
};
static cheb_series aph2_cs = {
  aph2_data,
  31,
  -1, 1,
  16
};

/*
 series for aph1 on the interval -1.25000e-01 to -1.56250e-02
					with weighted error   6.38e-17
					 log weighted error  16.20
			       significant figures required  14.91
				    decimal places required  16.87
*/
static double aph1_data[22] = {
  -0.1024172908077571694,
   0.0071697275146591248,
   0.0001209959363122329,
   0.0000073361512841220,
   0.0000007535382954272,
   0.0000001041478171741,
   0.0000000174358728519,
   0.0000000033399795033,
   0.0000000007073075174,
   0.0000000001619187515,
   0.0000000000394539982,
   0.0000000000101192282,
   0.0000000000027092778,
   0.0000000000007523806,
   0.0000000000002156369,
   0.0000000000000635283,
   0.0000000000000191757,
   0.0000000000000059143,
   0.0000000000000018597,
   0.0000000000000005950,
   0.0000000000000001934,
   0.0000000000000000638
};
static cheb_series aph1_cs = {
  aph1_data,
  21,
  -1, 1,
  10
};

/*
 series for aph0 on the interval -1.56250e-02 to  0.00000e+00
					with weighted error   2.29e-17
					 log weighted error  16.64
			       significant figures required  15.27
				    decimal places required  17.23
*/
static double aph0_data[15] = {
 -0.0855849241130933257,
  0.0011214378867065261,
  0.0000042721029353664,
  0.0000000817607381483,
  0.0000000033907645000,
  0.0000000002253264423,
  0.0000000000206284209,
  0.0000000000023858763,
  0.0000000000003301618,
  0.0000000000000527010,
  0.0000000000000094555,
  0.0000000000000018709,
  0.0000000000000004024,
  0.0000000000000000930,
  0.0000000000000000229
};
static cheb_series aph0_cs = {
  aph0_data,
  14,
  -1, 1,
  7
};


static
int
airy_deriv_mod_phase(const double x, gsl_mode_t mode,
                     gsl_sf_result * ampl, gsl_sf_result * phi)
{
  const double pi34 = 2.356194490192344928847;
  gsl_sf_result result_a;
  gsl_sf_result result_p;
  double a, p;
  double sqx;
  double x32;

  if(x <= -4.0) {
    double z = 128.0/(x*x*x) + 1.0;
    cheb_eval_mode_e(&an20_cs, z, mode, &result_a);
    cheb_eval_mode_e(&aph0_cs, z, mode, &result_p);
  }
  else if(x <= -2.0) {
    double z = (128.0/(x*x*x) + 9.0) / 7.0;
    cheb_eval_mode_e(&an21_cs, z, mode, &result_a);
    cheb_eval_mode_e(&aph1_cs, z, mode, &result_p);
  }
  else if(x <= -1.0) {
    double z = (16.0/(x*x*x) + 9.0) / 7.0;
    cheb_eval_mode_e(&an22_cs, z, mode, &result_a);
    cheb_eval_mode_e(&aph2_cs, z, mode, &result_p);
  }
  else {
    ampl->val = 0.0;
    ampl->err = 0.0;
    phi->val  = 0.0;
    phi->err  = 0.0;
    GSL_ERROR ("x is greater than 1.0", GSL_EDOM);
  }

  a =  0.3125 + result_a.val;
  p = -0.625  + result_p.val;
 
  sqx = sqrt(-x);
  x32   = x*sqx;

  ampl->val = sqrt(a * sqx);
  ampl->err = fabs(ampl->val) * (GSL_DBL_EPSILON + fabs(result_a.err/result_a.val));
  phi->val  = pi34 - x * sqx * p;
  phi->err = fabs(phi->val) * (GSL_DBL_EPSILON + fabs(result_p.err/result_p.val));

  return GSL_SUCCESS;
}


/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_airy_Ai_deriv_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x < -1.0) {
    gsl_sf_result a;
    gsl_sf_result p;
    int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);
    double c    = cos(p.val);
    result->val  = a.val * c;
    result->err  = fabs(result->val * p.err) + fabs(c * a.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return status_ap;
  }
  else if(x <= 1.0) {
    const double x3 = x*x*x;
    const double x2 = x*x;
    gsl_sf_result result_c0;
    gsl_sf_result result_c1;
    cheb_eval_mode_e(&aif_cs, x3, mode, &result_c0);
    cheb_eval_mode_e(&aig_cs, x3, mode, &result_c1);

    result->val  = (x2*(0.125 + result_c0.val) - result_c1.val) - 0.25;
    result->err  = fabs(x2*result_c0.val) + result_c1.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);

    if(x > GSL_ROOT3_DBL_EPSILON*GSL_ROOT3_DBL_EPSILON) {
      /* scale only if x is positive */
      double s = exp(2.0*x*sqrt(x)/3.0);
      result->val *= s;
      result->err *= s;
    }

    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double sqrtx = sqrt(x);
    const double z = (16.0/(x*sqrtx) - 9.0)/7.0;
    const double s = sqrt(sqrtx);
    gsl_sf_result result_c0;
    cheb_eval_mode_e(&aip1_cs, z, mode, &result_c0);
    result->val  = -(0.28125 + result_c0.val) * s;
    result->err  = result_c0.err * s;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double sqrtx = sqrt(x);
    const double z = 16.0/(x*sqrtx) - 1.0;
    const double s = sqrt(sqrtx);
    gsl_sf_result result_c0;
    cheb_eval_mode_e(&aip2_cs, z, mode, &result_c0);
    result->val  = -(0.28125 + result_c0.val) * s;
    result->err  = result_c0.err * s;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}


int
gsl_sf_airy_Ai_deriv_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x < -1.0) {
    gsl_sf_result a;
    gsl_sf_result p;
    int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);
    double c    = cos(p.val);
    result->val  = a.val * c;
    result->err  = fabs(result->val * p.err) + fabs(c * a.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return status_ap;
  }
  else if(x < 1.0) {
    const double x3 = x*x*x;
    gsl_sf_result result_c1;
    gsl_sf_result result_c2;
    cheb_eval_mode_e(&aif_cs, x3, mode, &result_c1);
    cheb_eval_mode_e(&aig_cs, x3, mode, &result_c2);
    result->val  = (x*x*(0.125 + result_c1.val) - result_c2.val) - 0.25;
    result->err  = fabs(x*x*result_c1.err) + result_c2.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x*x*x < 9.0/4.0 * GSL_LOG_DBL_MIN*GSL_LOG_DBL_MIN) {
    gsl_sf_result result_aps;
    const double arg = -2.0*x*sqrt(x)/3.0;
    const int stat_a = gsl_sf_airy_Ai_deriv_scaled_e(x, mode, &result_aps);
    const int stat_e = gsl_sf_exp_mult_err_e(arg, 1.5*fabs(arg*GSL_DBL_EPSILON),
                                                result_aps.val, result_aps.err,
						result);
    return GSL_ERROR_SELECT_2(stat_e, stat_a);
  }
  else {
    UNDERFLOW_ERROR(result);
  }
}


int
gsl_sf_airy_Bi_deriv_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  const double atr =  8.7506905708484345;   /* 16./(sqrt(8)-1) */
  const double btr = -2.0938363213560543;   /* -(sqrt(8)+1)/(sqrt(8)-1) */

  /* CHECK_POINTER(result) */

  if(x < -1.0) {
    gsl_sf_result a;
    gsl_sf_result p;
    int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);
    double s     = sin(p.val);
    result->val  = a.val * s;
    result->err  = fabs(result->val * p.err) + fabs(s * a.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return status_ap;
  }
  else if(x < 1.0) {
    const double x3 = x*x*x;
    const double x2 = x*x;
    gsl_sf_result result_c1;
    gsl_sf_result result_c2;
    cheb_eval_mode_e(&bif_cs, x3, mode, &result_c1);
    cheb_eval_mode_e(&big_cs, x3, mode, &result_c2);
    result->val  = x2 * (result_c1.val + 0.25) + result_c2.val + 0.5;
    result->err  = x2 * result_c1.err + result_c2.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);

    if(x > GSL_ROOT3_DBL_EPSILON*GSL_ROOT3_DBL_EPSILON) {
      /* scale only if x is positive */
      const double s = exp(-2.0*x*sqrt(x)/3.0);
      result->val *= s;
      result->err *= s;
    }

    return GSL_SUCCESS;
  }
  else if(x < 2.0) {
    const double z = (2.0*x*x*x - 9.0) / 7.0;
    const double s = exp(-2.0*x*sqrt(x)/3.0);
    gsl_sf_result result_c0;
    gsl_sf_result result_c1;
    cheb_eval_mode_e(&bif2_cs, z, mode, &result_c0);
    cheb_eval_mode_e(&big2_cs, z, mode, &result_c1);
    result->val  = s * (x*x * (0.25 + result_c0.val) + 0.5 + result_c1.val);
    result->err  = s * (x*x * result_c0.err + result_c1.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < 4.0) {
    const double sqrtx = sqrt(x);
    const double z = atr/(x*sqrtx) + btr;
    const double s = sqrt(sqrtx);
    gsl_sf_result result_c0;
    cheb_eval_mode_e(&bip1_cs, z, mode, &result_c0);
    result->val  = s * (0.625 + result_c0.val);
    result->err  = s * result_c0.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double sqrtx = sqrt(x);
    const double z = 16.0/(x*sqrtx) - 1.0;
    const double s = sqrt(sqrtx);
    gsl_sf_result result_c0;
    cheb_eval_mode_e(&bip2_cs, z, mode, &result_c0);
    result->val  = s * (0.625 + result_c0.val);
    result->err  = s * result_c0.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}


int
gsl_sf_airy_Bi_deriv_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x < -1.0) {
    gsl_sf_result a;
    gsl_sf_result p;
    int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);
    double s    = sin(p.val);
    result->val  = a.val * s;
    result->err  = fabs(result->val * p.err) + fabs(s * a.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return status_ap;
  }
  else if(x < 1.0) {
    const double x3 = x*x*x;
    const double x2 = x*x;
    gsl_sf_result result_c1;
    gsl_sf_result result_c2;
    cheb_eval_mode_e(&bif_cs, x3, mode, &result_c1);
    cheb_eval_mode_e(&big_cs, x3, mode, &result_c2);
    result->val  = x2 * (result_c1.val + 0.25) + result_c2.val + 0.5;
    result->err  = x2 * result_c1.err + result_c2.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < 2.0) {
    const double z = (2.0*x*x*x - 9.0) / 7.0;
    gsl_sf_result result_c1;
    gsl_sf_result result_c2;
    cheb_eval_mode_e(&bif2_cs, z, mode, &result_c1);
    cheb_eval_mode_e(&big2_cs, z, mode, &result_c2);
    result->val  = x*x * (result_c1.val + 0.25) + result_c2.val + 0.5;
    result->err  = x*x * result_c1.err + result_c2.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < GSL_ROOT3_DBL_MAX*GSL_ROOT3_DBL_MAX) {
    gsl_sf_result result_bps;
    const double arg = 2.0*(x*sqrt(x)/3.0);
    int stat_b = gsl_sf_airy_Bi_deriv_scaled_e(x, mode, &result_bps);
    int stat_e = gsl_sf_exp_mult_err_e(arg, 1.5*fabs(arg*GSL_DBL_EPSILON),
                                          result_bps.val, result_bps.err,
				          result);
    return GSL_ERROR_SELECT_2(stat_e, stat_b);
  }
  else {
    OVERFLOW_ERROR(result);
  }
}

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

#include "eval.h"

double gsl_sf_airy_Ai_deriv_scaled(const double x, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_airy_Ai_deriv_scaled_e(x, mode, &result));
}

double gsl_sf_airy_Ai_deriv(const double x, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_airy_Ai_deriv_e(x, mode, &result));
}

double gsl_sf_airy_Bi_deriv_scaled(const double x, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_airy_Bi_deriv_scaled_e(x, mode, &result));
}

double gsl_sf_airy_Bi_deriv(const double x, gsl_mode_t mode)
{
  EVAL_RESULT(gsl_sf_airy_Bi_deriv_e(x, mode, &result));
}