Mercurial > hg > octave-nkf

comparison liboctave/numeric/lo-specfun.cc @ 19552:c6437824681c stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improve Matlab compatibility for gamma function (bug #43551) * lo-specfun.cc (xgamma): Return Inf instead of NaN for negative integer arguments. * mappers.cc: Fix tests for gamma.
author John W. Eaton <jwe@octave.org>
date Thu, 20 Nov 2014 15:48:35 -0500
parents 49a5a4be04a1
children 998628b7963a
comparison
equal deleted inserted replaced
19544:043440fa7006 19552:c6437824681c
364 double 364 double
365 xgamma (double x) 365 xgamma (double x)
366 { 366 {
367 double result; 367 double result;
368 368
369 if (xisnan (x) || (x < 0 && (xisinf (x) || D_NINT (x) == x))) 369 // Special cases for (near) compatibility with Matlab instead of
370 // tgamma. Matlab does not have -0.
371
372 if (x == 0)
373 result = xnegative_sign (x) ? -octave_Inf : octave_Inf;
374 else if ((x < 0 && D_NINT (x) == x) || xisinf (x))
375 result = octave_Inf;
376 else if (xisnan (x))
370 result = octave_NaN; 377 result = octave_NaN;
371 else if (x == 0 && xnegative_sign (x))
372 result = -octave_Inf;
373 else if (x == 0 || xisinf (x))
374 result = octave_Inf;
375 else 378 else
376 { 379 {
377 #if defined (HAVE_TGAMMA) 380 #if defined (HAVE_TGAMMA)
378 result = tgamma (x); 381 result = tgamma (x);
379 #else 382 #else
380 F77_XFCN (xdgamma, XDGAMMA, (x, result)); 383 F77_XFCN (xdgamma, XDGAMMA, (x, result));
381 #endif 384 #endif
382 if (xisinf (result) && (static_cast<int> (gnulib::floor (x)) % 2))
383 result = -octave_Inf;
384 } 385 }
385 386
386 return result; 387 return result;
387 } 388 }
388 389
435 float 436 float
436 xgamma (float x) 437 xgamma (float x)
437 { 438 {
438 float result; 439 float result;
439 440
440 if (xisnan (x) || (x < 0 && (xisinf (x) || D_NINT (x) == x))) 441 // Special cases for (near) compatibility with Matlab instead of
442 // tgamma. Matlab does not have -0.
443
444 if (x == 0)
445 result = xnegative_sign (x) ? -octave_Float_Inf : octave_Float_Inf;
446 else if ((x < 0 && D_NINT (x) == x) || xisinf (x))
447 result = octave_Float_Inf;
448 else if (xisnan (x))
441 result = octave_Float_NaN; 449 result = octave_Float_NaN;
442 else if (x == 0 && xnegative_sign (x))
443 result = -octave_Float_Inf;
444 else if (x == 0 || xisinf (x))
445 result = octave_Float_Inf;
446 else 450 else
447 { 451 {
448 #if defined (HAVE_TGAMMA) 452 #if defined (HAVE_TGAMMA)
449 result = tgammaf (x); 453 result = tgammaf (x);
450 #else 454 #else
451 F77_XFCN (xgamma, XGAMMA, (x, result)); 455 F77_XFCN (xgamma, XGAMMA, (x, result));
452 #endif 456 #endif
453 if (xisinf (result) && (static_cast<int> (gnulib::floor (x)) % 2))
454 result = -octave_Float_Inf;
455 } 457 }
456 458
457 return result; 459 return result;
458 } 460 }
459 461