[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,
>=20
> On Thu, Feb 18, 2010 at 12:31 PM, Kenta Murata <[e-mail]> =
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.
>>=20
>> Do you think what does Math.lgamma(-Float::INFINITY) return?
>=20
> 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:
>>=20
>> (1) Pick up domain and pole error cases before calling actual C =
functions.
>> (2) Returns NAN if a domain error is picked up.
>=20
> 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.
>=20
> 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
>> ...
>=20
> The patch posted is a great step in the right direction.
>=20
> 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) =3D=3D -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;
=20
+RUBY_EXTERN VALUE rb_eMathDomainError;
+
 RUBY_EXTERN VALUE rb_stdin, rb_stdout, rb_stderr;
=20
 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]))
=20
 VALUE rb_mMath;
+VALUE rb_eMathDomainError;
=20
 extern VALUE rb_to_float(VALUE val);
 #define Need_Float(x) do {if (TYPE(x) !=3D T_FLOAT) {(x) =3D =
rb_to_float(x);}} while(0)
@@ -24,6 +25,9 @@ extern VALUE rb_to_float(VALUE val);
     Need_Float(y);\
 } while (0)
=20
+#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 =3D RFLOAT_VALUE(x);
+    dy =3D RFLOAT_VALUE(y);
+    if (dx =3D=3D 0.0 && dy =3D=3D 0.0) domain_error("atan2");
+    return DBL2NUM(atan2(dy, dx));
 }
=20
=20
@@ -153,6 +161,8 @@ math_acos(VALUE obj, VALUE x)
     Need_Float(x);
     errno =3D 0;
     d0 =3D RFLOAT_VALUE(x);
+    /* check for domain error */
+    if (d0 < -1.0 || 1.0 < d0) domain_error("acos");
     d =3D acos(d0);
     domain_check(d0, d, "acos");
     return DBL2NUM(d);
@@ -173,6 +183,8 @@ math_asin(VALUE obj, VALUE x)
     Need_Float(x);
     errno =3D 0;
     d0 =3D RFLOAT_VALUE(x);
+    /* check for domain error */
+    if (d0 < -1.0 || 1.0 < d0) domain_error("asin");
     d =3D asin(d0);
     domain_check(d0, d, "asin");
     return DBL2NUM(d);
@@ -276,6 +288,8 @@ math_acosh(VALUE obj, VALUE x)
     Need_Float(x);
     errno =3D 0;
     d0 =3D RFLOAT_VALUE(x);
+    /* check for domain error */
+    if (d0 < 1.0) domain_error("acosh");
     d =3D acosh(d0);
     domain_check(d0, d, "acosh");
     return DBL2NUM(d);
@@ -310,10 +324,11 @@ math_atanh(VALUE obj, VALUE x)
     Need_Float(x);
     errno =3D 0;
     d0 =3D RFLOAT_VALUE(x);
-    if (d0 =3D=3D 1.0 || d0 =3D=3D -1.0) {
-	errno =3D 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 =3D=3D -1.0) return DBL2NUM(-INFINITY);
+    if (d0 =3D=3D +1.0) return DBL2NUM(+INFINITY);
     d =3D 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 =3D 0;
     d0 =3D RFLOAT_VALUE(x);
+    /* check for domain error */
+    if (signbit(d0)) domain_error("log");
+    /* check for pole error */
+    if (d0 =3D=3D +0.0) return DBL2NUM(-INFINITY);
     d =3D log(d0);
     if (argc =3D=3D 2) {
 	Need_Float(base);
@@ -417,6 +436,10 @@ math_log2(VALUE obj, VALUE x)
     Need_Float(x);
     errno =3D 0;
     d0 =3D RFLOAT_VALUE(x);
+    /* check for domain error */
+    if (signbit(d0)) domain_error("log2");
+    /* check for pole error */
+    if (d0 =3D=3D +0.0) return DBL2NUM(-INFINITY);
     d =3D 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 =3D 0;
     d0 =3D RFLOAT_VALUE(x);
+    /* check for domain error */
+    if (signbit(d0)) domain_error("log10");
+    /* check for pole error */
+    if (d0 =3D=3D +0.0) return DBL2NUM(-INFINITY);
     d =3D 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 =3D 0;
     d0 =3D RFLOAT_VALUE(x);
+    /* check for domain error */
+    if (signbit(d0)) domain_error("sqrt");
     d =3D 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 =3D RFLOAT_VALUE(x);
+    /* check for domain error */
+    if (isinf(d0) && signbit(d0)) domain_error("gamma");
     fracpart =3D modf(d0, &intpart);
     if (fracpart =3D=3D 0.0 &&
         0 < intpart &&
@@ -719,6 +750,8 @@ math_lgamma(VALUE obj, VALUE x)
     Need_Float(x);
     errno =3D 0;
     d0 =3D 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 =3D rb_define_module("Math");
+    rb_eMathDomainError =3D rb_define_class_under(rb_mMath, =
"DomainError", rb_eArgError);
=20
 #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'
=20
 class TestMath < Test::Unit::TestCase
   def assert_infinity(a, *rest)
-    rest =3D ["not infinity"] if rest.empty?
+    rest =3D ["not infinity: #{a.inspect}"] if rest.empty?
     assert(!a.finite?, *rest)
   end
=20
+  def assert_nan(a, *rest)
+    rest =3D ["not nan: #{a.inspect}"] if rest.empty?
+    assert(a.nan?, *rest)
+  end
+
   def check(a, b)
     err =3D [Float::EPSILON * 4, [a.abs, b.abs].max * Float::EPSILON * =
256].max
     assert_in_delta(a, b, err)
   end
=20
   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
=20
   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
=20
   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
=20
   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
=20
   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
=20
@@ -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
=20
   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
=20
   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
=20
   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
=20
   def test_lgamma
@@ -241,6 +261,8 @@ class TestMath < Test::Unit::TestCase
     g, s =3D Math.lgamma(4)
     check(Math.log(6), g)
     assert_equal(s, 1)
+
+    assert_raise(Math::DomainError) { Math.lgamma(-Float::INFINITY) }
   end
=20
   def test_cbrt

--
Kenta Murata
OpenPGP FP =3D FA26 35D7 4F98 3498 0810 E0D5 F213 966F E9EB 0BCC

E-mail: [e-mail]
twitter: http://twitter.com/mrkn/
blog: http://d.hatena.ne.jp/mrkn/