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