|
Packit |
67cb25 |
/* specfunc/airy.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: G. Jungman */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_trig.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_airy.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "error.h"
|
|
Packit |
67cb25 |
#include "check.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "chebyshev.h"
|
|
Packit |
67cb25 |
#include "cheb_eval_mode.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* chebyshev expansions for Airy modulus and phase
|
|
Packit |
67cb25 |
based on SLATEC r9aimp()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Series for AM21 on the interval -1.25000D-01 to 0.
|
|
Packit |
67cb25 |
with weighted error 2.89E-17
|
|
Packit |
67cb25 |
log weighted error 16.54
|
|
Packit |
67cb25 |
significant figures required 14.15
|
|
Packit |
67cb25 |
decimal places required 17.34
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Series for ATH1 on the interval -1.25000D-01 to 0.
|
|
Packit |
67cb25 |
with weighted error 2.53E-17
|
|
Packit |
67cb25 |
log weighted error 16.60
|
|
Packit |
67cb25 |
significant figures required 15.15
|
|
Packit |
67cb25 |
decimal places required 17.38
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Series for AM22 on the interval -1.00000D+00 to -1.25000D-01
|
|
Packit |
67cb25 |
with weighted error 2.99E-17
|
|
Packit |
67cb25 |
log weighted error 16.52
|
|
Packit |
67cb25 |
significant figures required 14.57
|
|
Packit |
67cb25 |
decimal places required 17.28
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Series for ATH2 on the interval -1.00000D+00 to -1.25000D-01
|
|
Packit |
67cb25 |
with weighted error 2.57E-17
|
|
Packit |
67cb25 |
log weighted error 16.59
|
|
Packit |
67cb25 |
significant figures required 15.07
|
|
Packit |
67cb25 |
decimal places required 17.34
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double am21_data[37] = {
|
|
Packit |
67cb25 |
0.0065809191761485,
|
|
Packit |
67cb25 |
0.0023675984685722,
|
|
Packit |
67cb25 |
0.0001324741670371,
|
|
Packit |
67cb25 |
0.0000157600904043,
|
|
Packit |
67cb25 |
0.0000027529702663,
|
|
Packit |
67cb25 |
0.0000006102679017,
|
|
Packit |
67cb25 |
0.0000001595088468,
|
|
Packit |
67cb25 |
0.0000000471033947,
|
|
Packit |
67cb25 |
0.0000000152933871,
|
|
Packit |
67cb25 |
0.0000000053590722,
|
|
Packit |
67cb25 |
0.0000000020000910,
|
|
Packit |
67cb25 |
0.0000000007872292,
|
|
Packit |
67cb25 |
0.0000000003243103,
|
|
Packit |
67cb25 |
0.0000000001390106,
|
|
Packit |
67cb25 |
0.0000000000617011,
|
|
Packit |
67cb25 |
0.0000000000282491,
|
|
Packit |
67cb25 |
0.0000000000132979,
|
|
Packit |
67cb25 |
0.0000000000064188,
|
|
Packit |
67cb25 |
0.0000000000031697,
|
|
Packit |
67cb25 |
0.0000000000015981,
|
|
Packit |
67cb25 |
0.0000000000008213,
|
|
Packit |
67cb25 |
0.0000000000004296,
|
|
Packit |
67cb25 |
0.0000000000002284,
|
|
Packit |
67cb25 |
0.0000000000001232,
|
|
Packit |
67cb25 |
0.0000000000000675,
|
|
Packit |
67cb25 |
0.0000000000000374,
|
|
Packit |
67cb25 |
0.0000000000000210,
|
|
Packit |
67cb25 |
0.0000000000000119,
|
|
Packit |
67cb25 |
0.0000000000000068,
|
|
Packit |
67cb25 |
0.0000000000000039,
|
|
Packit |
67cb25 |
0.0000000000000023,
|
|
Packit |
67cb25 |
0.0000000000000013,
|
|
Packit |
67cb25 |
0.0000000000000008,
|
|
Packit |
67cb25 |
0.0000000000000005,
|
|
Packit |
67cb25 |
0.0000000000000003,
|
|
Packit |
67cb25 |
0.0000000000000001,
|
|
Packit |
67cb25 |
0.0000000000000001
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series am21_cs = {
|
|
Packit |
67cb25 |
am21_data,
|
|
Packit |
67cb25 |
36,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
20
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double ath1_data[36] = {
|
|
Packit |
67cb25 |
-0.07125837815669365,
|
|
Packit |
67cb25 |
-0.00590471979831451,
|
|
Packit |
67cb25 |
-0.00012114544069499,
|
|
Packit |
67cb25 |
-0.00000988608542270,
|
|
Packit |
67cb25 |
-0.00000138084097352,
|
|
Packit |
67cb25 |
-0.00000026142640172,
|
|
Packit |
67cb25 |
-0.00000006050432589,
|
|
Packit |
67cb25 |
-0.00000001618436223,
|
|
Packit |
67cb25 |
-0.00000000483464911,
|
|
Packit |
67cb25 |
-0.00000000157655272,
|
|
Packit |
67cb25 |
-0.00000000055231518,
|
|
Packit |
67cb25 |
-0.00000000020545441,
|
|
Packit |
67cb25 |
-0.00000000008043412,
|
|
Packit |
67cb25 |
-0.00000000003291252,
|
|
Packit |
67cb25 |
-0.00000000001399875,
|
|
Packit |
67cb25 |
-0.00000000000616151,
|
|
Packit |
67cb25 |
-0.00000000000279614,
|
|
Packit |
67cb25 |
-0.00000000000130428,
|
|
Packit |
67cb25 |
-0.00000000000062373,
|
|
Packit |
67cb25 |
-0.00000000000030512,
|
|
Packit |
67cb25 |
-0.00000000000015239,
|
|
Packit |
67cb25 |
-0.00000000000007758,
|
|
Packit |
67cb25 |
-0.00000000000004020,
|
|
Packit |
67cb25 |
-0.00000000000002117,
|
|
Packit |
67cb25 |
-0.00000000000001132,
|
|
Packit |
67cb25 |
-0.00000000000000614,
|
|
Packit |
67cb25 |
-0.00000000000000337,
|
|
Packit |
67cb25 |
-0.00000000000000188,
|
|
Packit |
67cb25 |
-0.00000000000000105,
|
|
Packit |
67cb25 |
-0.00000000000000060,
|
|
Packit |
67cb25 |
-0.00000000000000034,
|
|
Packit |
67cb25 |
-0.00000000000000020,
|
|
Packit |
67cb25 |
-0.00000000000000011,
|
|
Packit |
67cb25 |
-0.00000000000000007,
|
|
Packit |
67cb25 |
-0.00000000000000004,
|
|
Packit |
67cb25 |
-0.00000000000000002
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series ath1_cs = {
|
|
Packit |
67cb25 |
ath1_data,
|
|
Packit |
67cb25 |
35,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
15
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double am22_data[33] = {
|
|
Packit |
67cb25 |
-0.01562844480625341,
|
|
Packit |
67cb25 |
0.00778336445239681,
|
|
Packit |
67cb25 |
0.00086705777047718,
|
|
Packit |
67cb25 |
0.00015696627315611,
|
|
Packit |
67cb25 |
0.00003563962571432,
|
|
Packit |
67cb25 |
0.00000924598335425,
|
|
Packit |
67cb25 |
0.00000262110161850,
|
|
Packit |
67cb25 |
0.00000079188221651,
|
|
Packit |
67cb25 |
0.00000025104152792,
|
|
Packit |
67cb25 |
0.00000008265223206,
|
|
Packit |
67cb25 |
0.00000002805711662,
|
|
Packit |
67cb25 |
0.00000000976821090,
|
|
Packit |
67cb25 |
0.00000000347407923,
|
|
Packit |
67cb25 |
0.00000000125828132,
|
|
Packit |
67cb25 |
0.00000000046298826,
|
|
Packit |
67cb25 |
0.00000000017272825,
|
|
Packit |
67cb25 |
0.00000000006523192,
|
|
Packit |
67cb25 |
0.00000000002490471,
|
|
Packit |
67cb25 |
0.00000000000960156,
|
|
Packit |
67cb25 |
0.00000000000373448,
|
|
Packit |
67cb25 |
0.00000000000146417,
|
|
Packit |
67cb25 |
0.00000000000057826,
|
|
Packit |
67cb25 |
0.00000000000022991,
|
|
Packit |
67cb25 |
0.00000000000009197,
|
|
Packit |
67cb25 |
0.00000000000003700,
|
|
Packit |
67cb25 |
0.00000000000001496,
|
|
Packit |
67cb25 |
0.00000000000000608,
|
|
Packit |
67cb25 |
0.00000000000000248,
|
|
Packit |
67cb25 |
0.00000000000000101,
|
|
Packit |
67cb25 |
0.00000000000000041,
|
|
Packit |
67cb25 |
0.00000000000000017,
|
|
Packit |
67cb25 |
0.00000000000000007,
|
|
Packit |
67cb25 |
0.00000000000000002
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series am22_cs = {
|
|
Packit |
67cb25 |
am22_data,
|
|
Packit |
67cb25 |
32,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
15
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double ath2_data[32] = {
|
|
Packit |
67cb25 |
0.00440527345871877,
|
|
Packit |
67cb25 |
-0.03042919452318455,
|
|
Packit |
67cb25 |
-0.00138565328377179,
|
|
Packit |
67cb25 |
-0.00018044439089549,
|
|
Packit |
67cb25 |
-0.00003380847108327,
|
|
Packit |
67cb25 |
-0.00000767818353522,
|
|
Packit |
67cb25 |
-0.00000196783944371,
|
|
Packit |
67cb25 |
-0.00000054837271158,
|
|
Packit |
67cb25 |
-0.00000016254615505,
|
|
Packit |
67cb25 |
-0.00000005053049981,
|
|
Packit |
67cb25 |
-0.00000001631580701,
|
|
Packit |
67cb25 |
-0.00000000543420411,
|
|
Packit |
67cb25 |
-0.00000000185739855,
|
|
Packit |
67cb25 |
-0.00000000064895120,
|
|
Packit |
67cb25 |
-0.00000000023105948,
|
|
Packit |
67cb25 |
-0.00000000008363282,
|
|
Packit |
67cb25 |
-0.00000000003071196,
|
|
Packit |
67cb25 |
-0.00000000001142367,
|
|
Packit |
67cb25 |
-0.00000000000429811,
|
|
Packit |
67cb25 |
-0.00000000000163389,
|
|
Packit |
67cb25 |
-0.00000000000062693,
|
|
Packit |
67cb25 |
-0.00000000000024260,
|
|
Packit |
67cb25 |
-0.00000000000009461,
|
|
Packit |
67cb25 |
-0.00000000000003716,
|
|
Packit |
67cb25 |
-0.00000000000001469,
|
|
Packit |
67cb25 |
-0.00000000000000584,
|
|
Packit |
67cb25 |
-0.00000000000000233,
|
|
Packit |
67cb25 |
-0.00000000000000093,
|
|
Packit |
67cb25 |
-0.00000000000000037,
|
|
Packit |
67cb25 |
-0.00000000000000015,
|
|
Packit |
67cb25 |
-0.00000000000000006,
|
|
Packit |
67cb25 |
-0.00000000000000002
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series ath2_cs = {
|
|
Packit |
67cb25 |
ath2_data,
|
|
Packit |
67cb25 |
31,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
16
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Airy modulus and phase for x < -1 */
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
airy_mod_phase(const double x, gsl_mode_t mode, gsl_sf_result * mod, gsl_sf_result * phase)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_sf_result result_m;
|
|
Packit |
67cb25 |
gsl_sf_result result_p;
|
|
Packit |
67cb25 |
double m, p;
|
|
Packit |
67cb25 |
double sqx;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < -2.0) {
|
|
Packit |
67cb25 |
double z = 16.0/(x*x*x) + 1.0;
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&am21_cs, z, mode, &result_m);
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&ath1_cs, z, mode, &result_p);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= -1.0) {
|
|
Packit |
67cb25 |
double z = (16.0/(x*x*x) + 9.0)/7.0;
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&am22_cs, z, mode, &result_m);
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&ath2_cs, z, mode, &result_p);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
mod->val = 0.0;
|
|
Packit |
67cb25 |
mod->err = 0.0;
|
|
Packit |
67cb25 |
phase->val = 0.0;
|
|
Packit |
67cb25 |
phase->err = 0.0;
|
|
Packit |
67cb25 |
GSL_ERROR ("x is greater than 1.0", GSL_EDOM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m = 0.3125 + result_m.val;
|
|
Packit |
67cb25 |
p = -0.625 + result_p.val;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
sqx = sqrt(-x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
mod->val = sqrt(m/sqx);
|
|
Packit |
67cb25 |
mod->err = fabs(mod->val) * (GSL_DBL_EPSILON + fabs(result_m.err/result_m.val));
|
|
Packit |
67cb25 |
phase->val = M_PI_4 - x*sqx * p;
|
|
Packit |
67cb25 |
phase->err = fabs(phase->val) * (GSL_DBL_EPSILON + fabs(result_p.err/result_p.val));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Chebyshev series for Ai(x) with x in [-1,1]
|
|
Packit |
67cb25 |
based on SLATEC ai(x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
series for aif on the interval -1.00000d+00 to 1.00000d+00
|
|
Packit |
67cb25 |
with weighted error 1.09e-19
|
|
Packit |
67cb25 |
log weighted error 18.96
|
|
Packit |
67cb25 |
significant figures required 17.76
|
|
Packit |
67cb25 |
decimal places required 19.44
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
series for aig on the interval -1.00000d+00 to 1.00000d+00
|
|
Packit |
67cb25 |
with weighted error 1.51e-17
|
|
Packit |
67cb25 |
log weighted error 16.82
|
|
Packit |
67cb25 |
significant figures required 15.19
|
|
Packit |
67cb25 |
decimal places required 17.27
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double ai_data_f[9] = {
|
|
Packit |
67cb25 |
-0.03797135849666999750,
|
|
Packit |
67cb25 |
0.05919188853726363857,
|
|
Packit |
67cb25 |
0.00098629280577279975,
|
|
Packit |
67cb25 |
0.00000684884381907656,
|
|
Packit |
67cb25 |
0.00000002594202596219,
|
|
Packit |
67cb25 |
0.00000000006176612774,
|
|
Packit |
67cb25 |
0.00000000000010092454,
|
|
Packit |
67cb25 |
0.00000000000000012014,
|
|
Packit |
67cb25 |
0.00000000000000000010
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series aif_cs = {
|
|
Packit |
67cb25 |
ai_data_f,
|
|
Packit |
67cb25 |
8,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
8
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double ai_data_g[8] = {
|
|
Packit |
67cb25 |
0.01815236558116127,
|
|
Packit |
67cb25 |
0.02157256316601076,
|
|
Packit |
67cb25 |
0.00025678356987483,
|
|
Packit |
67cb25 |
0.00000142652141197,
|
|
Packit |
67cb25 |
0.00000000457211492,
|
|
Packit |
67cb25 |
0.00000000000952517,
|
|
Packit |
67cb25 |
0.00000000000001392,
|
|
Packit |
67cb25 |
0.00000000000000001
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series aig_cs = {
|
|
Packit |
67cb25 |
ai_data_g,
|
|
Packit |
67cb25 |
7,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
7
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Chebvyshev series for Bi(x) with x in [-1,1]
|
|
Packit |
67cb25 |
based on SLATEC bi(x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
series for bif on the interval -1.00000d+00 to 1.00000d+00
|
|
Packit |
67cb25 |
with weighted error 1.88e-19
|
|
Packit |
67cb25 |
log weighted error 18.72
|
|
Packit |
67cb25 |
significant figures required 17.74
|
|
Packit |
67cb25 |
decimal places required 19.20
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
series for big on the interval -1.00000d+00 to 1.00000d+00
|
|
Packit |
67cb25 |
with weighted error 2.61e-17
|
|
Packit |
67cb25 |
log weighted error 16.58
|
|
Packit |
67cb25 |
significant figures required 15.17
|
|
Packit |
67cb25 |
decimal places required 17.03
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double data_bif[9] = {
|
|
Packit |
67cb25 |
-0.01673021647198664948,
|
|
Packit |
67cb25 |
0.10252335834249445610,
|
|
Packit |
67cb25 |
0.00170830925073815165,
|
|
Packit |
67cb25 |
0.00001186254546774468,
|
|
Packit |
67cb25 |
0.00000004493290701779,
|
|
Packit |
67cb25 |
0.00000000010698207143,
|
|
Packit |
67cb25 |
0.00000000000017480643,
|
|
Packit |
67cb25 |
0.00000000000000020810,
|
|
Packit |
67cb25 |
0.00000000000000000018
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series bif_cs = {
|
|
Packit |
67cb25 |
data_bif,
|
|
Packit |
67cb25 |
8,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
8
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double data_big[8] = {
|
|
Packit |
67cb25 |
0.02246622324857452,
|
|
Packit |
67cb25 |
0.03736477545301955,
|
|
Packit |
67cb25 |
0.00044476218957212,
|
|
Packit |
67cb25 |
0.00000247080756363,
|
|
Packit |
67cb25 |
0.00000000791913533,
|
|
Packit |
67cb25 |
0.00000000001649807,
|
|
Packit |
67cb25 |
0.00000000000002411,
|
|
Packit |
67cb25 |
0.00000000000000002
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series big_cs = {
|
|
Packit |
67cb25 |
data_big,
|
|
Packit |
67cb25 |
7,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
7
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Chebyshev series for Bi(x) with x in [1,8]
|
|
Packit |
67cb25 |
based on SLATEC bi(x)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double data_bif2[10] = {
|
|
Packit |
67cb25 |
0.0998457269381604100,
|
|
Packit |
67cb25 |
0.4786249778630055380,
|
|
Packit |
67cb25 |
0.0251552119604330118,
|
|
Packit |
67cb25 |
0.0005820693885232645,
|
|
Packit |
67cb25 |
0.0000074997659644377,
|
|
Packit |
67cb25 |
0.0000000613460287034,
|
|
Packit |
67cb25 |
0.0000000003462753885,
|
|
Packit |
67cb25 |
0.0000000000014288910,
|
|
Packit |
67cb25 |
0.0000000000000044962,
|
|
Packit |
67cb25 |
0.0000000000000000111
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series bif2_cs = {
|
|
Packit |
67cb25 |
data_bif2,
|
|
Packit |
67cb25 |
9,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
9
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double data_big2[10] = {
|
|
Packit |
67cb25 |
0.033305662145514340,
|
|
Packit |
67cb25 |
0.161309215123197068,
|
|
Packit |
67cb25 |
0.0063190073096134286,
|
|
Packit |
67cb25 |
0.0001187904568162517,
|
|
Packit |
67cb25 |
0.0000013045345886200,
|
|
Packit |
67cb25 |
0.0000000093741259955,
|
|
Packit |
67cb25 |
0.0000000000474580188,
|
|
Packit |
67cb25 |
0.0000000000001783107,
|
|
Packit |
67cb25 |
0.0000000000000005167,
|
|
Packit |
67cb25 |
0.0000000000000000011
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series big2_cs = {
|
|
Packit |
67cb25 |
data_big2,
|
|
Packit |
67cb25 |
9,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
9
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* chebyshev for Ai(x) asymptotic factor
|
|
Packit |
67cb25 |
based on SLATEC aie()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Series for AIP on the interval 0. to 1.00000D+00
|
|
Packit |
67cb25 |
with weighted error 5.10E-17
|
|
Packit |
67cb25 |
log weighted error 16.29
|
|
Packit |
67cb25 |
significant figures required 14.41
|
|
Packit |
67cb25 |
decimal places required 17.06
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
[GJ] Sun Apr 19 18:14:31 EDT 1998
|
|
Packit |
67cb25 |
There was something wrong with these coefficients. I was getting
|
|
Packit |
67cb25 |
errors after 3 or 4 digits. So I recomputed this table. Now I get
|
|
Packit |
67cb25 |
double precision agreement with Mathematica. But it does not seem
|
|
Packit |
67cb25 |
possible that the small differences here would account for the
|
|
Packit |
67cb25 |
original discrepancy. There must have been something wrong with my
|
|
Packit |
67cb25 |
original usage...
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double data_aip[36] = {
|
|
Packit |
67cb25 |
-0.0187519297793867540198,
|
|
Packit |
67cb25 |
-0.0091443848250055004725,
|
|
Packit |
67cb25 |
0.0009010457337825074652,
|
|
Packit |
67cb25 |
-0.0001394184127221491507,
|
|
Packit |
67cb25 |
0.0000273815815785209370,
|
|
Packit |
67cb25 |
-0.0000062750421119959424,
|
|
Packit |
67cb25 |
0.0000016064844184831521,
|
|
Packit |
67cb25 |
-0.0000004476392158510354,
|
|
Packit |
67cb25 |
0.0000001334635874651668,
|
|
Packit |
67cb25 |
-0.0000000420735334263215,
|
|
Packit |
67cb25 |
0.0000000139021990246364,
|
|
Packit |
67cb25 |
-0.0000000047831848068048,
|
|
Packit |
67cb25 |
0.0000000017047897907465,
|
|
Packit |
67cb25 |
-0.0000000006268389576018,
|
|
Packit |
67cb25 |
0.0000000002369824276612,
|
|
Packit |
67cb25 |
-0.0000000000918641139267,
|
|
Packit |
67cb25 |
0.0000000000364278543037,
|
|
Packit |
67cb25 |
-0.0000000000147475551725,
|
|
Packit |
67cb25 |
0.0000000000060851006556,
|
|
Packit |
67cb25 |
-0.0000000000025552772234,
|
|
Packit |
67cb25 |
0.0000000000010906187250,
|
|
Packit |
67cb25 |
-0.0000000000004725870319,
|
|
Packit |
67cb25 |
0.0000000000002076969064,
|
|
Packit |
67cb25 |
-0.0000000000000924976214,
|
|
Packit |
67cb25 |
0.0000000000000417096723,
|
|
Packit |
67cb25 |
-0.0000000000000190299093,
|
|
Packit |
67cb25 |
0.0000000000000087790676,
|
|
Packit |
67cb25 |
-0.0000000000000040927557,
|
|
Packit |
67cb25 |
0.0000000000000019271068,
|
|
Packit |
67cb25 |
-0.0000000000000009160199,
|
|
Packit |
67cb25 |
0.0000000000000004393567,
|
|
Packit |
67cb25 |
-0.0000000000000002125503,
|
|
Packit |
67cb25 |
0.0000000000000001036735,
|
|
Packit |
67cb25 |
-0.0000000000000000509642,
|
|
Packit |
67cb25 |
0.0000000000000000252377,
|
|
Packit |
67cb25 |
-0.0000000000000000125793
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
-.0187519297793868
|
|
Packit |
67cb25 |
-.0091443848250055,
|
|
Packit |
67cb25 |
.0009010457337825,
|
|
Packit |
67cb25 |
-.0001394184127221,
|
|
Packit |
67cb25 |
.0000273815815785,
|
|
Packit |
67cb25 |
-.0000062750421119,
|
|
Packit |
67cb25 |
.0000016064844184,
|
|
Packit |
67cb25 |
-.0000004476392158,
|
|
Packit |
67cb25 |
.0000001334635874,
|
|
Packit |
67cb25 |
-.0000000420735334,
|
|
Packit |
67cb25 |
.0000000139021990,
|
|
Packit |
67cb25 |
-.0000000047831848,
|
|
Packit |
67cb25 |
.0000000017047897,
|
|
Packit |
67cb25 |
-.0000000006268389,
|
|
Packit |
67cb25 |
.0000000002369824,
|
|
Packit |
67cb25 |
-.0000000000918641,
|
|
Packit |
67cb25 |
.0000000000364278,
|
|
Packit |
67cb25 |
-.0000000000147475,
|
|
Packit |
67cb25 |
.0000000000060851,
|
|
Packit |
67cb25 |
-.0000000000025552,
|
|
Packit |
67cb25 |
.0000000000010906,
|
|
Packit |
67cb25 |
-.0000000000004725,
|
|
Packit |
67cb25 |
.0000000000002076,
|
|
Packit |
67cb25 |
-.0000000000000924,
|
|
Packit |
67cb25 |
.0000000000000417,
|
|
Packit |
67cb25 |
-.0000000000000190,
|
|
Packit |
67cb25 |
.0000000000000087,
|
|
Packit |
67cb25 |
-.0000000000000040,
|
|
Packit |
67cb25 |
.0000000000000019,
|
|
Packit |
67cb25 |
-.0000000000000009,
|
|
Packit |
67cb25 |
.0000000000000004,
|
|
Packit |
67cb25 |
-.0000000000000002,
|
|
Packit |
67cb25 |
.0000000000000001,
|
|
Packit |
67cb25 |
-.0000000000000000
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series aip_cs = {
|
|
Packit |
67cb25 |
data_aip,
|
|
Packit |
67cb25 |
35,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
17
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* chebyshev for Bi(x) asymptotic factor
|
|
Packit |
67cb25 |
based on SLATEC bie()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Series for BIP on the interval 1.25000D-01 to 3.53553D-01
|
|
Packit |
67cb25 |
with weighted error 1.91E-17
|
|
Packit |
67cb25 |
log weighted error 16.72
|
|
Packit |
67cb25 |
significant figures required 15.35
|
|
Packit |
67cb25 |
decimal places required 17.41
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Series for BIP2 on the interval 0. to 1.25000D-01
|
|
Packit |
67cb25 |
with weighted error 1.05E-18
|
|
Packit |
67cb25 |
log weighted error 17.98
|
|
Packit |
67cb25 |
significant figures required 16.74
|
|
Packit |
67cb25 |
decimal places required 18.71
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double data_bip[24] = {
|
|
Packit |
67cb25 |
-0.08322047477943447,
|
|
Packit |
67cb25 |
0.01146118927371174,
|
|
Packit |
67cb25 |
0.00042896440718911,
|
|
Packit |
67cb25 |
-0.00014906639379950,
|
|
Packit |
67cb25 |
-0.00001307659726787,
|
|
Packit |
67cb25 |
0.00000632759839610,
|
|
Packit |
67cb25 |
-0.00000042226696982,
|
|
Packit |
67cb25 |
-0.00000019147186298,
|
|
Packit |
67cb25 |
0.00000006453106284,
|
|
Packit |
67cb25 |
-0.00000000784485467,
|
|
Packit |
67cb25 |
-0.00000000096077216,
|
|
Packit |
67cb25 |
0.00000000070004713,
|
|
Packit |
67cb25 |
-0.00000000017731789,
|
|
Packit |
67cb25 |
0.00000000002272089,
|
|
Packit |
67cb25 |
0.00000000000165404,
|
|
Packit |
67cb25 |
-0.00000000000185171,
|
|
Packit |
67cb25 |
0.00000000000059576,
|
|
Packit |
67cb25 |
-0.00000000000012194,
|
|
Packit |
67cb25 |
0.00000000000001334,
|
|
Packit |
67cb25 |
0.00000000000000172,
|
|
Packit |
67cb25 |
-0.00000000000000145,
|
|
Packit |
67cb25 |
0.00000000000000049,
|
|
Packit |
67cb25 |
-0.00000000000000011,
|
|
Packit |
67cb25 |
0.00000000000000001
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series bip_cs = {
|
|
Packit |
67cb25 |
data_bip,
|
|
Packit |
67cb25 |
23,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
14
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double data_bip2[29] = {
|
|
Packit |
67cb25 |
-0.113596737585988679,
|
|
Packit |
67cb25 |
0.0041381473947881595,
|
|
Packit |
67cb25 |
0.0001353470622119332,
|
|
Packit |
67cb25 |
0.0000104273166530153,
|
|
Packit |
67cb25 |
0.0000013474954767849,
|
|
Packit |
67cb25 |
0.0000001696537405438,
|
|
Packit |
67cb25 |
-0.0000000100965008656,
|
|
Packit |
67cb25 |
-0.0000000167291194937,
|
|
Packit |
67cb25 |
-0.0000000045815364485,
|
|
Packit |
67cb25 |
0.0000000003736681366,
|
|
Packit |
67cb25 |
0.0000000005766930320,
|
|
Packit |
67cb25 |
0.0000000000621812650,
|
|
Packit |
67cb25 |
-0.0000000000632941202,
|
|
Packit |
67cb25 |
-0.0000000000149150479,
|
|
Packit |
67cb25 |
0.0000000000078896213,
|
|
Packit |
67cb25 |
0.0000000000024960513,
|
|
Packit |
67cb25 |
-0.0000000000012130075,
|
|
Packit |
67cb25 |
-0.0000000000003740493,
|
|
Packit |
67cb25 |
0.0000000000002237727,
|
|
Packit |
67cb25 |
0.0000000000000474902,
|
|
Packit |
67cb25 |
-0.0000000000000452616,
|
|
Packit |
67cb25 |
-0.0000000000000030172,
|
|
Packit |
67cb25 |
0.0000000000000091058,
|
|
Packit |
67cb25 |
-0.0000000000000009814,
|
|
Packit |
67cb25 |
-0.0000000000000016429,
|
|
Packit |
67cb25 |
0.0000000000000005533,
|
|
Packit |
67cb25 |
0.0000000000000002175,
|
|
Packit |
67cb25 |
-0.0000000000000001737,
|
|
Packit |
67cb25 |
-0.0000000000000000010
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series bip2_cs = {
|
|
Packit |
67cb25 |
data_bip2,
|
|
Packit |
67cb25 |
28,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
10
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* assumes x >= 1.0 */
|
|
Packit |
67cb25 |
inline static int
|
|
Packit |
67cb25 |
airy_aie(const double x, gsl_mode_t mode, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double sqx = sqrt(x);
|
|
Packit |
67cb25 |
double z = 2.0/(x*sqx) - 1.0;
|
|
Packit |
67cb25 |
double y = sqrt(sqx);
|
|
Packit |
67cb25 |
gsl_sf_result result_c;
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&aip_cs, z, mode, &result_c);
|
|
Packit |
67cb25 |
result->val = (0.28125 + result_c.val)/y;
|
|
Packit |
67cb25 |
result->err = result_c.err/y + GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* assumes x >= 2.0 */
|
|
Packit |
67cb25 |
static int airy_bie(const double x, gsl_mode_t mode, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double ATR = 8.7506905708484345;
|
|
Packit |
67cb25 |
const double BTR = -2.0938363213560543;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < 4.0) {
|
|
Packit |
67cb25 |
double sqx = sqrt(x);
|
|
Packit |
67cb25 |
double z = ATR/(x*sqx) + BTR;
|
|
Packit |
67cb25 |
double y = sqrt(sqx);
|
|
Packit |
67cb25 |
gsl_sf_result result_c;
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&bip_cs, z, mode, &result_c);
|
|
Packit |
67cb25 |
result->val = (0.625 + result_c.val)/y;
|
|
Packit |
67cb25 |
result->err = result_c.err/y + GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
double sqx = sqrt(x);
|
|
Packit |
67cb25 |
double z = 16.0/(x*sqx) - 1.0;
|
|
Packit |
67cb25 |
double y = sqrt(sqx);
|
|
Packit |
67cb25 |
gsl_sf_result result_c;
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&bip2_cs, z, mode, &result_c);
|
|
Packit |
67cb25 |
result->val = (0.625 + result_c.val)/y;
|
|
Packit |
67cb25 |
result->err = result_c.err/y + GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_airy_Ai_e(const double x, const gsl_mode_t mode, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < -1.0) {
|
|
Packit |
67cb25 |
gsl_sf_result mod;
|
|
Packit |
67cb25 |
gsl_sf_result theta;
|
|
Packit |
67cb25 |
gsl_sf_result cos_result;
|
|
Packit |
67cb25 |
int stat_mp = airy_mod_phase(x, mode, &mod, &theta);
|
|
Packit |
67cb25 |
int stat_cos = gsl_sf_cos_err_e(theta.val, theta.err, &cos_result);
|
|
Packit |
67cb25 |
result->val = mod.val * cos_result.val;
|
|
Packit |
67cb25 |
result->err = fabs(mod.val * cos_result.err) + fabs(cos_result.val * mod.err);
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_mp, stat_cos);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 1.0) {
|
|
Packit |
67cb25 |
const double z = x*x*x;
|
|
Packit |
67cb25 |
gsl_sf_result result_c0;
|
|
Packit |
67cb25 |
gsl_sf_result result_c1;
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&aif_cs, z, mode, &result_c0);
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&aig_cs, z, mode, &result_c1);
|
|
Packit |
67cb25 |
result->val = 0.375 + (result_c0.val - x*(0.25 + result_c1.val));
|
|
Packit |
67cb25 |
result->err = result_c0.err + fabs(x*result_c1.err);
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
double x32 = x * sqrt(x);
|
|
Packit |
67cb25 |
double s = exp(-2.0*x32/3.0);
|
|
Packit |
67cb25 |
gsl_sf_result result_aie;
|
|
Packit |
67cb25 |
int stat_aie = airy_aie(x, mode, &result_aie);
|
|
Packit |
67cb25 |
result->val = result_aie.val * s;
|
|
Packit |
67cb25 |
result->err = result_aie.err * s + result->val * x32 * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
CHECK_UNDERFLOW(result);
|
|
Packit |
67cb25 |
return stat_aie;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_airy_Ai_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < -1.0) {
|
|
Packit |
67cb25 |
gsl_sf_result mod;
|
|
Packit |
67cb25 |
gsl_sf_result theta;
|
|
Packit |
67cb25 |
gsl_sf_result cos_result;
|
|
Packit |
67cb25 |
int stat_mp = airy_mod_phase(x, mode, &mod, &theta);
|
|
Packit |
67cb25 |
int stat_cos = gsl_sf_cos_err_e(theta.val, theta.err, &cos_result);
|
|
Packit |
67cb25 |
result->val = mod.val * cos_result.val;
|
|
Packit |
67cb25 |
result->err = fabs(mod.val * cos_result.err) + fabs(cos_result.val * mod.err);
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_mp, stat_cos);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 1.0) {
|
|
Packit |
67cb25 |
const double z = x*x*x;
|
|
Packit |
67cb25 |
gsl_sf_result result_c0;
|
|
Packit |
67cb25 |
gsl_sf_result result_c1;
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&aif_cs, z, mode, &result_c0);
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&aig_cs, z, mode, &result_c1);
|
|
Packit |
67cb25 |
result->val = 0.375 + (result_c0.val - x*(0.25 + result_c1.val));
|
|
Packit |
67cb25 |
result->err = result_c0.err + fabs(x*result_c1.err);
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x > 0.0) {
|
|
Packit |
67cb25 |
const double scale = exp(2.0/3.0 * sqrt(z));
|
|
Packit |
67cb25 |
result->val *= scale;
|
|
Packit |
67cb25 |
result->err *= scale;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
return airy_aie(x, mode, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_airy_Bi_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
if(x < -1.0) {
|
|
Packit |
67cb25 |
gsl_sf_result mod;
|
|
Packit |
67cb25 |
gsl_sf_result theta;
|
|
Packit |
67cb25 |
gsl_sf_result sin_result;
|
|
Packit |
67cb25 |
int stat_mp = airy_mod_phase(x, mode, &mod, &theta);
|
|
Packit |
67cb25 |
int stat_sin = gsl_sf_sin_err_e(theta.val, theta.err, &sin_result);
|
|
Packit |
67cb25 |
result->val = mod.val * sin_result.val;
|
|
Packit |
67cb25 |
result->err = fabs(mod.val * sin_result.err) + fabs(sin_result.val * mod.err);
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_mp, stat_sin);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < 1.0) {
|
|
Packit |
67cb25 |
const double z = x*x*x;
|
|
Packit |
67cb25 |
gsl_sf_result result_c0;
|
|
Packit |
67cb25 |
gsl_sf_result result_c1;
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&bif_cs, z, mode, &result_c0);
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&big_cs, z, mode, &result_c1);
|
|
Packit |
67cb25 |
result->val = 0.625 + result_c0.val + x*(0.4375 + result_c1.val);
|
|
Packit |
67cb25 |
result->err = result_c0.err + fabs(x * result_c1.err);
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 2.0) {
|
|
Packit |
67cb25 |
const double z = (2.0*x*x*x - 9.0)/7.0;
|
|
Packit |
67cb25 |
gsl_sf_result result_c0;
|
|
Packit |
67cb25 |
gsl_sf_result result_c1;
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&bif2_cs, z, mode, &result_c0);
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&big2_cs, z, mode, &result_c1);
|
|
Packit |
67cb25 |
result->val = 1.125 + result_c0.val + x*(0.625 + result_c1.val);
|
|
Packit |
67cb25 |
result->err = result_c0.err + fabs(x * result_c1.err);
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double y = 2.0*x*sqrt(x)/3.0;
|
|
Packit |
67cb25 |
const double s = exp(y);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(y > GSL_LOG_DBL_MAX - 1.0) {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
gsl_sf_result result_bie;
|
|
Packit |
67cb25 |
int stat_bie = airy_bie(x, mode, &result_bie);
|
|
Packit |
67cb25 |
result->val = result_bie.val * s;
|
|
Packit |
67cb25 |
result->err = result_bie.err * s + fabs(1.5*y * (GSL_DBL_EPSILON * result->val));
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_bie;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_airy_Bi_scaled_e(const double x, gsl_mode_t mode, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < -1.0) {
|
|
Packit |
67cb25 |
gsl_sf_result mod;
|
|
Packit |
67cb25 |
gsl_sf_result theta;
|
|
Packit |
67cb25 |
gsl_sf_result sin_result;
|
|
Packit |
67cb25 |
int stat_mp = airy_mod_phase(x, mode, &mod, &theta);
|
|
Packit |
67cb25 |
int stat_sin = gsl_sf_sin_err_e(theta.val, theta.err, &sin_result);
|
|
Packit |
67cb25 |
result->val = mod.val * sin_result.val;
|
|
Packit |
67cb25 |
result->err = fabs(mod.val * sin_result.err) + fabs(sin_result.val * mod.err);
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_mp, stat_sin);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < 1.0) {
|
|
Packit |
67cb25 |
const double z = x*x*x;
|
|
Packit |
67cb25 |
gsl_sf_result result_c0;
|
|
Packit |
67cb25 |
gsl_sf_result result_c1;
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&bif_cs, z, mode, &result_c0);
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&big_cs, z, mode, &result_c1);
|
|
Packit |
67cb25 |
result->val = 0.625 + result_c0.val + x*(0.4375 + result_c1.val);
|
|
Packit |
67cb25 |
result->err = result_c0.err + fabs(x * result_c1.err);
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
if(x > 0.0) {
|
|
Packit |
67cb25 |
const double scale = exp(-2.0/3.0 * sqrt(z));
|
|
Packit |
67cb25 |
result->val *= scale;
|
|
Packit |
67cb25 |
result->err *= scale;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 2.0) {
|
|
Packit |
67cb25 |
const double x3 = x*x*x;
|
|
Packit |
67cb25 |
const double z = (2.0*x3 - 9.0)/7.0;
|
|
Packit |
67cb25 |
const double s = exp(-2.0/3.0 * sqrt(x3));
|
|
Packit |
67cb25 |
gsl_sf_result result_c0;
|
|
Packit |
67cb25 |
gsl_sf_result result_c1;
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&bif2_cs, z, mode, &result_c0);
|
|
Packit |
67cb25 |
cheb_eval_mode_e(&big2_cs, z, mode, &result_c1);
|
|
Packit |
67cb25 |
result->val = s * (1.125 + result_c0.val + x*(0.625 + result_c1.val));
|
|
Packit |
67cb25 |
result->err = s * (result_c0.err + fabs(x * result_c1.err));
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
return airy_bie(x, mode, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "eval.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_airy_Ai(const double x, gsl_mode_t mode)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_airy_Ai_e(x, mode, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_airy_Ai_scaled(const double x, gsl_mode_t mode)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_airy_Ai_scaled_e(x, mode, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_airy_Bi(const double x, gsl_mode_t mode)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_airy_Bi_e(x, mode, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_airy_Bi_scaled(const double x, gsl_mode_t mode)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_airy_Bi_scaled_e(x, mode, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|