[ruby-core:28291] Re: Is Math module a wrapper of libm?

Hi!

I like the direction we're taking.

I've been thinking about a general decision algorithm. It would look like:

Try to come up with a valid mathematical value (including +-inf but
excluding NaN):
- Check all possible values (excluding NaN). If there is exactly one,
you're done.
- If there are many possible values, are they all basically equal
(e.g. -0.0 vs +0.0, or +pi vs -pi for an angle)? If so, then return
the best one.
- If the possible values are +-Infinity, and the input is +-0.0 can
tell you which to pick, then return the right infinity (pole error
case)
- Otherwise return NaN or raise an error.

This algorithm agrees with your proposal that atan2(0,0) should raise
an error. While we're contradicting IEEE, have you thought of
atan2(+-Inf, +-Inf)? I feel it should also be an error, in the same
way that Inf/Inf doesn't return 1!

Here's an incremental patch that modifies the code according to this.
Note that infinity_check was already redundant, so I removed it.
domain_check was also redundant except potentially for some corner
cases. Indeed SUSv3 states that implementations may set the ERANGE
error flag in case of subnormal values. domain_check would raise an
error if this was the case. Since this behavior is platform dependant
(and far from useful), I removed domain_check also.

I included a line so that sqrt(-0.0) = +0.0 instead of -0.0, because I
feel it makes sense and is consistent with pow(-0.0, 0.5) == +0.0, but
I don't think it matters all that much and I definitely won't fight
over it.

The question that remains in my mind: is it better to return NaN or
raise a domain error? Why should 0.0/0.0 == NaN and atan2(0.0, 0.0)
raise an error? I personally would lean towards always returning NaN
instead of raising an error, but I am not confident in this
preference. Or never return NaN? Thoughts?


diff --git a/math.c b/math.c
index 5064d05..77a8757 100644
--- a/math.c
+++ b/math.c
@@ -11,7 +11,6 @@

 #include "ruby/ruby.h"
 #include <math.h>
-#include <errno.h>

 #define numberof(array) (int)(sizeof(array) / sizeof((array)[0]))

@@ -28,45 +27,6 @@ extern VALUE rb_to_float(VALUE val);
 #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)
-{
-    if (errno) {
-	if (isinf(y)) return;
-    }
-    else {
-	if (!isnan(y)) return;
-	else if (isnan(x)) return;
-	else {
-#if defined(EDOM)
-	    errno = EDOM;
-#else
-	    errno = ERANGE;
-#endif
-	}
-    }
-    rb_sys_fail(msg);
-}
-
-static void
-infinity_check(VALUE arg, double res, const char *msg)
-{
-    while(1) {
-	if (errno) {
-	    rb_sys_fail(msg);
-	}
-	if (isinf(res) && !isinf(RFLOAT_VALUE(arg))) {
-#if defined(EDOM)
-	    errno = EDOM;
-#elif defined(ERANGE)
-	    errno = ERANGE;
-#endif
-	    continue;
-	}
-	break;
-    }
-}
-
 /*
  *  call-seq:
  *     Math.atan2(y, x)  => float
@@ -95,6 +55,7 @@ math_atan2(VALUE obj, VALUE y, VALUE x)
     dx = RFLOAT_VALUE(x);
     dy = RFLOAT_VALUE(y);
     if (dx == 0.0 && dy == 0.0) domain_error("atan2");
+    if (isinf(dx) && isinf(dy)) domain_error("atan2");
     return DBL2NUM(atan2(dy, dx));
 }

@@ -159,12 +120,10 @@ math_acos(VALUE obj, VALUE x)
     double d0, d;

     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);
 }

@@ -181,12 +140,10 @@ math_asin(VALUE obj, VALUE x)
     double d0, d;

     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);
 }

@@ -286,12 +243,10 @@ math_acosh(VALUE obj, VALUE x)
     double d0, d;

     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);
 }

@@ -322,7 +277,6 @@ math_atanh(VALUE obj, VALUE x)
     double d0, d;

     Need_Float(x);
-    errno = 0;
     d0 = RFLOAT_VALUE(x);
     /* check for domain error */
     if (d0 <  -1.0 || +1.0 <  d0) domain_error("atanh");
@@ -330,8 +284,6 @@ math_atanh(VALUE obj, VALUE x)
     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");
     return DBL2NUM(d);
 }

@@ -387,19 +339,16 @@ math_log(int argc, VALUE *argv)

     rb_scan_args(argc, argv, "11", &x, &base);
     Need_Float(x);
-    errno = 0;
     d0 = RFLOAT_VALUE(x);
     /* check for domain error */
-    if (signbit(d0)) domain_error("log");
+    if (d0 < 0.0) domain_error("log");
     /* check for pole error */
-    if (d0 == +0.0) return DBL2NUM(-INFINITY);
+    if (d0 == 0.0) return DBL2NUM(-INFINITY);
     d = log(d0);
     if (argc == 2) {
 	Need_Float(base);
 	d /= log(RFLOAT_VALUE(base));
     }
-    domain_check(d0, d, "log");
-    infinity_check(x, d, "log");
     return DBL2NUM(d);
 }

@@ -434,15 +383,12 @@ math_log2(VALUE obj, VALUE x)
     double d0, d;

     Need_Float(x);
-    errno = 0;
     d0 = RFLOAT_VALUE(x);
     /* check for domain error */
-    if (signbit(d0)) domain_error("log2");
+    if (d0 < 0.0) domain_error("log2");
     /* check for pole error */
-    if (d0 == +0.0) return DBL2NUM(-INFINITY);
+    if (d0 == 0.0) return DBL2NUM(-INFINITY);
     d = log2(d0);
