Mercurial > octave-nkf

comparison liboctave/lo-specfun.cc @ 11586:12df7854fa7c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
strip trailing whitespace from source files
author John W. Eaton <jwe@octave.org>
date Thu, 20 Jan 2011 17:24:59 -0500
parents a83bad07f7e3
children 72c96de7a403
comparison
equal deleted inserted replaced
11585:1473d0cf86d2 11586:12df7854fa7c
331 { 331 {
332 double result; 332 double result;
333 333
334 #if defined (HAVE_LGAMMA_R) 334 #if defined (HAVE_LGAMMA_R)
335 int sgngam; 335 int sgngam;
336 result = lgamma_r (x, &sgngam); 336 result = lgamma_r (x, &sgngam);
337 #else 337 #else
338 double sgngam; 338 double sgngam;
339 339
340 if (xisnan (x)) 340 if (xisnan (x))
341 result = x; 341 result = x;
396 { 396 {
397 float result; 397 float result;
398 398
399 #if defined (HAVE_LGAMMAF_R) 399 #if defined (HAVE_LGAMMAF_R)
400 int sgngam; 400 int sgngam;
401 result = lgammaf_r (x, &sgngam); 401 result = lgammaf_r (x, &sgngam);
402 #else 402 #else
403 float sgngam; 403 float sgngam;
404 404
405 if (xisnan (x)) 405 if (xisnan (x))
406 result = x; 406 result = x;
429 { 429 {
430 ax /= 16; 430 ax /= 16;
431 431
432 // use Taylor series to calculate exp(x)-1. 432 // use Taylor series to calculate exp(x)-1.
433 double t = ax; 433 double t = ax;
434 double s = 0; 434 double s = 0;
435 for (int i = 2; i < 7; i++) 435 for (int i = 2; i < 7; i++)
436 s += (t *= ax/i); 436 s += (t *= ax/i);
437 s += ax; 437 s += ax;
438 438
439 // use the identity (a+1)^2-1 = a*(a+2) 439 // use the identity (a+1)^2-1 = a*(a+2)
451 451
452 return retval; 452 return retval;
453 } 453 }
454 #endif 454 #endif
455 455
456 Complex 456 Complex
457 expm1(const Complex& x) 457 expm1(const Complex& x)
458 { 458 {
459 Complex retval; 459 Complex retval;
460 460
461 if (std:: abs (x) < 1) 461 if (std:: abs (x) < 1)
484 { 484 {
485 ax /= 16; 485 ax /= 16;
486 486
487 // use Taylor series to calculate exp(x)-1. 487 // use Taylor series to calculate exp(x)-1.
488 float t = ax; 488 float t = ax;
489 float s = 0; 489 float s = 0;
490 for (int i = 2; i < 7; i++) 490 for (int i = 2; i < 7; i++)
491 s += (t *= ax/i); 491 s += (t *= ax/i);
492 s += ax; 492 s += ax;
493 493
494 // use the identity (a+1)^2-1 = a*(a+2) 494 // use the identity (a+1)^2-1 = a*(a+2)
506 506
507 return retval; 507 return retval;
508 } 508 }
509 #endif 509 #endif
510 510
511 FloatComplex 511 FloatComplex
512 expm1(const FloatComplex& x) 512 expm1(const FloatComplex& x)
513 { 513 {
514 FloatComplex retval; 514 FloatComplex retval;
515 515
516 if (std:: abs (x) < 1) 516 if (std:: abs (x) < 1)
549 549
550 return retval; 550 return retval;
551 } 551 }
552 #endif 552 #endif
553 553
554 Complex 554 Complex
555 log1p (const Complex& x) 555 log1p (const Complex& x)
556 { 556 {
557 Complex retval; 557 Complex retval;
558 558
559 double r = x.real (), i = x.imag(); 559 double r = x.real (), i = x.imag();
608 608
609 return retval; 609 return retval;
610 } 610 }
611 #endif 611 #endif
612 612
613 FloatComplex 613 FloatComplex
614 log1p (const FloatComplex& x) 614 log1p (const FloatComplex& x)
615 { 615 {
616 FloatComplex retval; 616 FloatComplex retval;
617 617
618 float r = x.real (), i = x.imag(); 618 float r = x.real (), i = x.imag();
713 713
714 F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr); 714 F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr);
715 715
716 if (kode != 2) 716 if (kode != 2)
717 { 717 {
718 double expz = exp (std::abs (zi)); 718 double expz = exp (std::abs (zi));
719 yr *= expz; 719 yr *= expz;
720 yi *= expz; 720 yi *= expz;
721 } 721 }
722 722
723 if (zi == 0.0 && zr >= 0.0) 723 if (zi == 0.0 && zr >= 0.0)
728 else if (is_integer_value (alpha)) 728 else if (is_integer_value (alpha))
729 { 729 {
730 // zbesy can overflow as z->0, and cause troubles for generic case below 730 // zbesy can overflow as z->0, and cause troubles for generic case below
731 alpha = -alpha; 731 alpha = -alpha;
732 Complex tmp = zbesj (z, alpha, kode, ierr); 732 Complex tmp = zbesj (z, alpha, kode, ierr);
733 if ((static_cast <long> (alpha)) & 1) 733 if ((static_cast <long> (alpha)) & 1)
734 tmp = - tmp; 734 tmp = - tmp;
735 retval = bessel_return_value (tmp, ierr); 735 retval = bessel_return_value (tmp, ierr);
736 } 736 }
737 else 737 else
738 { 738 {
798 else if (is_integer_value (alpha - 0.5)) 798 else if (is_integer_value (alpha - 0.5))
799 { 799 {
800 // zbesy can overflow as z->0, and cause troubles for generic case below 800 // zbesy can overflow as z->0, and cause troubles for generic case below
801 alpha = -alpha; 801 alpha = -alpha;
802 Complex tmp = zbesj (z, alpha, kode, ierr); 802 Complex tmp = zbesj (z, alpha, kode, ierr);
803 if ((static_cast <long> (alpha - 0.5)) & 1) 803 if ((static_cast <long> (alpha - 0.5)) & 1)
804 tmp = - tmp; 804 tmp = - tmp;
805 retval = bessel_return_value (tmp, ierr); 805 retval = bessel_return_value (tmp, ierr);
806 } 806 }
807 else 807 else
808 { 808 {
860 860
861 if (ierr == 0 || ierr == 3) 861 if (ierr == 0 || ierr == 3)
862 { 862 {
863 Complex tmp2 = (2.0 / M_PI) * sin (M_PI * alpha) 863 Complex tmp2 = (2.0 / M_PI) * sin (M_PI * alpha)
864 * zbesk (z, alpha, kode, ierr); 864 * zbesk (z, alpha, kode, ierr);
865 865
866 if (kode == 2) 866 if (kode == 2)
867 { 867 {
868 // Compensate for different scaling factor of besk. 868 // Compensate for different scaling factor of besk.
869 tmp2 *= exp(-z - std::abs(z.real())); 869 tmp2 *= exp(-z - std::abs(z.real()));
870 } 870 }
871 871
872 tmp += tmp2; 872 tmp += tmp2;
873 873
874 retval = bessel_return_value (tmp, ierr); 874 retval = bessel_return_value (tmp, ierr);
875 } 875 }
876 else 876 else
1350 else if (is_integer_value (alpha)) 1350 else if (is_integer_value (alpha))
1351 { 1351 {
1352 // zbesy can overflow as z->0, and cause troubles for generic case below 1352 // zbesy can overflow as z->0, and cause troubles for generic case below
1353 alpha = -alpha; 1353 alpha = -alpha;
1354 FloatComplex tmp = cbesj (z, alpha, kode, ierr); 1354 FloatComplex tmp = cbesj (z, alpha, kode, ierr);
1355 if ((static_cast <long> (alpha)) & 1) 1355 if ((static_cast <long> (alpha)) & 1)
1356 tmp = - tmp; 1356 tmp = - tmp;
1357 retval = bessel_return_value (tmp, ierr); 1357 retval = bessel_return_value (tmp, ierr);
1358 } 1358 }
1359 else 1359 else
1360 { 1360 {
1413 else if (is_integer_value (alpha - 0.5)) 1413 else if (is_integer_value (alpha - 0.5))
1414 { 1414 {
1415 // zbesy can overflow as z->0, and cause troubles for generic case below 1415 // zbesy can overflow as z->0, and cause troubles for generic case below
1416 alpha = -alpha; 1416 alpha = -alpha;
1417 FloatComplex tmp = cbesj (z, alpha, kode, ierr); 1417 FloatComplex tmp = cbesj (z, alpha, kode, ierr);
1418 if ((static_cast <long> (alpha - 0.5)) & 1) 1418 if ((static_cast <long> (alpha - 0.5)) & 1)
1419 tmp = - tmp; 1419 tmp = - tmp;
1420 retval = bessel_return_value (tmp, ierr); 1420 retval = bessel_return_value (tmp, ierr);
1421 } 1421 }
1422 else 1422 else
1423 { 1423 {
1470 1470
1471 if (ierr == 0 || ierr == 3) 1471 if (ierr == 0 || ierr == 3)
1472 { 1472 {
1473 FloatComplex tmp2 = static_cast<float> (2.0 / M_PI) * sinf (static_cast<float> (M_PI) * alpha) 1473 FloatComplex tmp2 = static_cast<float> (2.0 / M_PI) * sinf (static_cast<float> (M_PI) * alpha)
1474 * cbesk (z, alpha, kode, ierr); 1474 * cbesk (z, alpha, kode, ierr);
1475 1475
1476 if (kode == 2) 1476 if (kode == 2)
1477 { 1477 {
1478 // Compensate for different scaling factor of besk. 1478 // Compensate for different scaling factor of besk.
1479 tmp2 *= exp(-z - std::abs(z.real())); 1479 tmp2 *= exp(-z - std::abs(z.real()));
1480 } 1480 }
1481 1481
2257 for (octave_idx_type i = 0; i < nel; i++) 2257 for (octave_idx_type i = 0; i < nel; i++)
2258 retval (i) = betainc (x, a(i), b(i)); 2258 retval (i) = betainc (x, a(i), b(i));
2259 } 2259 }
2260 else 2260 else
2261 gripe_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ()); 2261 gripe_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ());
2262 2262
2263 return retval; 2263 return retval;
2264 } 2264 }
2265 2265
2266 2266
2267 Matrix 2267 Matrix
2386 for (octave_idx_type i = 0; i < nel; i++) 2386 for (octave_idx_type i = 0; i < nel; i++)
2387 retval (i) = betainc (x(i), a, b(i)); 2387 retval (i) = betainc (x(i), a, b(i));
2388 } 2388 }
2389 else 2389 else
2390 gripe_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ()); 2390 gripe_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ());
2391 2391
2392 return retval; 2392 return retval;
2393 } 2393 }
2394 2394
2395 NDArray 2395 NDArray
2396 betainc (const NDArray& x, const NDArray& a, double b) 2396 betainc (const NDArray& x, const NDArray& a, double b)
2407 for (octave_idx_type i = 0; i < nel; i++) 2407 for (octave_idx_type i = 0; i < nel; i++)
2408 retval (i) = betainc (x(i), a(i), b); 2408 retval (i) = betainc (x(i), a(i), b);
2409 } 2409 }
2410 else 2410 else
2411 gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0)); 2411 gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0));
2412 2412
2413 return retval; 2413 return retval;
2414 } 2414 }
2415 2415
2416 NDArray 2416 NDArray
2417 betainc (const NDArray& x, const NDArray& a, const NDArray& b) 2417 betainc (const NDArray& x, const NDArray& a, const NDArray& b)
2540 for (octave_idx_type i = 0; i < nel; i++) 2540 for (octave_idx_type i = 0; i < nel; i++)
2541 retval (i) = betainc (x, a(i), b(i)); 2541 retval (i) = betainc (x, a(i), b(i));
2542 } 2542 }
2543 else 2543 else
2544 gripe_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ()); 2544 gripe_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ());
2545 2545
2546 return retval; 2546 return retval;
2547 } 2547 }
2548 2548
2549 2549
2550 FloatMatrix 2550 FloatMatrix
2669 for (octave_idx_type i = 0; i < nel; i++) 2669 for (octave_idx_type i = 0; i < nel; i++)
2670 retval (i) = betainc (x(i), a, b(i)); 2670 retval (i) = betainc (x(i), a, b(i));
2671 } 2671 }
2672 else 2672 else
2673 gripe_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ()); 2673 gripe_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ());
2674 2674
2675 return retval; 2675 return retval;
2676 } 2676 }
2677 2677
2678 FloatNDArray 2678 FloatNDArray
2679 betainc (const FloatNDArray& x, const FloatNDArray& a, float b) 2679 betainc (const FloatNDArray& x, const FloatNDArray& a, float b)
2690 for (octave_idx_type i = 0; i < nel; i++) 2690 for (octave_idx_type i = 0; i < nel; i++)
2691 retval (i) = betainc (x(i), a(i), b); 2691 retval (i) = betainc (x(i), a(i), b);
2692 } 2692 }
2693 else 2693 else
2694 gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0)); 2694 gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0));
2695 2695
2696 return retval; 2696 return retval;
2697 } 2697 }
2698 2698
2699 FloatNDArray 2699 FloatNDArray
2700 betainc (const FloatNDArray& x, const FloatNDArray& a, const FloatNDArray& b) 2700 betainc (const FloatNDArray& x, const FloatNDArray& a, const FloatNDArray& b)
2900 bool err; 2900 bool err;
2901 2901
2902 for (octave_idx_type i = 0; i < nel; i++) 2902 for (octave_idx_type i = 0; i < nel; i++)
2903 { 2903 {
2904 result (i) = gammainc (x(i), a(i), err); 2904 result (i) = gammainc (x(i), a(i), err);
2905 2905
2906 if (err) 2906 if (err)
2907 goto done; 2907 goto done;
2908 } 2908 }
2909 2909
2910 retval = result; 2910 retval = result;
3105 bool err; 3105 bool err;
3106 3106
3107 for (octave_idx_type i = 0; i < nel; i++) 3107 for (octave_idx_type i = 0; i < nel; i++)
3108 { 3108 {
3109 result (i) = gammainc (x(i), a(i), err); 3109 result (i) = gammainc (x(i), a(i), err);
3110 3110
3111 if (err) 3111 if (err)
3112 goto done; 3112 goto done;
3113 } 3113 }
3114 3114
3115 retval = result; 3115 retval = result;
3150 // faster evaluation. 3150 // faster evaluation.
3151 3151
3152 static double do_erfinv (double x, bool refine) 3152 static double do_erfinv (double x, bool refine)
3153 { 3153 {
3154 // Coefficients of rational approximation. 3154 // Coefficients of rational approximation.
3155 static const double a[] = 3155 static const double a[] =
3156 { -2.806989788730439e+01, 1.562324844726888e+02, 3156 { -2.806989788730439e+01, 1.562324844726888e+02,
3157 -1.951109208597547e+02, 9.783370457507161e+01, 3157 -1.951109208597547e+02, 9.783370457507161e+01,
3158 -2.168328665628878e+01, 1.772453852905383e+00 }; 3158 -2.168328665628878e+01, 1.772453852905383e+00 };
3159 static const double b[] = 3159 static const double b[] =
3160 { -5.447609879822406e+01, 1.615858368580409e+02, 3160 { -5.447609879822406e+01, 1.615858368580409e+02,
3161 -1.556989798598866e+02, 6.680131188771972e+01, 3161 -1.556989798598866e+02, 6.680131188771972e+01,
3162 -1.328068155288572e+01 }; 3162 -1.328068155288572e+01 };
3163 static const double c[] = 3163 static const double c[] =
3164 { -5.504751339936943e-03, -2.279687217114118e-01, 3164 { -5.504751339936943e-03, -2.279687217114118e-01,
3165 -1.697592457770869e+00, -1.802933168781950e+00, 3165 -1.697592457770869e+00, -1.802933168781950e+00,
3166 3.093354679843505e+00, 2.077595676404383e+00 }; 3166 3.093354679843505e+00, 2.077595676404383e+00 };
3167 static const double d[] = 3167 static const double d[] =
3168 { 7.784695709041462e-03, 3.224671290700398e-01, 3168 { 7.784695709041462e-03, 3.224671290700398e-01,
3169 2.445134137142996e+00, 3.754408661907416e+00 }; 3169 2.445134137142996e+00, 3.754408661907416e+00 };
3170 3170
3171 static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2. 3171 static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2.
3172 static const double pbreak = 0.95150; 3172 static const double pbreak = 0.95150;
3209 return do_erfinv (x, true); 3209 return do_erfinv (x, true);
3210 } 3210 }
3211 3211
3212 float erfinv (float x) 3212 float erfinv (float x)
3213 { 3213 {
3214 return do_erfinv (x, false); 3214 return do_erfinv (x, false);
3215 } 3215 }
3216 3216
3217 // Implementation based on the Fortran code by W.J.Cody 3217 // Implementation based on the Fortran code by W.J.Cody
3218 // see http://www.netlib.org/specfun/erf. 3218 // see http://www.netlib.org/specfun/erf.
3219 // Templatized and simplified workflow. 3219 // Templatized and simplified workflow.
3220 3220
3221 // FIXME: Maybe this should be globally visible. 3221 // FIXME: Maybe this should be globally visible.
3222 static inline float erfc (float x) { return erfcf (x); } 3222 static inline float erfc (float x) { return erfcf (x); }
3223 3223
3224 template <class T> 3224 template <class T>
3225 static T 3225 static T
3226 erfcx_impl (T x) 3226 erfcx_impl (T x)
3227 { 3227 {
3228 static const T c[] = 3228 static const T c[] =
3229 { 3229 {
3230 5.64188496988670089e-1,8.88314979438837594, 3230 5.64188496988670089e-1,8.88314979438837594,
3231 6.61191906371416295e+1,2.98635138197400131e+2, 3231 6.61191906371416295e+1,2.98635138197400131e+2,
3232 8.81952221241769090e+2,1.71204761263407058e+3, 3232 8.81952221241769090e+2,1.71204761263407058e+3,
3233 2.05107837782607147e+3,1.23033935479799725e+3, 3233 2.05107837782607147e+3,1.23033935479799725e+3,
3234 2.15311535474403846e-8 3234 2.15311535474403846e-8
3235 }; 3235 };
3236 3236
3237 static const T d[] = 3237 static const T d[] =
3238 { 3238 {
3239 1.57449261107098347e+1,1.17693950891312499e+2, 3239 1.57449261107098347e+1,1.17693950891312499e+2,
3240 5.37181101862009858e+2,1.62138957456669019e+3, 3240 5.37181101862009858e+2,1.62138957456669019e+3,
3241 3.29079923573345963e+3,4.36261909014324716e+3, 3241 3.29079923573345963e+3,4.36261909014324716e+3,
3242 3.43936767414372164e+3,1.23033935480374942e+3 3242 3.43936767414372164e+3,1.23033935480374942e+3
3243 }; 3243 };
3244 3244
3245 static const T p[] = 3245 static const T p[] =
3246 { 3246 {
3247 3.05326634961232344e-1,3.60344899949804439e-1, 3247 3.05326634961232344e-1,3.60344899949804439e-1,
3248 1.25781726111229246e-1,1.60837851487422766e-2, 3248 1.25781726111229246e-1,1.60837851487422766e-2,
3249 6.58749161529837803e-4,1.63153871373020978e-2 3249 6.58749161529837803e-4,1.63153871373020978e-2
3250 }; 3250 };
3251 3251
3252 static const T q[] = 3252 static const T q[] =
3253 { 3253 {
3254 2.56852019228982242,1.87295284992346047, 3254 2.56852019228982242,1.87295284992346047,
3255 5.27905102951428412e-1,6.05183413124413191e-2, 3255 5.27905102951428412e-1,6.05183413124413191e-2,
3256 2.33520497626869185e-3 3256 2.33520497626869185e-3
3257 }; 3257 };