/* specfunc/airy_der.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * 
 * 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 Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_airy.h>

#include "error.h"

#include "chebyshev.h"
#include "cheb_eval_mode.c"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/


/* based on SLATEC aide.f, bide.f, aid.f, bid.f, r9admp.f */
 
/* 
 series for aif on the interval -1.00000e+00 to  1.00000e+00
                                        with weighted error   5.22e-18
                                         log weighted error  17.28
                               significant figures required  16.01
                                    decimal places required  17.73
*/
static double aif_data[8] = {
   0.10527461226531408809,
   0.01183613628152997844,
   0.00012328104173225664,
   0.00000062261225638140,
   0.00000000185298887844,
   0.00000000000363328873,
   0.00000000000000504622,
   0.00000000000000000522
};
static cheb_series aif_cs = {
  aif_data,
  7,
  -1, 1,
  7
};

/*
 series for aig on the interval -1.00000e+00 to  1.00000e+00
                                        with weighted error   3.14e-19
                                         log weighted error  18.50
                               significant figures required  17.44
                                    decimal places required  18.98
*/
static double aig_data[9] = {
   0.021233878150918666852,
   0.086315930335214406752,
   0.001797594720383231358,
   0.000014265499875550693,
   0.000000059437995283683,
   0.000000000152403366479,
   0.000000000000264587660,
   0.000000000000000331562,
   0.000000000000000000314
};
static cheb_series aig_cs = {
  aig_data,
  8,
  -1, 1,
  8
};

/*
 series for aip2 on the interval  0.00000e+00 to  1.25000e-01
                                        with weighted error   2.15e-17
                                         log weighted error  16.67
                               significant figures required  14.27
                                    decimal places required  17.26
*/
static double aip2_data[15] = {
    0.0065457691989713757,
    0.0023833724120774592,
   -0.0000430700770220586,
    0.0000015629125858629,
   -0.0000000815417186163,
    0.0000000054103738057,
   -0.0000000004284130883,
    0.0000000000389497963,
   -0.0000000000039623161,
    0.0000000000004428184,
   -0.0000000000000536297,
    0.0000000000000069650,
   -0.0000000000000009620,
    0.0000000000000001403,
   -0.0000000000000000215
};
static cheb_series aip2_cs = {
  aip2_data,
  14,
  -1, 1,
  9
};

/*
 series for aip1 on the interval  1.25000e-01 to  1.00000e+00
                                        with weighted error   2.60e-17
                                         log weighted error  16.58
                               significant figures required  14.91
                                    decimal places required  17.28
*/
static double aip1_data[25] = {
    0.0358865097808301538,
    0.0114668575627764899,
   -0.0007592073583861400,
    0.0000869517610893841,
   -0.0000128237294298592,
    0.0000022062695681038,
   -0.0000004222295185921,
    0.0000000874686415726,
   -0.0000000192773588418,
    0.0000000044668460054,
   -0.0000000010790108052,
    0.0000000002700029447,
   -0.0000000000696480108,
    0.0000000000184489907,
   -0.0000000000050027817,
    0.0000000000013852243,
   -0.0000000000003908218,
    0.0000000000001121536,
   -0.0000000000000326862,
    0.0000000000000096619,
   -0.0000000000000028935,
    0.0000000000000008770,
   -0.0000000000000002688,
    0.0000000000000000832,
   -0.0000000000000000260
};
static cheb_series aip1_cs = {
  aip1_data,
  24,
  -1, 1,
  14
};


/*
 series for bif on the interval -1.00000e+00 to  1.00000e+00
                                        with weighted error   9.05e-18
                                         log weighted error  17.04
                               significant figures required  15.83
                                    decimal places required  17.49
*/
static double bif_data[8] = {
   0.1153536790828570243,
   0.0205007894049192875,
   0.0002135290278902876,
   0.0000010783960614677,
   0.0000000032094708833,
   0.0000000000062930407,
   0.0000000000000087403,
   0.0000000000000000090
};
static cheb_series bif_cs = {
  bif_data,
  7,
  -1, 1,
  7
};

/*
 series for big on the interval -1.00000e+00 to  1.00000e+00
                                        with weighted error   5.44e-19
                                         log weighted error  18.26
                               significant figures required  17.46
                                    decimal places required  18.74
*/
static double big_data[9] = {
   -0.097196440416443537390,
    0.149503576843167066571,
    0.003113525387121326042,
    0.000024708570579821297,
    0.000000102949627731379,
    0.000000000263970373987,
    0.000000000000458279271,
    0.000000000000000574283,
    0.000000000000000000544
};
static cheb_series big_cs = {
  big_data,
  8,
  -1, 1,
  8
};

/*
 series for bif2 on the interval  1.00000e+00 to  8.00000e+00
                                        with weighted error   3.82e-19
                                         log weighted error  18.42
                               significant figures required  17.68
                                    decimal places required  18.92
*/
static double bif2_data[10] = {
   0.323493987603522033521,
   0.086297871535563559139,
   0.002994025552655397426,
   0.000051430528364661637,
   0.000000525840250036811,
   0.000000003561751373958,
   0.000000000017146864007,
   0.000000000000061663520,
   0.000000000000000171911,
   0.000000000000000000382
};
static cheb_series bif2_cs = {
  bif2_data,
  9,
  -1, 1,
  9
};