-    domain_check(d0, d, "log2");
-    infinity_check(x, d, "log2");
     return DBL2NUM(d);
 }

@@ -464,15 +410,12 @@ math_log10(VALUE obj, VALUE x)
     double d0, d;

     Need_Float(x);
-    errno = 0;
     d0 = RFLOAT_VALUE(x);
     /* check for domain error */
-    if (signbit(d0)) domain_error("log10");
+    if (d0 < 0.0) domain_error("log10");
     /* check for pole error */
-    if (d0 == +0.0) return DBL2NUM(-INFINITY);
+    if (d0 == 0.0) return DBL2NUM(-INFINITY);
     d = log10(d0);
-    domain_check(d0, d, "log10");
-    infinity_check(x, d, "log10");
     return DBL2NUM(d);
 }

@@ -506,12 +449,11 @@ math_sqrt(VALUE obj, VALUE x)
     double d0, d;

     Need_Float(x);
-    errno = 0;
     d0 = RFLOAT_VALUE(x);
     /* check for domain error */
-    if (signbit(d0)) domain_error("sqrt");
+    if (d0 < 0) domain_error("sqrt");
+    if (d0 == 0) return DBL2NUM(0);
     d = sqrt(d0);
-    domain_check(d0, d, "sqrt");
     return DBL2NUM(d);
 }

@@ -718,14 +660,14 @@ math_gamma(VALUE obj, VALUE x)
     /* check for domain error */
     if (isinf(d0) && signbit(d0)) domain_error("gamma");
     fracpart = modf(d0, &intpart);
-    if (fracpart == 0.0 &&
-        0 < intpart &&
-        intpart - 1 < (double)numberof(fact_table)) {
-        return DBL2NUM(fact_table[(int)intpart - 1]);
+    if (fracpart == 0.0) {
+	if (intpart < 0) domain_error("gamma");
+        if (0 < intpart &&
+            intpart - 1 < (double)numberof(fact_table)) {
+            return DBL2NUM(fact_table[(int)intpart - 1]);
+	}
     }
-    errno = 0;
     d = tgamma(d0);
-    domain_check(d0, d, "gamma");
     return DBL2NUM(d);
 }

@@ -748,7 +690,6 @@ math_lgamma(VALUE obj, VALUE x)
     int sign=1;
     VALUE v;
     Need_Float(x);
-    errno = 0;
     d0 = RFLOAT_VALUE(x);
     /* check for domain error */
     if (isinf(d0) && signbit(d0)) domain_error("lgamma");
@@ -756,7 +697,6 @@ math_lgamma(VALUE obj, VALUE x)
 	return rb_assoc_new(DBL2NUM(INFINITY), INT2FIX(1));
     }
     d = lgamma_r(d0, &sign);
-    domain_check(d0, d, "lgamma");
     v = DBL2NUM(d);
     return rb_assoc_new(v, INT2FIX(sign));
 }
diff --git a/test/ruby/test_math.rb b/test/ruby/test_math.rb
index 3963e91..2f654c1 100644
--- a/test/ruby/test_math.rb
+++ b/test/ruby/test_math.rb
@@ -18,6 +18,7 @@ class TestMath < Test::Unit::TestCase

   def test_atan2
     assert_raise(Math::DomainError) { Math.atan2(0, 0) }
+    assert_raise(Math::DomainError) { Math.atan2(1.0 / 0.0, 1.0 / 0.0) }
     check(0, Math.atan2(0, 1))
     check(Math::PI / 4, Math.atan2(1, 1))
     check(Math::PI / 2, Math.atan2(1, 0))
@@ -131,7 +132,7 @@ class TestMath < Test::Unit::TestCase
     check(2, Math.log(100, 10))
     assert_equal(1.0/0, Math.log(1.0/0))
     assert_nothing_raised { assert_infinity -Math.log(+0.0) }
-    assert_raise(Math::DomainError) { Math.log(-0.0) } # TODO:
[ruby-core:28265]
+    assert_nothing_raised { assert_infinity -Math.log(-0.0) }
     assert_raise(Math::DomainError) { Math.log(-1.0) }
     assert_raise(TypeError) { Math.log(1,nil) }
   end
@@ -142,7 +143,7 @@ class TestMath < Test::Unit::TestCase
     check(2, Math.log2(4))
     assert_equal(1.0/0, Math.log2(1.0/0))
     assert_nothing_raised { assert_infinity -Math.log2(+0.0) }
-    assert_raise(Math::DomainError) { Math.log2(-0.0) } # TODO:
[ruby-core:28265]
+    assert_nothing_raised { assert_infinity -Math.log2(-0.0) }
     assert_raise(Math::DomainError) { Math.log2(-1.0) }
   end

@@ -152,7 +153,7 @@ class TestMath < Test::Unit::TestCase
     check(2, Math.log10(100))
     assert_equal(1.0/0, Math.log10(1.0/0))
     assert_nothing_raised { assert_infinity -Math.log10(+0.0) }
-    assert_raise(Math::DomainError) { Math.log10(-0.0) } # TODO:
[ruby-core:28265]
+    assert_nothing_raised { assert_infinity -Math.log10(-0.0) }
     assert_raise(Math::DomainError) { Math.log10(-1.0) }
   end

@@ -161,7 +162,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(Math::DomainError) { Math.sqrt(-0.0) } # TODO:
[ruby-core:28265]
+    assert_equal("0.0", Math.sqrt(-0.0).to_s) # insure it is +0.0, not -0.0
     assert_raise(Math::DomainError) { Math.sqrt(-1.0) }
   end