|
Packit |
67cb25 |
/* specfunc/trig.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_log.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_trig.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "error.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "chebyshev.h"
|
|
Packit |
67cb25 |
#include "cheb_eval.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* sinh(x) series
|
|
Packit |
67cb25 |
* double-precision for |x| < 1.0
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
inline
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
sinh_series(const double x, double * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double y = x*x;
|
|
Packit |
67cb25 |
const double c0 = 1.0/6.0;
|
|
Packit |
67cb25 |
const double c1 = 1.0/120.0;
|
|
Packit |
67cb25 |
const double c2 = 1.0/5040.0;
|
|
Packit |
67cb25 |
const double c3 = 1.0/362880.0;
|
|
Packit |
67cb25 |
const double c4 = 1.0/39916800.0;
|
|
Packit |
67cb25 |
const double c5 = 1.0/6227020800.0;
|
|
Packit |
67cb25 |
const double c6 = 1.0/1307674368000.0;
|
|
Packit |
67cb25 |
const double c7 = 1.0/355687428096000.0;
|
|
Packit |
67cb25 |
*result = x*(1.0 + y*(c0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*(c6+y*c7))))))));
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* cosh(x)-1 series
|
|
Packit |
67cb25 |
* double-precision for |x| < 1.0
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
inline
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
cosh_m1_series(const double x, double * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double y = x*x;
|
|
Packit |
67cb25 |
const double c0 = 0.5;
|
|
Packit |
67cb25 |
const double c1 = 1.0/24.0;
|
|
Packit |
67cb25 |
const double c2 = 1.0/720.0;
|
|
Packit |
67cb25 |
const double c3 = 1.0/40320.0;
|
|
Packit |
67cb25 |
const double c4 = 1.0/3628800.0;
|
|
Packit |
67cb25 |
const double c5 = 1.0/479001600.0;
|
|
Packit |
67cb25 |
const double c6 = 1.0/87178291200.0;
|
|
Packit |
67cb25 |
const double c7 = 1.0/20922789888000.0;
|
|
Packit |
67cb25 |
const double c8 = 1.0/6402373705728000.0;
|
|
Packit |
67cb25 |
*result = y*(c0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*(c6+y*(c7+y*c8))))))));
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Chebyshev expansion for f(t) = sinc((t+1)/2), -1 < t < 1
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double sinc_data[17] = {
|
|
Packit |
67cb25 |
1.133648177811747875422,
|
|
Packit |
67cb25 |
-0.532677564732557348781,
|
|
Packit |
67cb25 |
-0.068293048346633177859,
|
|
Packit |
67cb25 |
0.033403684226353715020,
|
|
Packit |
67cb25 |
0.001485679893925747818,
|
|
Packit |
67cb25 |
-0.000734421305768455295,
|
|
Packit |
67cb25 |
-0.000016837282388837229,
|
|
Packit |
67cb25 |
0.000008359950146618018,
|
|
Packit |
67cb25 |
0.000000117382095601192,
|
|
Packit |
67cb25 |
-0.000000058413665922724,
|
|
Packit |
67cb25 |
-0.000000000554763755743,
|
|
Packit |
67cb25 |
0.000000000276434190426,
|
|
Packit |
67cb25 |
0.000000000001895374892,
|
|
Packit |
67cb25 |
-0.000000000000945237101,
|
|
Packit |
67cb25 |
-0.000000000000004900690,
|
|
Packit |
67cb25 |
0.000000000000002445383,
|
|
Packit |
67cb25 |
0.000000000000000009925
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series sinc_cs = {
|
|
Packit |
67cb25 |
sinc_data,
|
|
Packit |
67cb25 |
16,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
10
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Chebyshev expansion for f(t) = g((t+1)Pi/8), -1
|
|
Packit |
67cb25 |
* g(x) = (sin(x)/x - 1)/(x*x)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double sin_data[12] = {
|
|
Packit |
67cb25 |
-0.3295190160663511504173,
|
|
Packit |
67cb25 |
0.0025374284671667991990,
|
|
Packit |
67cb25 |
0.0006261928782647355874,
|
|
Packit |
67cb25 |
-4.6495547521854042157541e-06,
|
|
Packit |
67cb25 |
-5.6917531549379706526677e-07,
|
|
Packit |
67cb25 |
3.7283335140973803627866e-09,
|
|
Packit |
67cb25 |
3.0267376484747473727186e-10,
|
|
Packit |
67cb25 |
-1.7400875016436622322022e-12,
|
|
Packit |
67cb25 |
-1.0554678305790849834462e-13,
|
|
Packit |
67cb25 |
5.3701981409132410797062e-16,
|
|
Packit |
67cb25 |
2.5984137983099020336115e-17,
|
|
Packit |
67cb25 |
-1.1821555255364833468288e-19
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series sin_cs = {
|
|
Packit |
67cb25 |
sin_data,
|
|
Packit |
67cb25 |
11,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
11
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Chebyshev expansion for f(t) = g((t+1)Pi/8), -1
|
|
Packit |
67cb25 |
* g(x) = (2(cos(x) - 1)/(x^2) + 1) / x^2
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double cos_data[11] = {
|
|
Packit |
67cb25 |
0.165391825637921473505668118136,
|
|
Packit |
67cb25 |
-0.00084852883845000173671196530195,
|
|
Packit |
67cb25 |
-0.000210086507222940730213625768083,
|
|
Packit |
67cb25 |
1.16582269619760204299639757584e-6,
|
|
Packit |
67cb25 |
1.43319375856259870334412701165e-7,
|
|
Packit |
67cb25 |
-7.4770883429007141617951330184e-10,
|
|
Packit |
67cb25 |
-6.0969994944584252706997438007e-11,
|
|
Packit |
67cb25 |
2.90748249201909353949854872638e-13,
|
|
Packit |
67cb25 |
1.77126739876261435667156490461e-14,
|
|
Packit |
67cb25 |
-7.6896421502815579078577263149e-17,
|
|
Packit |
67cb25 |
-3.7363121133079412079201377318e-18
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series cos_cs = {
|
|
Packit |
67cb25 |
cos_data,
|
|
Packit |
67cb25 |
10,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
10
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* I would have prefered just using the library sin() function.
|
|
Packit |
67cb25 |
* But after some experimentation I decided that there was
|
|
Packit |
67cb25 |
* no good way to understand the error; library sin() is just a black box.
|
|
Packit |
67cb25 |
* So we have to roll our own.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_sin_e(double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double P1 = 7.85398125648498535156e-1;
|
|
Packit |
67cb25 |
const double P2 = 3.77489470793079817668e-8;
|
|
Packit |
67cb25 |
const double P3 = 2.69515142907905952645e-15;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double sgn_x = GSL_SIGN(x);
|
|
Packit |
67cb25 |
const double abs_x = fabs(x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(abs_x < GSL_ROOT4_DBL_EPSILON) {
|
|
Packit |
67cb25 |
const double x2 = x*x;
|
|
Packit |
67cb25 |
result->val = x * (1.0 - x2/6.0);
|
|
Packit |
67cb25 |
result->err = fabs(x*x2*x2 / 100.0);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
double sgn_result = sgn_x;
|
|
Packit |
67cb25 |
double y = floor(abs_x/(0.25*M_PI));
|
|
Packit |
67cb25 |
int octant = y - ldexp(floor(ldexp(y,-3)),3);
|
|
Packit |
67cb25 |
int stat_cs;
|
|
Packit |
67cb25 |
double z;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(GSL_IS_ODD(octant)) {
|
|
Packit |
67cb25 |
octant += 1;
|
|
Packit |
67cb25 |
octant &= 07;
|
|
Packit |
67cb25 |
y += 1.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(octant > 3) {
|
|
Packit |
67cb25 |
octant -= 4;
|
|
Packit |
67cb25 |
sgn_result = -sgn_result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
z = ((abs_x - y * P1) - y * P2) - y * P3;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(octant == 0) {
|
|
Packit |
67cb25 |
gsl_sf_result sin_cs_result;
|
|
Packit |
67cb25 |
const double t = 8.0*fabs(z)/M_PI - 1.0;
|
|
Packit |
67cb25 |
stat_cs = cheb_eval_e(&sin_cs, t, &sin_cs_result);
|
|
Packit |
67cb25 |
result->val = z * (1.0 + z*z * sin_cs_result.val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else { /* octant == 2 */
|
|
Packit |
67cb25 |
gsl_sf_result cos_cs_result;
|
|
Packit |
67cb25 |
const double t = 8.0*fabs(z)/M_PI - 1.0;
|
|
Packit |
67cb25 |
stat_cs = cheb_eval_e(&cos_cs, t, &cos_cs_result);
|
|
Packit |
67cb25 |
result->val = 1.0 - 0.5*z*z * (1.0 - z*z * cos_cs_result.val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val *= sgn_result;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(abs_x > 1.0/GSL_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->err = fabs(result->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(abs_x > 100.0/GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->err = 2.0 * abs_x * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(abs_x > 0.1/GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_SQRT_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return stat_cs;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_cos_e(double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double P1 = 7.85398125648498535156e-1;
|
|
Packit |
67cb25 |
const double P2 = 3.77489470793079817668e-8;
|
|
Packit |
67cb25 |
const double P3 = 2.69515142907905952645e-15;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double abs_x = fabs(x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(abs_x < GSL_ROOT4_DBL_EPSILON) {
|
|
Packit |
67cb25 |
const double x2 = x*x;
|
|
Packit |
67cb25 |
result->val = 1.0 - 0.5*x2;
|
|
Packit |
67cb25 |
result->err = fabs(x2*x2/12.0);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
double sgn_result = 1.0;
|
|
Packit |
67cb25 |
double y = floor(abs_x/(0.25*M_PI));
|
|
Packit |
67cb25 |
int octant = y - ldexp(floor(ldexp(y,-3)),3);
|
|
Packit |
67cb25 |
int stat_cs;
|
|
Packit |
67cb25 |
double z;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(GSL_IS_ODD(octant)) {
|
|
Packit |
67cb25 |
octant += 1;
|
|
Packit |
67cb25 |
octant &= 07;
|
|
Packit |
67cb25 |
y += 1.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(octant > 3) {
|
|
Packit |
67cb25 |
octant -= 4;
|
|
Packit |
67cb25 |
sgn_result = -sgn_result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(octant > 1) {
|
|
Packit |
67cb25 |
sgn_result = -sgn_result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
z = ((abs_x - y * P1) - y * P2) - y * P3;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(octant == 0) {
|
|
Packit |
67cb25 |
gsl_sf_result cos_cs_result;
|
|
Packit |
67cb25 |
const double t = 8.0*fabs(z)/M_PI - 1.0;
|
|
Packit |
67cb25 |
stat_cs = cheb_eval_e(&cos_cs, t, &cos_cs_result);
|
|
Packit |
67cb25 |
result->val = 1.0 - 0.5*z*z * (1.0 - z*z * cos_cs_result.val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else { /* octant == 2 */
|
|
Packit |
67cb25 |
gsl_sf_result sin_cs_result;
|
|
Packit |
67cb25 |
const double t = 8.0*fabs(z)/M_PI - 1.0;
|
|
Packit |
67cb25 |
stat_cs = cheb_eval_e(&sin_cs, t, &sin_cs_result);
|
|
Packit |
67cb25 |
result->val = z * (1.0 + z*z * sin_cs_result.val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val *= sgn_result;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(abs_x > 1.0/GSL_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->err = fabs(result->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(abs_x > 100.0/GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->err = 2.0 * abs_x * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(abs_x > 0.1/GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_SQRT_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return stat_cs;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_hypot_e(const double x, const double y, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x == 0.0 && y == 0.0) {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double a = fabs(x);
|
|
Packit |
67cb25 |
const double b = fabs(y);
|
|
Packit |
67cb25 |
const double min = GSL_MIN_DBL(a,b);
|
|
Packit |
67cb25 |
const double max = GSL_MAX_DBL(a,b);
|
|
Packit |
67cb25 |
const double rat = min/max;
|
|
Packit |
67cb25 |
const double root_term = sqrt(1.0 + rat*rat);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(max < GSL_DBL_MAX/root_term) {
|
|
Packit |
67cb25 |
result->val = max * root_term;
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_complex_sin_e(const double zr, const double zi,
|
|
Packit |
67cb25 |
gsl_sf_result * szr, gsl_sf_result * szi)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(szr) */
|
|
Packit |
67cb25 |
/* CHECK_POINTER(szi) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(zi) < 1.0) {
|
|
Packit |
67cb25 |
double ch_m1, sh;
|
|
Packit |
67cb25 |
sinh_series(zi, &sh);
|
|
Packit |
67cb25 |
cosh_m1_series(zi, &ch_m1);
|
|
Packit |
67cb25 |
szr->val = sin(zr)*(ch_m1 + 1.0);
|
|
Packit |
67cb25 |
szi->val = cos(zr)*sh;
|
|
Packit |
67cb25 |
szr->err = 2.0 * GSL_DBL_EPSILON * fabs(szr->val);
|
|
Packit |
67cb25 |
szi->err = 2.0 * GSL_DBL_EPSILON * fabs(szi->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(zi) < GSL_LOG_DBL_MAX) {
|
|
Packit |
67cb25 |
double ex = exp(zi);
|
|
Packit |
67cb25 |
double ch = 0.5*(ex+1.0/ex);
|
|
Packit |
67cb25 |
double sh = 0.5*(ex-1.0/ex);
|
|
Packit |
67cb25 |
szr->val = sin(zr)*ch;
|
|
Packit |
67cb25 |
szi->val = cos(zr)*sh;
|
|
Packit |
67cb25 |
szr->err = 2.0 * GSL_DBL_EPSILON * fabs(szr->val);
|
|
Packit |
67cb25 |
szi->err = 2.0 * GSL_DBL_EPSILON * fabs(szi->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR_2(szr, szi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_complex_cos_e(const double zr, const double zi,
|
|
Packit |
67cb25 |
gsl_sf_result * czr, gsl_sf_result * czi)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(czr) */
|
|
Packit |
67cb25 |
/* CHECK_POINTER(czi) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(zi) < 1.0) {
|
|
Packit |
67cb25 |
double ch_m1, sh;
|
|
Packit |
67cb25 |
sinh_series(zi, &sh);
|
|
Packit |
67cb25 |
cosh_m1_series(zi, &ch_m1);
|
|
Packit |
67cb25 |
czr->val = cos(zr)*(ch_m1 + 1.0);
|
|
Packit |
67cb25 |
czi->val = -sin(zr)*sh;
|
|
Packit |
67cb25 |
czr->err = 2.0 * GSL_DBL_EPSILON * fabs(czr->val);
|
|
Packit |
67cb25 |
czi->err = 2.0 * GSL_DBL_EPSILON * fabs(czi->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(zi) < GSL_LOG_DBL_MAX) {
|
|
Packit |
67cb25 |
double ex = exp(zi);
|
|
Packit |
67cb25 |
double ch = 0.5*(ex+1.0/ex);
|
|
Packit |
67cb25 |
double sh = 0.5*(ex-1.0/ex);
|
|
Packit |
67cb25 |
czr->val = cos(zr)*ch;
|
|
Packit |
67cb25 |
czi->val = -sin(zr)*sh;
|
|
Packit |
67cb25 |
czr->err = 2.0 * GSL_DBL_EPSILON * fabs(czr->val);
|
|
Packit |
67cb25 |
czi->err = 2.0 * GSL_DBL_EPSILON * fabs(czi->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR_2(czr,czi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_complex_logsin_e(const double zr, const double zi,
|
|
Packit |
67cb25 |
gsl_sf_result * lszr, gsl_sf_result * lszi)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(lszr) */
|
|
Packit |
67cb25 |
/* CHECK_POINTER(lszi) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(zi > 60.0) {
|
|
Packit |
67cb25 |
lszr->val = -M_LN2 + zi;
|
|
Packit |
67cb25 |
lszi->val = 0.5*M_PI - zr;
|
|
Packit |
67cb25 |
lszr->err = 2.0 * GSL_DBL_EPSILON * fabs(lszr->val);
|
|
Packit |
67cb25 |
lszi->err = 2.0 * GSL_DBL_EPSILON * fabs(lszi->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(zi < -60.0) {
|
|
Packit |
67cb25 |
lszr->val = -M_LN2 - zi;
|
|
Packit |
67cb25 |
lszi->val = -0.5*M_PI + zr;
|
|
Packit |
67cb25 |
lszr->err = 2.0 * GSL_DBL_EPSILON * fabs(lszr->val);
|
|
Packit |
67cb25 |
lszi->err = 2.0 * GSL_DBL_EPSILON * fabs(lszi->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
gsl_sf_result sin_r, sin_i;
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
gsl_sf_complex_sin_e(zr, zi, &sin_r, &sin_i); /* ok by construction */
|
|
Packit |
67cb25 |
status = gsl_sf_complex_log_e(sin_r.val, sin_i.val, lszr, lszi);
|
|
Packit |
67cb25 |
if(status == GSL_EDOM) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR_2(lszr, lszi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return gsl_sf_angle_restrict_symm_e(&(lszi->val));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_lnsinh_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x <= 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(x) < 1.0) {
|
|
Packit |
67cb25 |
double eps;
|
|
Packit |
67cb25 |
sinh_series(x, &eps);
|
|
Packit |
67cb25 |
result->val = log(eps);
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -0.5*GSL_LOG_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->val = x + log(0.5*(1.0 - exp(-2.0*x)));
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = -M_LN2 + x;
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_lncosh_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(x) < 1.0) {
|
|
Packit |
67cb25 |
double eps;
|
|
Packit |
67cb25 |
cosh_m1_series(x, &eps);
|
|
Packit |
67cb25 |
return gsl_sf_log_1plusx_e(eps, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(x) < -0.5*GSL_LOG_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->val = fabs(x) + log(0.5*(1.0 + exp(-2.0*fabs(x))));
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = -M_LN2 + fabs(x);
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
inline int gsl_sf_sincos_e(const double theta, double * s, double * c)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double tan_half = tan(0.5 * theta);
|
|
Packit |
67cb25 |
double den = 1. + tan_half*tan_half;
|
|
Packit |
67cb25 |
double cos_theta = (1.0 - tan_half*tan_half) / den;
|
|
Packit |
67cb25 |
double sin_theta = 2.0 * tan_half / den;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_polar_to_rect(const double r, const double theta,
|
|
Packit |
67cb25 |
gsl_sf_result * x, gsl_sf_result * y)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double t = theta;
|
|
Packit |
67cb25 |
int status = gsl_sf_angle_restrict_symm_e(&t);
|
|
Packit |
67cb25 |
double c = cos(t);
|
|
Packit |
67cb25 |
double s = sin(t);
|
|
Packit |
67cb25 |
x->val = r * cos(t);
|
|
Packit |
67cb25 |
y->val = r * sin(t);
|
|
Packit |
67cb25 |
x->err = r * fabs(s * GSL_DBL_EPSILON * t);
|
|
Packit |
67cb25 |
x->err += 2.0 * GSL_DBL_EPSILON * fabs(x->val);
|
|
Packit |
67cb25 |
y->err = r * fabs(c * GSL_DBL_EPSILON * t);
|
|
Packit |
67cb25 |
y->err += 2.0 * GSL_DBL_EPSILON * fabs(y->val);
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_rect_to_polar(const double x, const double y,
|
|
Packit |
67cb25 |
gsl_sf_result * r, gsl_sf_result * theta)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int stat_h = gsl_sf_hypot_e(x, y, r);
|
|
Packit |
67cb25 |
if(r->val > 0.0) {
|
|
Packit |
67cb25 |
theta->val = atan2(y, x);
|
|
Packit |
67cb25 |
theta->err = 2.0 * GSL_DBL_EPSILON * fabs(theta->val);
|
|
Packit |
67cb25 |
return stat_h;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(theta);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_angle_restrict_symm_err_e(const double theta, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* synthetic extended precision constants */
|
|
Packit |
67cb25 |
const double P1 = 4 * 7.8539812564849853515625e-01;
|
|
Packit |
67cb25 |
const double P2 = 4 * 3.7748947079307981766760e-08;
|
|
Packit |
67cb25 |
const double P3 = 4 * 2.6951514290790594840552e-15;
|
|
Packit |
67cb25 |
const double TwoPi = 2*(P1 + P2 + P3);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double y = GSL_SIGN(theta) * 2 * floor(fabs(theta)/TwoPi);
|
|
Packit |
67cb25 |
double r = ((theta - y*P1) - y*P2) - y*P3;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(r > M_PI) { r = (((r-2*P1)-2*P2)-2*P3); } /* r-TwoPi */
|
|
Packit |
67cb25 |
else if (r < -M_PI) r = (((r+2*P1)+2*P2)+2*P3); /* r+TwoPi */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = r;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(theta) > 0.0625/GSL_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->val = GSL_NAN;
|
|
Packit |
67cb25 |
result->err = GSL_NAN;
|
|
Packit |
67cb25 |
GSL_ERROR ("error", GSL_ELOSS);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(theta) > 0.0625/GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val - theta);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
double delta = fabs(result->val - theta);
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * ((delta < M_PI) ? delta : M_PI);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_angle_restrict_pos_err_e(const double theta, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* synthetic extended precision constants */
|
|
Packit |
67cb25 |
const double P1 = 4 * 7.85398125648498535156e-01;
|
|
Packit |
67cb25 |
const double P2 = 4 * 3.77489470793079817668e-08;
|
|
Packit |
67cb25 |
const double P3 = 4 * 2.69515142907905952645e-15;
|
|
Packit |
67cb25 |
const double TwoPi = 2*(P1 + P2 + P3);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double y = 2*floor(theta/TwoPi);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double r = ((theta - y*P1) - y*P2) - y*P3;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(r > TwoPi) {r = (((r-2*P1)-2*P2)-2*P3); } /* r-TwoPi */
|
|
Packit |
67cb25 |
else if (r < 0) { /* may happen due to FP rounding */
|
|
Packit |
67cb25 |
r = (((r+2*P1)+2*P2)+2*P3); /* r+TwoPi */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = r;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(theta) > 0.0625/GSL_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->val = GSL_NAN;
|
|
Packit |
67cb25 |
result->err = fabs(result->val);
|
|
Packit |
67cb25 |
GSL_ERROR ("error", GSL_ELOSS);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(theta) > 0.0625/GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(result->val - theta);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
double delta = fabs(result->val - theta);
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * ((delta < M_PI) ? delta : M_PI);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_angle_restrict_symm_e(double * theta)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_sf_result r;
|
|
Packit |
67cb25 |
int stat = gsl_sf_angle_restrict_symm_err_e(*theta, &r);
|
|
Packit |
67cb25 |
*theta = r.val;
|
|
Packit |
67cb25 |
return stat;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_angle_restrict_pos_e(double * theta)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_sf_result r;
|
|
Packit |
67cb25 |
int stat = gsl_sf_angle_restrict_pos_err_e(*theta, &r);
|
|
Packit |
67cb25 |
*theta = r.val;
|
|
Packit |
67cb25 |
return stat;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_sin_err_e(const double x, const double dx, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int stat_s = gsl_sf_sin_e(x, result);
|
|
Packit |
67cb25 |
result->err += fabs(cos(x) * dx);
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_cos_err_e(const double x, const double dx, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int stat_c = gsl_sf_cos_e(x, result);
|
|
Packit |
67cb25 |
result->err += fabs(sin(x) * dx);
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_c;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#if 0
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_sin_pi_x_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(-100.0 < x && x < 100.0) {
|
|
Packit |
67cb25 |
result->val = sin(M_PI * x) / (M_PI * x);
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double N = floor(x + 0.5);
|
|
Packit |
67cb25 |
const double f = x - N;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(N < INT_MAX && N > INT_MIN) {
|
|
Packit |
67cb25 |
/* Make it an integer if we can. Saves another
|
|
Packit |
67cb25 |
* call to floor().
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const int intN = (int)N;
|
|
Packit |
67cb25 |
const double sign = ( GSL_IS_ODD(intN) ? -1.0 : 1.0 );
|
|
Packit |
67cb25 |
result->val = sign * sin(M_PI * f);
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(N > 2.0/GSL_DBL_EPSILON || N < -2.0/GSL_DBL_EPSILON) {
|
|
Packit |
67cb25 |
/* All integer-valued floating point numbers
|
|
Packit |
67cb25 |
* bigger than 2/eps=2^53 are actually even.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double resN = N - 2.0*floor(0.5*N); /* 0 for even N, 1 for odd N */
|
|
Packit |
67cb25 |
const double sign = ( fabs(resN) > 0.5 ? -1.0 : 1.0 );
|
|
Packit |
67cb25 |
result->val = sign * sin(M_PI*f);
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_sinc_e(double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double ax = fabs(x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(ax < 0.8) {
|
|
Packit |
67cb25 |
/* Do not go to the limit of the fit since
|
|
Packit |
67cb25 |
* there is a zero there and the Chebyshev
|
|
Packit |
67cb25 |
* accuracy will go to zero.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return cheb_eval_e(&sinc_cs, 2.0*ax-1.0, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(ax < 100.0) {
|
|
Packit |
67cb25 |
/* Small arguments are no problem.
|
|
Packit |
67cb25 |
* We trust the library sin() to
|
|
Packit |
67cb25 |
* roughly machine precision.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
result->val = sin(M_PI * ax)/(M_PI * ax);
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* Large arguments must be handled separately.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const double r = M_PI*ax;
|
|
Packit |
67cb25 |
gsl_sf_result s;
|
|
Packit |
67cb25 |
int stat_s = gsl_sf_sin_e(r, &s);
|
|
Packit |
67cb25 |
result->val = s.val/r;
|
|
Packit |
67cb25 |
result->err = s.err/r + 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
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_sin(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_sin_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_cos(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_cos_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_hypot(const double x, const double y)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_hypot_e(x, y, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_lnsinh(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_lnsinh_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_lncosh(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_lncosh_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_angle_restrict_symm(const double theta)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double result = theta;
|
|
Packit |
67cb25 |
EVAL_DOUBLE(gsl_sf_angle_restrict_symm_e(&result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_angle_restrict_pos(const double theta)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double result = theta;
|
|
Packit |
67cb25 |
EVAL_DOUBLE(gsl_sf_angle_restrict_pos_e(&result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#if 0
|
|
Packit |
67cb25 |
double gsl_sf_sin_pi_x(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_sin_pi_x_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_sinc(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_sinc_e(x, &result));
|
|
Packit |
67cb25 |
}
|