Mercurial > hg > octave-lyh

comparison liboctave/lo-specfun.cc @ 4911:14027e0bafa4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2004-07-22 19:58:06 by jwe]
author jwe
date Thu, 22 Jul 2004 19:58:06 +0000
parents 9f7ef92b50b0
children 23b37da9fd5b
comparison
equal deleted inserted replaced
4910:1242acab4246 4911:14027e0bafa4
221 } 221 }
222 222
223 return retval; 223 return retval;
224 } 224 }
225 225
226 static inline bool
227 is_integer_value (double x)
228 {
229 return x == static_cast<double> (static_cast<long> (x));
230 }
231
226 static inline Complex 232 static inline Complex
227 zbesj (const Complex& z, double alpha, int kode, int& ierr) 233 zbesj (const Complex& z, double alpha, int kode, int& ierr)
228 { 234 {
229 Complex retval; 235 Complex retval;
230 236
249 255
250 if (zi == 0.0 && zr >= 0.0) 256 if (zi == 0.0 && zr >= 0.0)
251 yi = 0.0; 257 yi = 0.0;
252 258
253 retval = bessel_return_value (Complex (yr, yi), ierr); 259 retval = bessel_return_value (Complex (yr, yi), ierr);
260 }
261 else if (is_integer_value (alpha))
262 {
263 // zbesy can overflow as z->0, and cause troubles for generic case below
264 alpha = -alpha;
265 Complex tmp = zbesj (z, alpha, kode, ierr);
266 if ((static_cast <long> (alpha)) & 1)
267 tmp = - tmp;
268 retval = bessel_return_value (tmp, ierr);
254 } 269 }
255 else 270 else
256 { 271 {
257 alpha = -alpha; 272 alpha = -alpha;
258 273
311 yi = 0.0; 326 yi = 0.0;
312 } 327 }
313 328
314 return bessel_return_value (Complex (yr, yi), ierr); 329 return bessel_return_value (Complex (yr, yi), ierr);
315 } 330 }
331 else if (is_integer_value (alpha - 0.5))
332 {
333 // zbesy can overflow as z->0, and cause troubles for generic case below
334 alpha = -alpha;
335 Complex tmp = zbesj (z, alpha, kode, ierr);
336 if ((static_cast <long> (alpha - 0.5)) & 1)
337 tmp = - tmp;
338 retval = bessel_return_value (tmp, ierr);
339 }
316 else 340 else
317 { 341 {
318 alpha = -alpha; 342 alpha = -alpha;
319 343
320 Complex tmp = cos (M_PI * alpha) * zbesy (z, alpha, kode, ierr); 344 Complex tmp = cos (M_PI * alpha) * zbesy (z, alpha, kode, ierr);
367 391
368 Complex tmp = zbesi (z, alpha, kode, ierr); 392 Complex tmp = zbesi (z, alpha, kode, ierr);
369 393
370 if (ierr == 0 || ierr == 3) 394 if (ierr == 0 || ierr == 3)
371 { 395 {
372 tmp += (2.0 / M_PI) * sin (M_PI * alpha) 396 if (! is_integer_value (alpha - 0.5))
373 * zbesk (z, alpha, kode, ierr); 397 tmp += (2.0 / M_PI) * sin (M_PI * alpha)
398 * zbesk (z, alpha, kode, ierr);
374 399
375 retval = bessel_return_value (tmp, ierr); 400 retval = bessel_return_value (tmp, ierr);
376 } 401 }
377 else 402 else
378 retval = Complex (octave_NaN, octave_NaN); 403 retval = Complex (octave_NaN, octave_NaN);