From: Kenta Murata Date: 2010-02-22T03:32:55+09:00 Subject: [ruby-core:28286] Re: Is Math module a wrapper of libm? Hi First, I apologise to confuse you in the previous mail. I ideally hope that Math functions raise appropriate errors for their domain errors. On 2010/02/20, at 13:11, Marc-Andre Lafortune wrote: > Hi, > > On Thu, Feb 18, 2010 at 12:31 PM, Kenta Murata wrote: >> lgamma(x) is defined as log(fabs(tgamma(x))). >> SUSv3 defines that lgamma(-INFINITY) is +INFINITY >> though tgamma(-INFINITY) is defined a domain error. >> >> Do you think what does Math.lgamma(-Float::INFINITY) return? > > Indeed, it makes no mathematical sense that > Math.lgamma(-Float::INFINITY) return +INFINITY. Since lgamma > oscillates wildly for negative arguments, it should return NaN or > raise an error, same as tgamma. So I want lgamma raises an error for -INFINITY input. >> I try to write a patch for atanh, log, log2, log10, and sqrt >> according to the following policies: >> >> (1) Pick up domain and pole error cases before calling actual C functions. >> (2) Returns NAN if a domain error is picked up. > > Are you proposing to never raise an error for any domain errors and > instead returning NaN? > I am a bit confused here, since you suggest a bit afterwardsto create > a new error class Math::DomainError. As above description, I hope Math functions raise errors for their domain errors. >> (3) Returns INFINITY with suitable flag if a pole error is picked up. > > Returning +-INFINITY is great. What "suitable flag" did you have in mind? I made a mistake, "flag" is not right but "sign". >> ... >> diff --git a/math.c b/math.c >> ... > > The patch posted is a great step in the right direction. > > Some care has to be taken with -0.0. I completely understand why you > feel that Log(-0.0) should return NaN. IEEE decided it was wiser to > return -Infinity, though. It can be argued that having an actual value > that could be exact is more useful than a NaN. On the other hand, > don't ask me why IEEE decided that Sqrt(-0.0) == -0.0 instead of > +0.0, I have no clue! When I wrote the previous mail, I thought that -0.0 includes negative underflowed values, so I defined Log(-0.0) as NaN. But my thought isn't useful for users who are unfamiliar with characteristics of float values. It is a difficult issue. I think it needs more discussions. Now I refine the previous patch. Math::DomainError is introduced now and used by the following cases: - Math.acos(|x| > 1) - Math.acosh(x < 1) - Math.asin(|x| > 1) - Math.atan2(+-0, +-0) # needs more discussions - Math.atanh(|x| > 1) - Math.log(x < 0) - Math.log2(x < 0) - Math.log10(x < 0) - Math.log(-0.0) # needs more discussions - Math.log2(-0.0) # needs more discussions - Math.log10(-0.0) # needs more discussions - Math.sqrt(x < 0) - Math.sqrt(-0.0) # needs more discussions - Math.gamma(-inf) - Math.lgamma(-inf) Please review it. diff --git a/include/ruby/ruby.h b/include/ruby/ruby.h index 5d87c0d..d56d42c 100644 --- a/include/ruby/ruby.h +++ b/include/ruby/ruby.h @@ -1234,6 +1234,8 @@ RUBY_EXTERN VALUE rb_eNameError; RUBY_EXTERN VALUE rb_eSyntaxError; RUBY_EXTERN VALUE rb_eLoadError; +RUBY_EXTERN VALUE rb_eMathDomainError; + RUBY_EXTERN VALUE rb_stdin, rb_stdout, rb_stderr; static inline VALUE diff --git a/math.c b/math.c index 51a89c3..5064d05 100644 --- a/math.c +++ b/math.c @@ -16,6 +16,7 @@ #define numberof(array) (int)(sizeof(array) / sizeof((array)[0])) VALUE rb_mMath; +VALUE rb_eMathDomainError; extern VALUE rb_to_float(VALUE val); #define Need_Float(x) do {if (TYPE(x) != T_FLOAT) {(x) = rb_to_float(x);}} while(0) @@ -24,6 +25,9 @@ extern VALUE rb_to_float(VALUE val); Need_Float(y);\ } while (0) +#define domain_error(msg) \ + rb_raise(rb_eMathDomainError, "Numerical argument is out of domain - " #msg); + static void domain_check(double x, double y, const char *msg) { @@ -86,8 +90,12 @@ infinity_check(VALUE arg, double res, const char *msg) static VALUE math_atan2(VALUE obj, VALUE y, VALUE x) { + double dx, dy; Need_Float2(y, x); - return DBL2NUM(atan2(RFLOAT_VALUE(y), RFLOAT_VALUE(x))); + dx = RFLOAT_VALUE(x); + dy = RFLOAT_VALUE(y); + if (dx == 0.0 && dy == 0.0) domain_error("atan2"); + return DBL2NUM(atan2(dy, dx)); } @@ -153,6 +161,8 @@ math_acos(VALUE obj, VALUE x) Need_Float(x); errno = 0; d0 = RFLOAT_VALUE(x); + /* check for domain error */ + if (d0 < -1.0 || 1.0 < d0) domain_error("acos"); d = acos(d0); domain_check(d0, d, "acos"); return DBL2NUM(d); @@ -173,6 +183,8 @@ math_asin(VALUE obj, VALUE x) Need_Float(x); errno = 0; d0 = RFLOAT_VALUE(x); + /* check for domain error */ + if (d0 < -1.0 || 1.0 < d0) domain_error("asin"); d = asin(d0); domain_check(d0, d, "asin"); return DBL2NUM(d); @@ -276,6 +288,8 @@ math_acosh(VALUE obj, VALUE x) Need_Float(x); errno = 0; d0 = RFLOAT_VALUE(x); + /* check for domain error */ + if (d0 < 1.0) domain_error("acosh"); d = acosh(d0); domain_check(d0, d, "acosh"); return DBL2NUM(d); @@ -310,10 +324,11 @@ math_atanh(VALUE obj, VALUE x) Need_Float(x); errno = 0; d0 = RFLOAT_VALUE(x); - if (d0 == 1.0 || d0 == -1.0) { - errno = ERANGE; - rb_sys_fail("atanh"); - } + /* check for domain error */ + if (d0 < -1.0 || +1.0 < d0) domain_error("atanh"); + /* check for pole error */ + if (d0 == -1.0) return DBL2NUM(-INFINITY); + if (d0 == +1.0) return DBL2NUM(+INFINITY); d = atanh(d0); domain_check(d0, d, "atanh"); infinity_check(x, d, "atanh"); @@ -374,6 +389,10 @@ math_log(int argc, VALUE *argv) Need_Float(x); errno = 0; d0 = RFLOAT_VALUE(x); + /* check for domain error */ + if (signbit(d0)) domain_error("log"); + /* check for pole error */ + if (d0 == +0.0) return DBL2NUM(-INFINITY); d = log(d0); if (argc == 2) { Need_Float(base); @@ -417,6 +436,10 @@ math_log2(VALUE obj, VALUE x) Need_Float(x); errno = 0; d0 = RFLOAT_VALUE(x); + /* check for domain error */ + if (signbit(d0)) domain_error("log2"); + /* check for pole error */ + if (d0 == +0.0) return DBL2NUM(-INFINITY); d = log2(d0); domain_check(d0, d, "log2"); infinity_check(x, d, "log2"); @@ -443,6 +466,10 @@ math_log10(VALUE obj, VALUE x) Need_Float(x); errno = 0; d0 = RFLOAT_VALUE(x); + /* check for domain error */ + if (signbit(d0)) domain_error("log10"); + /* check for pole error */ + if (d0 == +0.0) return DBL2NUM(-INFINITY); d = log10(d0); domain_check(d0, d, "log10"); infinity_check(x, d, "log10"); @@ -481,6 +508,8 @@ math_sqrt(VALUE obj, VALUE x) Need_Float(x); errno = 0; d0 = RFLOAT_VALUE(x); + /* check for domain error */ + if (signbit(d0)) domain_error("sqrt"); d = sqrt(d0); domain_check(d0, d, "sqrt"); return DBL2NUM(d); @@ -686,6 +715,8 @@ math_gamma(VALUE obj, VALUE x) double intpart, fracpart; Need_Float(x); d0 = RFLOAT_VALUE(x); + /* check for domain error */ + if (isinf(d0) && signbit(d0)) domain_error("gamma"); fracpart = modf(d0, &intpart); if (fracpart == 0.0 && 0 < intpart && @@ -719,6 +750,8 @@ math_lgamma(VALUE obj, VALUE x) Need_Float(x); errno = 0; d0 = RFLOAT_VALUE(x); + /* check for domain error */ + if (isinf(d0) && signbit(d0)) domain_error("lgamma"); if (isinf(d0)) { return rb_assoc_new(DBL2NUM(INFINITY), INT2FIX(1)); } @@ -772,6 +805,7 @@ void Init_Math(void) { rb_mMath = rb_define_module("Math"); + rb_eMathDomainError = rb_define_class_under(rb_mMath, "DomainError", rb_eArgError); #ifdef M_PI rb_define_const(rb_mMath, "PI", DBL2NUM(M_PI)); diff --git a/test/ruby/test_math.rb b/test/ruby/test_math.rb index 5b49291..7270511 100644 --- a/test/ruby/test_math.rb +++ b/test/ruby/test_math.rb @@ -2,16 +2,22 @@ require 'test/unit' class TestMath < Test::Unit::TestCase def assert_infinity(a, *rest) - rest = ["not infinity"] if rest.empty? + rest = ["not infinity: #{a.inspect}"] if rest.empty? assert(!a.finite?, *rest) end + def assert_nan(a, *rest) + rest = ["not nan: #{a.inspect}"] if rest.empty? + assert(a.nan?, *rest) + end + def check(a, b) err = [Float::EPSILON * 4, [a.abs, b.abs].max * Float::EPSILON * 256].max assert_in_delta(a, b, err) end def test_atan2 + assert_raise(Math::DomainError) { Math.atan2(0, 0) } check(0, Math.atan2(0, 1)) check(Math::PI / 4, Math.atan2(1, 1)) check(Math::PI / 2, Math.atan2(1, 0)) @@ -46,7 +52,9 @@ class TestMath < Test::Unit::TestCase check(1 * Math::PI / 4, Math.acos( 1.0 / Math.sqrt(2))) check(2 * Math::PI / 4, Math.acos( 0.0)) check(4 * Math::PI / 4, Math.acos(-1.0)) - assert_raise(Errno::EDOM, Errno::ERANGE) { Math.acos(2.0) } + assert_raise(Math::DomainError) { Math.acos(+1.0 + Float::EPSILON) } + assert_raise(Math::DomainError) { Math.acos(-1.0 - Float::EPSILON) } + assert_raise(Math::DomainError) { Math.acos(2.0) } end def test_asin @@ -54,7 +62,9 @@ class TestMath < Test::Unit::TestCase check( 1 * Math::PI / 4, Math.asin( 1.0 / Math.sqrt(2))) check( 2 * Math::PI / 4, Math.asin( 1.0)) check(-2 * Math::PI / 4, Math.asin(-1.0)) - assert_raise(Errno::EDOM, Errno::ERANGE) { Math.asin(2.0) } + assert_raise(Math::DomainError) { Math.asin(+1.0 + Float::EPSILON) } + assert_raise(Math::DomainError) { Math.asin(-1.0 - Float::EPSILON) } + assert_raise(Math::DomainError) { Math.asin(2.0) } end def test_atan @@ -86,7 +96,8 @@ class TestMath < Test::Unit::TestCase check(0, Math.acosh(1)) check(1, Math.acosh((Math::E ** 1 + Math::E ** -1) / 2)) check(2, Math.acosh((Math::E ** 2 + Math::E ** -2) / 2)) - assert_raise(Errno::EDOM, Errno::ERANGE) { Math.acosh(0) } + assert_raise(Math::DomainError) { Math.acosh(1.0 - Float::EPSILON) } + assert_raise(Math::DomainError) { Math.acosh(0) } end def test_asinh @@ -99,7 +110,10 @@ class TestMath < Test::Unit::TestCase check(0, Math.atanh(Math.sinh(0) / Math.cosh(0))) check(1, Math.atanh(Math.sinh(1) / Math.cosh(1))) check(2, Math.atanh(Math.sinh(2) / Math.cosh(2))) - assert_raise(Errno::EDOM, Errno::ERANGE) { Math.atanh(-1) } + assert_nothing_raised { assert_infinity Math.atanh(1) } + assert_nothing_raised { assert_infinity -Math.atanh(-1) } + assert_raise(Math::DomainError) { Math.atanh(+1.0 + Float::EPSILON) } + assert_raise(Math::DomainError) { Math.atanh(-1.0 - Float::EPSILON) } end def test_exp @@ -116,8 +130,9 @@ class TestMath < Test::Unit::TestCase check(1, Math.log(10, 10)) check(2, Math.log(100, 10)) assert_equal(1.0/0, Math.log(1.0/0)) - assert_raise(Errno::EDOM, Errno::ERANGE) { Math.log(0) } - assert_raise(Errno::EDOM, Errno::ERANGE) { Math.log(-1) } + assert_nothing_raised { assert_infinity -Math.log(+0.0) } + assert_raise(Math::DomainError) { Math.log(-0.0) } # TODO: [ruby-core:28265] + assert_raise(Math::DomainError) { Math.log(-1.0) } assert_raise(TypeError) { Math.log(1,nil) } end @@ -126,8 +141,9 @@ class TestMath < Test::Unit::TestCase check(1, Math.log2(2)) check(2, Math.log2(4)) assert_equal(1.0/0, Math.log2(1.0/0)) - assert_raise(Errno::EDOM, Errno::ERANGE) { Math.log2(0) } - assert_raise(Errno::EDOM, Errno::ERANGE) { Math.log2(-1) } + assert_nothing_raised { assert_infinity -Math.log2(+0.0) } + assert_raise(Math::DomainError) { Math.log2(-0.0) } # TODO: [ruby-core:28265] + assert_raise(Math::DomainError) { Math.log2(-1.0) } end def test_log10 @@ -135,8 +151,9 @@ class TestMath < Test::Unit::TestCase check(1, Math.log10(10)) check(2, Math.log10(100)) assert_equal(1.0/0, Math.log10(1.0/0)) - assert_raise(Errno::EDOM, Errno::ERANGE) { Math.log10(0) } - assert_raise(Errno::EDOM, Errno::ERANGE) { Math.log10(-1) } + assert_nothing_raised { assert_infinity -Math.log10(+0.0) } + assert_raise(Math::DomainError) { Math.log10(-0.0) } # TODO: [ruby-core:28265] + assert_raise(Math::DomainError) { Math.log10(-1.0) } end def test_sqrt @@ -144,7 +161,8 @@ class TestMath < Test::Unit::TestCase check(1, Math.sqrt(1)) check(2, Math.sqrt(4)) assert_equal(1.0/0, Math.sqrt(1.0/0)) - assert_raise(Errno::EDOM, Errno::ERANGE) { Math.sqrt(-1) } + assert_raise(Math::DomainError) { Math.sqrt(-0.0) } # TODO: [ruby-core:28265] + assert_raise(Math::DomainError) { Math.sqrt(-1.0) } end def test_frexp @@ -201,6 +219,8 @@ class TestMath < Test::Unit::TestCase assert_infinity(Math.gamma(i), "Math.gamma(#{i}) should be INF") assert_infinity(Math.gamma(i-1), "Math.gamma(#{i-1}) should be INF") end + + assert_raise(Math::DomainError) { Math.gamma(-Float::INFINITY) } end def test_lgamma @@ -241,6 +261,8 @@ class TestMath < Test::Unit::TestCase g, s = Math.lgamma(4) check(Math.log(6), g) assert_equal(s, 1) + + assert_raise(Math::DomainError) { Math.lgamma(-Float::INFINITY) } end def test_cbrt -- Kenta Murata OpenPGP FP = FA26 35D7 4F98 3498 0810 E0D5 F213 966F E9EB 0BCC E-mail: mrkn@mrkn.jp twitter: http://twitter.com/mrkn/ blog: http://d.hatena.ne.jp/mrkn/