| /* specfunc/transport.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_transport.h> |
| |
| #include "error.h" |
| #include "check.h" |
| |
| #include "chebyshev.h" |
| #include "cheb_eval.c" |
| |
| static double transport2_data[18] = { |
| 1.671760446434538503, |
| -0.147735359946794490, |
| 0.148213819946936338e-01, |
| -0.14195330326305613e-02, |
| 0.1306541324415708e-03, |
| -0.117155795867579e-04, |
| 0.10333498445756e-05, |
| -0.901911304223e-07, |
| 0.78177169833e-08, |
| -0.6744565684e-09, |
| 0.579946394e-10, |
| -0.49747619e-11, |
| 0.425961e-12, |
| -0.36422e-13, |
| 0.3111e-14, |
| -0.265e-15, |
| 0.23e-16, |
| -0.19e-17 |
| }; |
| static cheb_series transport2_cs = { |
| transport2_data, |
| 17, |
| -1, 1, |
| 9 |
| }; |
| |
| static double transport3_data[18] = { |
| 0.762012543243872007, |
| -0.105674387705058533, |
| 0.119778084819657810e-01, |
| -0.12144015203698307e-02, |
| 0.1155099769392855e-03, |
| -0.105815992124423e-04, |
| 0.9474663385302e-06, |
| -0.836221212858e-07, |
| 0.73109099278e-08, |
| -0.6350594779e-09, |
| 0.549118282e-10, |
| -0.47321395e-11, |
| 0.4067695e-12, |
| -0.348971e-13, |
| 0.29892e-14, |
| -0.256e-15, |
| 0.219e-16, |
| -0.19e-17 |
| }; |
| static cheb_series transport3_cs = { |
| transport3_data, |
| 17, |
| -1, 1, |
| 9 |
| }; |
| |
| |
| static double transport4_data[18] = { |
| 0.4807570994615110579, |
| -0.8175378810321083956e-01, |
| 0.1002700665975162973e-01, |
| -0.10599339359820151e-02, |
| 0.1034506245030405e-03, |
| -0.96442705485899e-05, |
| 0.8745544408515e-06, |
| -0.779321207981e-07, |
| 0.68649886141e-08, |
| -0.5999571076e-09, |
| 0.521366241e-10, |
| -0.45118382e-11, |
| 0.3892159e-12, |
| -0.334936e-13, |
| 0.28767e-14, |
| -0.2467e-15, |
| 0.211e-16, |
| -0.18e-17 |
| }; |
| static cheb_series transport4_cs = { |
| transport4_data, |
| 17, |
| -1, 1, |
| 9 |
| }; |
| |
| |
| static double transport5_data[18] = { |
| 0.347777777133910789, |
| -0.66456988976050428e-01, |
| 0.8611072656883309e-02, |
| -0.9396682223755538e-03, |
| 0.936324806081513e-04, |
| -0.88571319340833e-05, |
| 0.811914989145e-06, |
| -0.72957654233e-07, |
| 0.646971455e-08, |
| -0.568490283e-09, |
| 0.49625598e-10, |
| -0.4310940e-11, |
| 0.373100e-12, |
| -0.32198e-13, |
| 0.2772e-14, |
| -0.238e-15, |
| 0.21e-16, |
| -0.18e-17 |
| }; |
| static cheb_series transport5_cs = { |
| transport5_data, |
| 17, |
| -1, 1, |
| 9 |
| }; |
| |
| |
| static |
| double |
| transport_sumexp(const int numexp, const int order, const double t, double x) |
| { |
| double rk = (double)numexp; |
| double sumexp = 0.0; |
| int k; |
| for(k=1; k<=numexp; k++) { |
| double sum2 = 1.0; |
| double xk = 1.0/(rk*x); |
| double xk1 = 1.0; |
| int j; |
| for(j=1; j<=order; j++) { |
| sum2 = sum2*xk1*xk + 1.0; |
| xk1 += 1.0; |
| } |
| sumexp *= t; |
| sumexp += sum2; |
| rk -= 1.0; |
| } |
| return sumexp; |
| } |
| |
| |
| /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ |
| |
| int |
| gsl_sf_transport_2_e(const double x, gsl_sf_result * result) |
| { |
| const double val_infinity = 3.289868133696452873; |
| |
| /* CHECK_POINTER(result) */ |
| |
| if(x < 0.0) { |
| DOMAIN_ERROR(result); |
| } |
| else if(x < 3.0*GSL_SQRT_DBL_EPSILON) { |
| result->val = x; |
| result->err = GSL_DBL_EPSILON*fabs(x) + x*x/2.0; |
| return GSL_SUCCESS; |
| } |
| else if(x <= 4.0) { |
| double t = (x*x/8.0 - 0.5) - 0.5; |
| gsl_sf_result result_c; |
| cheb_eval_e(&transport2_cs, t, &result_c); |
| result->val = x * result_c.val; |
| result->err = x * result_c.err; |
| result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
| return GSL_SUCCESS; |
| } |
| else if(x < -GSL_LOG_DBL_EPSILON) { |
| const int numexp = (int)((-GSL_LOG_DBL_EPSILON)/x) + 1; |
| const double sumexp = transport_sumexp(numexp, 2, exp(-x), x); |
| const double t = 2.0 * log(x) - x + log(sumexp); |
| if(t < GSL_LOG_DBL_EPSILON) { |
| result->val = val_infinity; |
| result->err = 2.0 * GSL_DBL_EPSILON * val_infinity; |
| } |
| else { |
| const double et = exp(t); |
| result->val = val_infinity - et; |
| result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + fabs(t) * et); |
| } |
| return GSL_SUCCESS; |
| } |
| else if(x < 2.0/GSL_DBL_EPSILON) { |
| const int numexp = 1; |
| const double sumexp = transport_sumexp(numexp, 2, 1.0, x); |
| const double t = 2.0 * log(x) - x + log(sumexp); |
| if(t < GSL_LOG_DBL_EPSILON) { |
| result->val = val_infinity; |
| result->err = 2.0 * GSL_DBL_EPSILON * val_infinity; |
| } |
| else { |
| const double et = exp(t); |
| result->val = val_infinity - et; |
| result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et); |
| } |
| return GSL_SUCCESS; |
| } |
| else { |
| const double t = 2.0 * log(x) - x; |
| if(t < GSL_LOG_DBL_EPSILON) { |
| result->val = val_infinity; |
| result->err = 2.0 * GSL_DBL_EPSILON * val_infinity; |
| } |
| else { |
| const double et = exp(t); |
| result->val = val_infinity - et; |
| result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et); |
| } |
| return GSL_SUCCESS; |
| } |
| } |
| |
| |
| int |
| gsl_sf_transport_3_e(const double x, gsl_sf_result * result) |
| { |
| const double val_infinity = 7.212341418957565712; |
| |
| /* CHECK_POINTER(result) */ |
| |
| if(x < 0.0) { |
| DOMAIN_ERROR(result); |
| } |
| else if(x == 0.0) { |
| result->val = 0.0; |
| result->err = 0.0; |
| return GSL_SUCCESS; |
| } |
| else if(x < 3.0*GSL_SQRT_DBL_EPSILON) { |
| result->val = 0.5*x*x; |
| result->err = 2.0 * GSL_DBL_EPSILON * result->val; |
| CHECK_UNDERFLOW(result); |
| return GSL_SUCCESS; |
| } |
| else if(x <= 4.0) { |
| const double x2 = x*x; |
| const double t = (x2/8.0 - 0.5) - 0.5; |
| gsl_sf_result result_c; |
| cheb_eval_e(&transport3_cs, t, &result_c); |
| result->val = x2 * result_c.val; |
| result->err = x2 * result_c.err; |
| result->err += GSL_DBL_EPSILON * fabs(result->val); |
| return GSL_SUCCESS; |
| } |
| else if(x < -GSL_LOG_DBL_EPSILON) { |
| const int numexp = (int)((-GSL_LOG_DBL_EPSILON)/x) + 1; |
| const double sumexp = transport_sumexp(numexp, 3, exp(-x), x); |
| const double t = 3.0 * log(x) - x + log(sumexp); |
| if(t < GSL_LOG_DBL_EPSILON) { |
| result->val = val_infinity; |
| result->err = 2.0 * GSL_DBL_EPSILON * val_infinity; |
| } |
| else { |
| const double et = exp(t); |
| result->val = val_infinity - et; |
| result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + fabs(t) * et); |
| } |
| return GSL_SUCCESS; |
| } |
| else if(x < 3.0/GSL_DBL_EPSILON) { |
| const int numexp = 1; |
| const double sumexp = transport_sumexp(numexp, 3, 1.0, x); |
| const double t = 3.0 * log(x) - x + log(sumexp); |
| if(t < GSL_LOG_DBL_EPSILON) { |
| result->val = val_infinity; |
| result->err = 2.0 * GSL_DBL_EPSILON * val_infinity; |
| } |
| else { |
| const double et = exp(t); |
| result->val = val_infinity - et; |
| result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et); |
| } |
| return GSL_SUCCESS; |
| } |
| else { |
| const double t = 3.0 * log(x) - x; |
| if(t < GSL_LOG_DBL_EPSILON) { |
| result->val = val_infinity; |
| result->err = 2.0 * GSL_DBL_EPSILON * val_infinity; |
| } |
| else { |
| const double et = exp(t); |
| result->val = val_infinity - et; |
| result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et); |
| } |
| return GSL_SUCCESS; |
| } |
| } |
| |
| |
| int |
| gsl_sf_transport_4_e(const double x, gsl_sf_result * result) |
| { |
| const double val_infinity = 25.97575760906731660; |
| |
| /* CHECK_POINTER(result) */ |
| |
| if(x < 0.0) { |
| DOMAIN_ERROR(result); |
| } |
| else if(x == 0.0) { |
| result->val = 0.0; |
| result->err = 0.0; |
| return GSL_SUCCESS; |
| } |
| else if(x < 3.0*GSL_SQRT_DBL_EPSILON) { |
| result->val = x*x*x/3.0; |
| result->err = 3.0 * GSL_DBL_EPSILON * result->val; |
| CHECK_UNDERFLOW(result); |
| return GSL_SUCCESS; |
| } |
| else if(x <= 4.0) { |
| const double x2 = x*x; |
| const double t = (x2/8.0 - 0.5) - 0.5; |
| gsl_sf_result result_c; |
| cheb_eval_e(&transport4_cs, t, &result_c); |
| result->val = x2*x * result_c.val; |
| result->err = x2*x * result_c.err; |
| result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
| return GSL_SUCCESS; |
| } |
| else if(x < -GSL_LOG_DBL_EPSILON) { |
| const int numexp = (int)((-GSL_LOG_DBL_EPSILON)/x) + 1; |
| const double sumexp = transport_sumexp(numexp, 4, exp(-x), x); |
| const double t = 4.0 * log(x) - x + log(sumexp); |
| if(t < GSL_LOG_DBL_EPSILON) { |
| result->val = val_infinity; |
| result->err = 2.0 * GSL_DBL_EPSILON * val_infinity; |
| } |
| else { |
| const double et = exp(t); |
| result->val = val_infinity - et; |
| result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et); |
| } |
| return GSL_SUCCESS; |
| } |
| else if(x < 3.0/GSL_DBL_EPSILON) { |
| const int numexp = 1; |
| const double sumexp = transport_sumexp(numexp, 4, 1.0, x); |
| const double t = 4.0 * log(x) - x + log(sumexp); |
| if(t < GSL_LOG_DBL_EPSILON) { |
| result->val = val_infinity; |
| result->err = 2.0 * GSL_DBL_EPSILON * val_infinity; |
| } |
| else { |
| const double et = exp(t); |
| result->val = val_infinity - et; |
| result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et); |
| } |
| return GSL_SUCCESS; |
| } |
| else { |
| const double t = 4.0 * log(x) - x; |
| if(t < GSL_LOG_DBL_EPSILON) { |
| result->val = val_infinity; |
| result->err = 2.0 * GSL_DBL_EPSILON * val_infinity; |
| } |
| else { |
| const double et = exp(t); |
| result->val = val_infinity - et; |
| result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et); |
| } |
| return GSL_SUCCESS; |
| } |
| } |
| |
| |
| int |
| gsl_sf_transport_5_e(const double x, gsl_sf_result * result) |
| { |
| const double val_infinity = 124.4313306172043912; |
| |
| /* CHECK_POINTER(result) */ |
| |
| if(x < 0.0) { |
| DOMAIN_ERROR(result); |
| } |
| else if(x == 0.0) { |
| result->val = 0.0; |
| result->err = 0.0; |
| return GSL_SUCCESS; |
| } |
| else if(x < 3.0*GSL_SQRT_DBL_EPSILON) { |
| result->val = x*x*x*x/4.0; |
| result->err = 4.0 * GSL_DBL_EPSILON * result->val; |
| CHECK_UNDERFLOW(result); |
| return GSL_SUCCESS; |
| } |
| else if(x <= 4.0) { |
| const double x2 = x*x; |
| const double t = (x2/8.0 - 0.5) - 0.5; |
| gsl_sf_result result_c; |
| cheb_eval_e(&transport5_cs, t, &result_c); |
| result->val = x2*x2 * result_c.val; |
| result->err = x2*x2 * result_c.err; |
| result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
| return GSL_SUCCESS; |
| } |
| else if(x < -GSL_LOG_DBL_EPSILON) { |
| const int numexp = (int)((-GSL_LOG_DBL_EPSILON)/x) + 1; |
| const double sumexp = transport_sumexp(numexp, 5, exp(-x), x); |
| const double t = 5.0 * log(x) - x + log(sumexp); |
| if(t < GSL_LOG_DBL_EPSILON) { |
| result->val = val_infinity; |
| result->err = 2.0 * GSL_DBL_EPSILON * val_infinity; |
| } |
| else { |
| const double et = exp(t); |
| result->val = val_infinity - et; |
| result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et); |
| } |
| return GSL_SUCCESS; |
| } |
| else if(x < 3.0/GSL_DBL_EPSILON) { |
| const int numexp = 1; |
| const double sumexp = transport_sumexp(numexp, 5, 1.0, x); |
| const double t = 5.0 * log(x) - x + log(sumexp); |
| if(t < GSL_LOG_DBL_EPSILON) { |
| result->val = val_infinity; |
| result->err = 2.0 * GSL_DBL_EPSILON * val_infinity; |
| } |
| else { |
| const double et = exp(t); |
| result->val = val_infinity - et; |
| result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et); |
| } |
| return GSL_SUCCESS; |
| } |
| else { |
| const double t = 5.0 * log(x) - x; |
| if(t < GSL_LOG_DBL_EPSILON) { |
| result->val = val_infinity; |
| result->err = 2.0 * GSL_DBL_EPSILON * val_infinity; |
| } |
| else { |
| const double et = exp(t); |
| result->val = val_infinity - et; |
| result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et); |
| } |
| return GSL_SUCCESS; |
| } |
| } |
| |
| /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ |
| |
| #include "eval.h" |
| |
| double gsl_sf_transport_2(const double x) |
| { |
| EVAL_RESULT(gsl_sf_transport_2_e(x, &result)); |
| } |
| |
| double gsl_sf_transport_3(const double x) |
| { |
| EVAL_RESULT(gsl_sf_transport_3_e(x, &result)); |
| } |
| |
| double gsl_sf_transport_4(const double x) |
| { |
| EVAL_RESULT(gsl_sf_transport_4_e(x, &result)); |
| } |
| |
| double gsl_sf_transport_5(const double x) |
| { |
| EVAL_RESULT(gsl_sf_transport_5_e(x, &result)); |
| } |