/*
 series for big2 on the interval  1.00000e+00 to  8.00000e+00
                                        with weighted error   3.35e-17
                                         log weighted error  16.48
                               significant figures required  16.52
                                    decimal places required  16.98
*/
static double big2_data[10] = {
   1.6062999463621294578,
   0.7449088819876088652,
   0.0470138738610277380,
   0.0012284422062548239,
   0.0000173222412256624,
   0.0000001521901652368,
   0.0000000009113560249,
   0.0000000000039547918,
   0.0000000000000130017,
   0.0000000000000000335
};
static cheb_series big2_cs = {
  big2_data,
  9,
  -1, 1,
  9
};

/*
 series for bip2 on the interval  0.00000e+00 to  1.25000e-01
                                        with weighted error   2.07e-18
                                         log weighted error  17.69
                               significant figures required  16.51
                                    decimal places required  18.42
*/
static double bip2_data[29] = {
    -0.13269705443526630495,
    -0.00568443626045977481,
    -0.00015643601119611610,
    -0.00001136737203679562,
    -0.00000143464350991284,
    -0.00000018098531185164,
     0.00000000926177343611,
     0.00000001710005490721,
     0.00000000476698163504,
    -0.00000000035195022023,
    -0.00000000058890614316,
    -0.00000000006678499608,
     0.00000000006395565102,
     0.00000000001554529427,
    -0.00000000000792397000,
    -0.00000000000258326243,
     0.00000000000121655048,
     0.00000000000038707207,
    -0.00000000000022487045,
    -0.00000000000004953477,
     0.00000000000004563782,
     0.00000000000000332998,
    -0.00000000000000921750,
     0.00000000000000094157,
     0.00000000000000167154,
    -0.00000000000000055134,
    -0.00000000000000022369,
     0.00000000000000017487,
     0.00000000000000000207
};
static cheb_series bip2_cs = {
  bip2_data,
  28,
  -1, 1,
  14
};

/*
 series for bip1 on the interval  1.25000e-01 to  3.53553e-01
                                        with weighted error   1.86e-17
                                         log weighted error  16.73
                               significant figures required  15.67
                                    decimal places required  17.42
*/
static double bip1_data[24] = {
   -0.1729187351079553719,
   -0.0149358492984694364,
   -0.0005471104951678566,
    0.0001537966292958408,
    0.0000154353476192179,
   -0.0000065434113851906,
    0.0000003728082407879,
    0.0000002072078388189,
   -0.0000000658173336470,
    0.0000000074926746354,
    0.0000000011101336884,
   -0.0000000007265140553,
    0.0000000001782723560,
   -0.0000000000217346352,
   -0.0000000000020302035,
    0.0000000000019311827,
   -0.0000000000006044953,
    0.0000000000001209450,
   -0.0000000000000125109,
   -0.0000000000000019917,
    0.0000000000000015154,
   -0.0000000000000004977,
    0.0000000000000001155,
   -0.0000000000000000186
};
static cheb_series bip1_cs = {
  bip1_data,
  23,
  -1, 1,
  13
};

/*
 series for an22 on the interval -1.00000e+00 to -1.25000e-01
                                        with weighted error   3.30e-17
                                         log weighted error  16.48
                               significant figures required  14.95
                                    decimal places required  17.24
*/
static double an22_data[33] = {
    0.0537418629629794329,
   -0.0126661435859883193,
   -0.0011924334106593007,
   -0.0002032327627275655,
   -0.0000446468963075164,
   -0.0000113359036053123,
   -0.0000031641352378546,
   -0.0000009446708886149,
   -0.0000002966562236472,
   -0.0000000969118892024,
   -0.0000000326822538653,
   -0.0000000113144618964,
   -0.0000000040042691002,
   -0.0000000014440333684,
   -0.0000000005292853746,
   -0.0000000001967763374,
   -0.0000000000740800096,
   -0.0000000000282016314,
   -0.0000000000108440066,
   -0.0000000000042074801,
   -0.0000000000016459150,
   -0.0000000000006486827,
   -0.0000000000002574095,
   -0.0000000000001027889,
   -0.0000000000000412846,
   -0.0000000000000166711,
   -0.0000000000000067657,
   -0.0000000000000027585,
   -0.0000000000000011296,
   -0.0000000000000004645,
   -0.0000000000000001917,
   -0.0000000000000000794,
   -0.0000000000000000330
};
static cheb_series an22_cs = {
  an22_data,
  32,
  -1, 1,
  18
};

/*
 series for an21 on the interval -1.25000e-01 to -1.56250e-02
                                        with weighted error   3.43e-17
                                         log weighted error  16.47
                               significant figures required  14.48
                                    decimal places required  17.16
*/
static double an21_data[24] = {
    0.0198313155263169394,
   -0.0029376249067087533,
   -0.0001136260695958196,
   -0.0000100554451087156,
   -0.0000013048787116563,
   -0.0000002123881993151,
   -0.0000000402270833384,
   -0.0000000084996745953,
   -0.0000000019514839426,
   -0.0000000004783865344,
   -0.0000000001236733992,
   -0.0000000000334137486,
   -0.0000000000093702824,
   -0.0000000000027130128,
   -0.0000000000008075954,
   -0.0000000000002463214,
   -0.0000000000000767656,
   -0.0000000000000243883,
   -0.0000000000000078831,
   -0.0000000000000025882,
   -0.0000000000000008619,
   -0.0000000000000002908,
   -0.0000000000000000993,
   -0.0000000000000000343
};
static cheb_series an21_cs = {
  an21_data,
  23,
  -1, 1,
  12
};

/*
 series for an20 on the interval -1.56250e-02 to  0.00000e+00
                                        with weighted error   4.41e-17
                                         log weighted error  16.36
                               significant figures required  14.16
                                    decimal places required  16.96
*/
static double an20_data[16] = {
    0.0126732217145738027,
   -0.0005212847072615621,
   -0.0000052672111140370,
   -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));
}