Mercurial > hg > octave-lyh
diff liboctave/lo-specfun.cc @ 4911:14027e0bafa4
[project @ 2004-07-22 19:58:06 by jwe]
| author | jwe |
|---|---|
| date | Thu, 22 Jul 2004 19:58:06 +0000 |
| parents | 9f7ef92b50b0 |
| children | 23b37da9fd5b |
line wrap: on
line diff
--- a/liboctave/lo-specfun.cc +++ b/liboctave/lo-specfun.cc @@ -223,6 +223,12 @@ return retval; } +static inline bool +is_integer_value (double x) +{ + return x == static_cast<double> (static_cast<long> (x)); +} + static inline Complex zbesj (const Complex& z, double alpha, int kode, int& ierr) { @@ -252,6 +258,15 @@ retval = bessel_return_value (Complex (yr, yi), ierr); } + else if (is_integer_value (alpha)) + { + // zbesy can overflow as z->0, and cause troubles for generic case below + alpha = -alpha; + Complex tmp = zbesj (z, alpha, kode, ierr); + if ((static_cast <long> (alpha)) & 1) + tmp = - tmp; + retval = bessel_return_value (tmp, ierr); + } else { alpha = -alpha; @@ -313,6 +328,15 @@ return bessel_return_value (Complex (yr, yi), ierr); } + else if (is_integer_value (alpha - 0.5)) + { + // zbesy can overflow as z->0, and cause troubles for generic case below + alpha = -alpha; + Complex tmp = zbesj (z, alpha, kode, ierr); + if ((static_cast <long> (alpha - 0.5)) & 1) + tmp = - tmp; + retval = bessel_return_value (tmp, ierr); + } else { alpha = -alpha; @@ -369,8 +393,9 @@ if (ierr == 0 || ierr == 3) { - tmp += (2.0 / M_PI) * sin (M_PI * alpha) - * zbesk (z, alpha, kode, ierr); + if (! is_integer_value (alpha - 0.5)) + tmp += (2.0 / M_PI) * sin (M_PI * alpha) + * zbesk (z, alpha, kode, ierr); retval = bessel_return_value (tmp, ierr); }