diff --git a/Eigen/src/Core/SpecialFunctions.h b/Eigen/src/Core/SpecialFunctions.h index 4a61325d4..567c02c61 100644 --- a/Eigen/src/Core/SpecialFunctions.h +++ b/Eigen/src/Core/SpecialFunctions.h @@ -520,14 +520,16 @@ struct igammac_impl { Copyright 1985, 1987, 1992 by Stephen L. Moshier Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */ + using std::log; const Scalar zero = 0; const Scalar one = 1; const Scalar two = 2; const Scalar machep = igamma_helper::machep(); - const Scalar maxlog = ::log(NumTraits::highest()); + const Scalar maxlog = log(NumTraits::highest()); const Scalar big = igamma_helper::big(); const Scalar biginv = 1 / big; const Scalar nan = NumTraits::quiet_NaN(); + const Scalar inf = NumTraits::infinity(); Scalar ans, ax, c, yc, r, t, y, z; Scalar pk, pkm1, pkm2, qk, qkm1, qkm2; @@ -541,11 +543,9 @@ struct igammac_impl { return (one - igamma_impl::run(a, x)); } - using std::isinf; - if ((isinf)(x)) return zero; + if (x == inf) return zero; // std::isinf crashes on CUDA /* Compute x**a * exp(-x) / gamma(a) */ - using std::log; ax = a * log(x) - x - lgamma_impl::run(a); if (ax < -maxlog) { // underflow return zero; @@ -564,7 +564,7 @@ struct igammac_impl { ans = pkm1 / qkm1; using std::abs; - do { + while (true) { c += one; y += one; z += two; @@ -588,7 +588,8 @@ struct igammac_impl { qkm2 *= biginv; qkm1 *= biginv; } - } while (t > machep); + if (t <= machep) break; + } return (ans * ax); } @@ -683,10 +684,11 @@ struct igamma_impl { * k=0 | (a+k+1) * */ + using std::log; const Scalar zero = 0; const Scalar one = 1; const Scalar machep = igamma_helper::machep(); - const Scalar maxlog = ::log(NumTraits::highest()); + const Scalar maxlog = log(NumTraits::highest()); const Scalar nan = NumTraits::quiet_NaN(); double ans, ax, c, r; @@ -702,7 +704,6 @@ struct igamma_impl { } /* Compute x**a * exp(-x) / gamma(a) */ - using std::log; ax = a * log(x) - x - lgamma_impl::run(a); if (ax < -maxlog) { // underflow @@ -716,11 +717,12 @@ struct igamma_impl { c = one; ans = one; - do { + while (true) { r += one; c *= x/r; ans += c; - } while (c/ans > machep); + if (c/ans <= machep) break; + } return (ans * ax / a); } diff --git a/unsupported/test/cxx11_tensor_cuda.cu b/unsupported/test/cxx11_tensor_cuda.cu index 348271e4b..1964d9e07 100644 --- a/unsupported/test/cxx11_tensor_cuda.cu +++ b/unsupported/test/cxx11_tensor_cuda.cu @@ -632,7 +632,6 @@ void test_cuda_igamma() Tensor a(6, 6); Tensor x(6, 6); Tensor out(6, 6); - Tensor expected_out(6, 6); out.setZero(); Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; @@ -641,7 +640,7 @@ void test_cuda_igamma() for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) { a(i, j) = a_s[i]; - x(i, j) = x_s[i]; + x(i, j) = x_s[j]; } } @@ -686,8 +685,7 @@ void test_cuda_igamma() for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) { if ((std::isnan)(igamma_s[i][j])) { - printf("got: %g\n", out(i, j)); - //VERIFY((std::isnan)(out(i, j))); + VERIFY((std::isnan)(out(i, j))); } else { VERIFY_IS_APPROX(out(i, j), igamma_s[i][j]); } @@ -701,7 +699,6 @@ void test_cuda_igammac() Tensor a(6, 6); Tensor x(6, 6); Tensor out(6, 6); - Tensor expected_out(6, 6); out.setZero(); Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; @@ -710,7 +707,7 @@ void test_cuda_igammac() for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) { a(i, j) = a_s[i]; - x(i, j) = x_s[i]; + x(i, j) = x_s[j]; } } @@ -754,8 +751,7 @@ void test_cuda_igammac() for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j) { if ((std::isnan)(igammac_s[i][j])) { - printf("got: %g\n", out(i, j)); - //VERIFY((std::isnan)(out(i, j))); + VERIFY((std::isnan)(out(i, j))); } else { VERIFY_IS_APPROX(out(i, j), igammac_s[i][j]); }