htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/exp.c

   1 /* specfunc/exp.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_gamma.h>
  26 #include <gsl/gsl_sf_exp.h>
  27
  28 #include "error.h"
  29
  30 /* Evaluate the continued fraction for exprel.
  31  * [Abramowitz+Stegun, 4.2.41]
  32  */
  33 static
  34 int
  35 exprel_n_CF(const int N, const double x, gsl_sf_result * result)
  36 {
  37   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  38   const int maxiter = 5000;
  39   int n = 1;
  40   double Anm2 = 1.0;
  41   double Bnm2 = 0.0;
  42   double Anm1 = 0.0;
  43   double Bnm1 = 1.0;
  44   double a1 = 1.0;
  45   double b1 = 1.0;
  46   double a2 = -x;
  47   double b2 = N+1;
  48   double an, bn;
  49
  50   double fn;
  51
  52   double An = b1*Anm1 + a1*Anm2;   /* A1 */
  53   double Bn = b1*Bnm1 + a1*Bnm2;   /* B1 */
  54
  55   /* One explicit step, before we get to the main pattern. */
  56   n++;
  57   Anm2 = Anm1;
  58   Bnm2 = Bnm1;
  59   Anm1 = An;
  60   Bnm1 = Bn;
  61   An = b2*Anm1 + a2*Anm2;   /* A2 */
  62   Bn = b2*Bnm1 + a2*Bnm2;   /* B2 */
  63
  64   fn = An/Bn;
  65
  66   while(n < maxiter) {
  67     double old_fn;
  68     double del;
  69     n++;
  70     Anm2 = Anm1;
  71     Bnm2 = Bnm1;
  72     Anm1 = An;
  73     Bnm1 = Bn;
  74     an = ( GSL_IS_ODD(n) ? ((n-1)/2)*x : -(N+(n/2)-1)*x );
  75     bn = N + n - 1;
  76     An = bn*Anm1 + an*Anm2;
  77     Bn = bn*Bnm1 + an*Bnm2;
  78
  79     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
  80       An /= RECUR_BIG;
  81       Bn /= RECUR_BIG;
  82       Anm1 /= RECUR_BIG;
  83       Bnm1 /= RECUR_BIG;
  84       Anm2 /= RECUR_BIG;
  85       Bnm2 /= RECUR_BIG;
  86     }
  87
  88     old_fn = fn;
  89     fn = An/Bn;
  90     del = old_fn/fn;
  91
  92     if(fabs(del - 1.0) < 2.0*GSL_DBL_EPSILON) break;
  93   }
  94
  95   result->val = fn;
  96   result->err = 2.0*(n+1.0)*GSL_DBL_EPSILON*fabs(fn);
  97
  98   if(n == maxiter)
  99     GSL_ERROR ("error", GSL_EMAXITER);
 100   else
 101     return GSL_SUCCESS;
 102 }
 103
 104
 105 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 106
 107 int gsl_sf_exp_e(const double x, gsl_sf_result * result)
 108 {
 109   if(x > GSL_LOG_DBL_MAX) {
 110     OVERFLOW_ERROR(result);
 111   }
 112   else if(x < GSL_LOG_DBL_MIN) {
 113     UNDERFLOW_ERROR(result);
 114   }
 115   else {
 116     result->val = exp(x);
 117     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 118     return GSL_SUCCESS;
 119   }
 120 }
 121
 122 int gsl_sf_exp_e10_e(const double x, gsl_sf_result_e10 * result)
 123 {
 124   if(x > INT_MAX-1) {
 125     OVERFLOW_ERROR_E10(result);
 126   }
 127   else if(x < INT_MIN+1) {
 128     UNDERFLOW_ERROR_E10(result);
 129   }
 130   else {
 131     const int N = (x > GSL_LOG_DBL_MAX || x < GSL_LOG_DBL_MIN) ? (int) floor(x/M_LN10) : 0;
 132     result->val = exp(x-N*M_LN10);
 133     result->err = 2.0 * (fabs(x)+1.0) * GSL_DBL_EPSILON * fabs(result->val);
 134     result->e10 = N;
 135     return GSL_SUCCESS;
 136   }
 137 }
 138
 139
 140 int gsl_sf_exp_mult_e(const double x, const double y, gsl_sf_result * result)
 141 {
 142   const double ay  = fabs(y);
 143
 144   if(y == 0.0) {
 145     result->val = 0.0;
 146     result->err = 0.0;
 147     return GSL_SUCCESS;
 148   }
 149   else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
 150           && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
 151     ) {
 152     const double ex = exp(x);
 153     result->val = y * ex;
 154     result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
 155     return GSL_SUCCESS;
 156   }
 157   else {
 158     const double ly  = log(ay);
 159     const double lnr = x + ly;
 160
 161     if(lnr > GSL_LOG_DBL_MAX - 0.01) {
 162       OVERFLOW_ERROR(result);
 163     }
 164     else if(lnr < GSL_LOG_DBL_MIN + 0.01) {
 165       UNDERFLOW_ERROR(result);
 166     }
 167     else {
 168       const double sy   = GSL_SIGN(y);
 169       const double M    = floor(x);
 170       const double N    = floor(ly);
 171       const double a    = x  - M;
 172       const double b    = ly - N;
 173       const double berr = 2.0 * GSL_DBL_EPSILON * (fabs(ly) + fabs(N));
 174       result->val  = sy * exp(M+N) * exp(a+b);
 175       result->err  = berr * fabs(result->val);
 176       result->err += 2.0 * GSL_DBL_EPSILON * (M + N + 1.0) * fabs(result->val);
 177       return GSL_SUCCESS;
 178     }
 179   }
 180 }
 181
 182
 183 int gsl_sf_exp_mult_e10_e(const double x, const double y, gsl_sf_result_e10 * result)
 184 {
 185   const double ay  = fabs(y);
 186
 187   if(y == 0.0) {
 188     result->val = 0.0;
 189     result->err = 0.0;
 190     result->e10 = 0;
 191     return GSL_SUCCESS;
 192   }
 193   else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
 194           && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
 195     ) {
 196     const double ex = exp(x);
 197     result->val = y * ex;
 198     result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
 199     result->e10 = 0;
 200     return GSL_SUCCESS;
 201   }
 202   else {
 203     const double ly  = log(ay);
 204     const double l10_val = (x + ly)/M_LN10;
 205
 206     if(l10_val > INT_MAX-1) {
 207       OVERFLOW_ERROR_E10(result);
 208     }
 209     else if(l10_val < INT_MIN+1) {
 210       UNDERFLOW_ERROR_E10(result);
 211     }
 212     else {
 213       const double sy  = GSL_SIGN(y);
 214       const int    N   = (int) floor(l10_val);
 215       const double arg_val = (l10_val - N) * M_LN10;
 216       const double arg_err = 2.0 * GSL_DBL_EPSILON * (fabs(x) + fabs(ly) + M_LN10*fabs(N));
 217
 218       result->val  = sy * exp(arg_val);
 219       result->err  = arg_err * fabs(result->val);
 220       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 221       result->e10 = N;
 222
 223       return GSL_SUCCESS;
 224     }
 225   }
 226 }
 227
 228
 229 int gsl_sf_exp_mult_err_e(const double x, const double dx,
 230                              const double y, const double dy,
 231                              gsl_sf_result * result)
 232 {
 233   const double ay  = fabs(y);
 234
 235   if(y == 0.0) {
 236     result->val = 0.0;
 237     result->err = fabs(dy * exp(x));
 238     return GSL_SUCCESS;
 239   }
 240   else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
 241           && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
 242     ) {
 243     double ex = exp(x);
 244     result->val  = y * ex;
 245     result->err  = ex * (fabs(dy) + fabs(y*dx));
 246     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 247     return GSL_SUCCESS;
 248   }
 249   else {
 250     const double ly  = log(ay);
 251     const double lnr = x + ly;
 252
 253     if(lnr > GSL_LOG_DBL_MAX - 0.01) {
 254       OVERFLOW_ERROR(result);
 255     }
 256     else if(lnr < GSL_LOG_DBL_MIN + 0.01) {
 257       UNDERFLOW_ERROR(result);
 258     }
 259     else {
 260       const double sy  = GSL_SIGN(y);
 261       const double M   = floor(x);
 262       const double N   = floor(ly);
 263       const double a   = x  - M;
 264       const double b   = ly - N;
 265       const double eMN = exp(M+N);
 266       const double eab = exp(a+b);
 267       result->val  = sy * eMN * eab;
 268       result->err  = eMN * eab * 2.0*GSL_DBL_EPSILON;
 269       result->err += eMN * eab * fabs(dy/y);
 270       result->err += eMN * eab * fabs(dx);
 271       return GSL_SUCCESS;
 272     }
 273   }
 274 }
 275
 276
 277 int gsl_sf_exp_mult_err_e10_e(const double x, const double dx,
 278                              const double y, const double dy,
 279                              gsl_sf_result_e10 * result)
 280 {
 281   const double ay  = fabs(y);
 282
 283   if(y == 0.0) {
 284     result->val = 0.0;
 285     result->err = fabs(dy * exp(x));
 286     result->e10 = 0;
 287     return GSL_SUCCESS;
 288   }
 289   else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
 290           && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
 291     ) {
 292     const double ex = exp(x);
 293     result->val  = y * ex;
 294     result->err  = ex * (fabs(dy) + fabs(y*dx));
 295     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 296     result->e10 = 0;
 297     return GSL_SUCCESS;
 298   }
 299   else {
 300     const double ly  = log(ay);
 301     const double l10_val = (x + ly)/M_LN10;
 302
 303     if(l10_val > INT_MAX-1) {
 304       OVERFLOW_ERROR_E10(result);
 305     }
 306     else if(l10_val < INT_MIN+1) {
 307       UNDERFLOW_ERROR_E10(result);
 308     }
 309     else {
 310       const double sy  = GSL_SIGN(y);
 311       const int    N   = (int) floor(l10_val);
 312       const double arg_val = (l10_val - N) * M_LN10;
 313       const double arg_err = dy/fabs(y) + dx + 2.0*GSL_DBL_EPSILON*fabs(arg_val);
 314
 315       result->val  = sy * exp(arg_val);
 316       result->err  = arg_err * fabs(result->val);
 317       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 318       result->e10 = N;
 319
 320       return GSL_SUCCESS;
 321     }
 322   }
 323 }
 324
 325
 326 int gsl_sf_expm1_e(const double x, gsl_sf_result * result)
 327 {
 328   const double cut = 0.002;
 329
 330   if(x < GSL_LOG_DBL_MIN) {
 331     result->val = -1.0;
 332     result->err = GSL_DBL_EPSILON;
 333     return GSL_SUCCESS;
 334   }
 335   else if(x < -cut) {
 336     result->val = exp(x) - 1.0;
 337     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 338     return GSL_SUCCESS;
 339   }
 340   else if(x < cut) {
 341     result->val = x * (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x))));
 342     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 343     return GSL_SUCCESS;
 344   }
 345   else if(x < GSL_LOG_DBL_MAX) {
 346     result->val = exp(x) - 1.0;
 347     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 348     return GSL_SUCCESS;
 349   }
 350   else {
 351     OVERFLOW_ERROR(result);
 352   }
 353 }
 354
 355
 356 int gsl_sf_exprel_e(const double x, gsl_sf_result * result)
 357 {
 358   const double cut = 0.002;
 359
 360   if(x < GSL_LOG_DBL_MIN) {
 361     result->val = -1.0/x;
 362     result->err = GSL_DBL_EPSILON * fabs(result->val);
 363     return GSL_SUCCESS;
 364   }
 365   else if(x < -cut) {
 366     result->val = (exp(x) - 1.0)/x;
 367     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 368     return GSL_SUCCESS;
 369   }
 370   else if(x < cut) {
 371     result->val = (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x))));
 372     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 373     return GSL_SUCCESS;
 374   }
 375   else if(x < GSL_LOG_DBL_MAX) {
 376     result->val = (exp(x) - 1.0)/x;
 377     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 378     return GSL_SUCCESS;
 379   }
 380   else {
 381     OVERFLOW_ERROR(result);
 382   }
 383 }
 384
 385
 386 int gsl_sf_exprel_2_e(double x, gsl_sf_result * result)
 387 {
 388   const double cut = 0.002;
 389
 390   if(x < GSL_LOG_DBL_MIN) {
 391     result->val = -2.0/x*(1.0 + 1.0/x);
 392     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 393     return GSL_SUCCESS;
 394   }
 395   else if(x < -cut) {
 396     result->val = 2.0*(exp(x) - 1.0 - x)/(x*x);
 397     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 398     return GSL_SUCCESS;
 399   }
 400   else if(x < cut) {
 401     result->val = (1.0 + 1.0/3.0*x*(1.0 + 0.25*x*(1.0 + 0.2*x*(1.0 + 1.0/6.0*x))));
 402     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 403     return GSL_SUCCESS;
 404   }
 405   else if(x < GSL_LOG_DBL_MAX) {
 406     result->val = 2.0*(exp(x) - 1.0 - x)/(x*x);
 407     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 408     return GSL_SUCCESS;
 409   }
 410   else {
 411     OVERFLOW_ERROR(result);
 412   }
 413 }
 414
 415
 416 int
 417 gsl_sf_exprel_n_e(const int N, const double x, gsl_sf_result * result)
 418 {
 419   if(N < 0) {
 420     DOMAIN_ERROR(result);
 421   }
 422   else if(x == 0.0) {
 423     result->val = 1.0;
 424     result->err = 0.0;
 425     return GSL_SUCCESS;
 426   }
 427   else if(fabs(x) < GSL_ROOT3_DBL_EPSILON * N) {
 428     result->val = 1.0 + x/(N+1) * (1.0 + x/(N+2));
 429     result->err = 2.0 * GSL_DBL_EPSILON;
 430     return GSL_SUCCESS;
 431   }
 432   else if(N == 0) {
 433     return gsl_sf_exp_e(x, result);
 434   }
 435   else if(N == 1) {
 436     return gsl_sf_exprel_e(x, result);
 437   }
 438   else if(N == 2) {
 439     return gsl_sf_exprel_2_e(x, result);
 440   }
 441   else {
 442     if(x > N && (-x + N*(1.0 + log(x/N)) < GSL_LOG_DBL_EPSILON)) {
 443       /* x is much larger than n.
 444        * Ignore polynomial part, so
 445        * exprel_N(x) ~= e^x N!/x^N
 446        */
 447       gsl_sf_result lnf_N;
 448       double lnr_val;
 449       double lnr_err;
 450       double lnterm;
 451       gsl_sf_lnfact_e(N, &lnf_N);
 452       lnterm = N*log(x);
 453       lnr_val  = x + lnf_N.val - lnterm;
 454       lnr_err  = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(lnterm));
 455       lnr_err += lnf_N.err;
 456       return gsl_sf_exp_err_e(lnr_val, lnr_err, result);
 457     }
 458     else if(x > N) {
 459       /* Write the identity
 460        *   exprel_n(x) = e^x n! / x^n (1 - Gamma[n,x]/Gamma[n])
 461        * then use the asymptotic expansion
 462        * Gamma[n,x] ~ x^(n-1) e^(-x) (1 + (n-1)/x + (n-1)(n-2)/x^2 + ...)
 463        */
 464       double ln_x = log(x);
 465       gsl_sf_result lnf_N;
 466       double lg_N;
 467       double lnpre_val;
 468       double lnpre_err;
 469       gsl_sf_lnfact_e(N, &lnf_N);    /* log(N!)       */
 470       lg_N  = lnf_N.val - log(N);       /* log(Gamma(N)) */
 471       lnpre_val  = x + lnf_N.val - N*ln_x;
 472       lnpre_err  = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(N*ln_x));
 473       lnpre_err += lnf_N.err;
 474       if(lnpre_val < GSL_LOG_DBL_MAX - 5.0) {
 475         int stat_eG;
 476         gsl_sf_result bigG_ratio;
 477         gsl_sf_result pre;
 478         int stat_ex = gsl_sf_exp_err_e(lnpre_val, lnpre_err, &pre);
 479         double ln_bigG_ratio_pre = -x + (N-1)*ln_x - lg_N;
 480         double bigGsum = 1.0;
 481         double term = 1.0;
 482         int k;
 483         for(k=1; k<N; k++) {
 484           term *= (N-k)/x;
 485           bigGsum += term;
 486         }
 487         stat_eG = gsl_sf_exp_mult_e(ln_bigG_ratio_pre, bigGsum, &bigG_ratio);
 488         if(stat_eG == GSL_SUCCESS) {
 489           result->val  = pre.val * (1.0 - bigG_ratio.val);
 490           result->err  = pre.val * (2.0*GSL_DBL_EPSILON + bigG_ratio.err);
 491           result->err += pre.err * fabs(1.0 - bigG_ratio.val);
 492           result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 493           return stat_ex;
 494         }
 495         else {
 496           result->val = 0.0;
 497           result->err = 0.0;
 498           return stat_eG;
 499         }
 500       }
 501       else {
 502         OVERFLOW_ERROR(result);
 503       }
 504     }
 505     else if(x > -10.0*N) {
 506       return exprel_n_CF(N, x, result);
 507     }
 508     else {
 509       /* x -> -Inf asymptotic:
 510        * exprel_n(x) ~ e^x n!/x^n - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...)
 511        *             ~ - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...)
 512        */
 513       double sum  = 1.0;
 514       double term = 1.0;
 515       int k;
 516       for(k=1; k<N; k++) {
 517         term *= (N-k)/x;
 518         sum  += term;
 519       }
 520       result->val = -N/x * sum;
 521       result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 522       return GSL_SUCCESS;
 523     }
 524   }
 525 }
 526
 527
 528 int
 529 gsl_sf_exp_err_e(const double x, const double dx, gsl_sf_result * result)
 530 {
 531   const double adx = fabs(dx);
 532
 533   /* CHECK_POINTER(result) */
 534
 535   if(x + adx > GSL_LOG_DBL_MAX) {
 536     OVERFLOW_ERROR(result);
 537   }
 538   else if(x - adx < GSL_LOG_DBL_MIN) {
 539     UNDERFLOW_ERROR(result);
 540   }
 541   else {
 542     const double ex  = exp(x);
 543     const double edx = exp(adx);
 544     result->val  = ex;
 545     result->err  = ex * GSL_MAX_DBL(GSL_DBL_EPSILON, edx - 1.0/edx);
 546     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 547     return GSL_SUCCESS;
 548   }
 549 }
 550
 551
 552 int
 553 gsl_sf_exp_err_e10_e(const double x, const double dx, gsl_sf_result_e10 * result)
 554 {
 555   const double adx = fabs(dx);
 556
 557   /* CHECK_POINTER(result) */
 558
 559   if(x + adx > INT_MAX - 1) {
 560     OVERFLOW_ERROR_E10(result);
 561   }
 562   else if(x - adx < INT_MIN + 1) {
 563     UNDERFLOW_ERROR_E10(result);
 564   }
 565   else {
 566     const int    N  = (int)floor(x/M_LN10);
 567     const double ex = exp(x-N*M_LN10);
 568     result->val = ex;
 569     result->err = ex * (2.0 * GSL_DBL_EPSILON * (fabs(x) + 1.0) + adx);
 570     result->e10 = N;
 571     return GSL_SUCCESS;
 572   }
 573 }
 574
 575
 576 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 577
 578 #include "eval.h"
 579
 580 double gsl_sf_exp(const double x)
 581 {
 582   EVAL_RESULT(gsl_sf_exp_e(x, &result));
 583 }
 584
 585 double gsl_sf_exp_mult(const double x, const double y)
 586 {
 587   EVAL_RESULT(gsl_sf_exp_mult_e(x, y, &result));
 588 }
 589
 590 double gsl_sf_expm1(const double x)
 591 {
 592   EVAL_RESULT(gsl_sf_expm1_e(x, &result));
 593 }
 594
 595 double gsl_sf_exprel(const double x)
 596 {
 597   EVAL_RESULT(gsl_sf_exprel_e(x, &result));
 598 }
 599
 600 double gsl_sf_exprel_2(const double x)
 601 {
 602   EVAL_RESULT(gsl_sf_exprel_2_e(x, &result));
 603 }
 604
 605 double gsl_sf_exprel_n(const int n, const double x)
 606 {
 607   EVAL_RESULT(gsl_sf_exprel_n_e(n, x, &result));
 608 }