From: Kenta Murata <muraken@...>
Date: 2010-02-19T02:31:28+09:00
Subject: [ruby-core:28233] Re: Is Math module a wrapper of libm?

Hi,

On 2010/02/19, at 0:25, Marc-Andre Lafortune wrote:

> On Thu, Feb 18, 2010 at 7:25 AM, Yusuke ENDOH <mame@tsg.ne.jp> wrote:
>> At first, I'll soon change the behaviors of atanh(1) and lgamma(inf)
>> which are currently platform-dependent.  According to SUSv3, the former
>> should raise ERANGE (currently raises EDOM on Linux), and the latter
>> should return inf (currently raises EDOM on darwin).
> 
> I disagree with this reading. SUSv3 states that atanh(1) should return
> Infinity and (only when the application has requested that behavior)
> set the ERANGE error flag. Ruby must choose to do one or the other (or
> else add a way to do both). Mathematically speaking, there is no doubt
> that it should return Infinity and that this is more helpful than
> raising an error.

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?

> The question is also valid for
> Math.log(0.0)
> Math.atanh(1.0)
> Math.gamma(0.0)
> 0.0**(-n)
> 
> Look at the graph of 1/x. Look at the graph of Log(x). Now explain to
> me why 1/0.0 is +Inf and Log(x) would not be -Inf?


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.
(3) Returns INFINITY with suitable flag if a pole error is picked up.

But, I hope that a suitable exception is raised when picked up domain errors.
Unfortunately the existence of EDOM depends on platform.
So I believe that we need an error class like Math::DomainError.

--
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/


diff --git a/math.c b/math.c
index 027932d..ebde179 100644
--- a/math.c
+++ b/math.c
@@ -310,6 +310,11 @@ math_atanh(VALUE obj, VALUE x)
     Need_Float(x);
     errno = 0;
     d0 = RFLOAT_VALUE(x);
+    /* domain error check */
+    if (d0 <  -1.0 || +1.0 <  d0) return DBL2NUM(NAN);
+    /* pole error cases */
+    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");
@@ -370,6 +375,10 @@ math_log(int argc, VALUE *argv)
     Need_Float(x);
     errno = 0;
     d0 = RFLOAT_VALUE(x);
+    /* domain error case */
+    if (signbit(d0)) return DBL2NUM(NAN);
+    /* pole error case */
+    if (d0 == +0.0) return DBL2NUM(-INFINITY);
     d = log(d0);
     if (argc == 2) {
 	Need_Float(base);
@@ -413,6 +422,10 @@ math_log2(VALUE obj, VALUE x)
     Need_Float(x);
     errno = 0;
     d0 = RFLOAT_VALUE(x);
+    /* domain error case */
+    if (signbit(d0)) return DBL2NUM(NAN);
+    /* pole error case */
+    if (d0 == +0.0) return DBL2NUM(-INFINITY);
     d = log2(d0);
     domain_check(d0, d, "log2");
     infinity_check(x, d, "log2");
@@ -439,6 +452,10 @@ math_log10(VALUE obj, VALUE x)
     Need_Float(x);
     errno = 0;
     d0 = RFLOAT_VALUE(x);
+    /* domain error case */
+    if (signbit(d0)) return DBL2NUM(NAN);
+    /* pole error case */
+    if (d0 == +0.0) return DBL2NUM(-INFINITY);
     d = log10(d0);
     domain_check(d0, d, "log10");
     infinity_check(x, d, "log10");
@@ -477,6 +494,8 @@ math_sqrt(VALUE obj, VALUE x)
     Need_Float(x);
     errno = 0;
     d0 = RFLOAT_VALUE(x);
+    /* domain error case */
+    if (signbit(d0)) return DBL2NUM(NAN);
     d = sqrt(d0);
     domain_check(d0, d, "sqrt");
     return DBL2NUM(d);
diff --git a/test/ruby/test_math.rb b/test/ruby/test_math.rb
index 5b49291..8189b31 100644
--- a/test/ruby/test_math.rb
+++ b/test/ruby/test_math.rb
@@ -2,10 +2,15 @@ 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)
@@ -99,7 +104,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_nothing_raised { assert_nan Math.atanh(1 + Float::EPSILON) }
+    assert_nothing_raised { assert_nan Math.atanh(-1 - Float::EPSILON) }
   end
 
   def test_exp
@@ -116,8 +124,8 @@ 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_nothing_raised { assert_nan Math.log(-0.0) }
     assert_raise(TypeError) { Math.log(1,nil) }
   end
 
@@ -126,8 +134,8 @@ 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_nothing_raised { assert_nan Math.log2(-0.0) }
   end
 
   def test_log10
@@ -135,8 +143,8 @@ 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_nothing_raised { assert_nan Math.log10(-0.0) }
   end
 
   def test_sqrt
@@ -144,7 +152,7 @@ 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_nothing_raised { assert_nan Math.sqrt(-0.0) }
   end
 
   def test_frexp