diff --git a/numeric.c b/numeric.c index 60c8ab6..35aa8e9 100644 --- a/numeric.c +++ b/numeric.c @@ -1480,6 +1480,119 @@ flo_is_finite_p(VALUE num) /* * call-seq: + * float.next_float -> float + * + * Returns the next representable floating-point number. + * + * Float::MAX.next_float and Float::INFINITY.next_float is Float::INFINITY. + * + * Float::NAN.next_float is Float::NAN. + * + * For example: + * + * p 0.01.next_float #=> 0.010000000000000002 + * p 1.0.next_float #=> 1.0000000000000002 + * p 100.0.next_float #=> 100.00000000000001 + * + * p 0.01.next_float - 0.01 #=> 1.734723475976807e-18 + * p 1.0.next_float - 1.0 #=> 2.220446049250313e-16 + * p 100.0.next_float - 100.0 #=> 1.4210854715202004e-14 + * + * f = 0.01; 20.times { printf "%-20a %s\n", f, f.to_s; f = f.next_float } + * #=> 0x1.47ae147ae147bp-7 0.01 + * # 0x1.47ae147ae147cp-7 0.010000000000000002 + * # 0x1.47ae147ae147dp-7 0.010000000000000004 + * # 0x1.47ae147ae147ep-7 0.010000000000000005 + * # 0x1.47ae147ae147fp-7 0.010000000000000007 + * # 0x1.47ae147ae148p-7 0.010000000000000009 + * # 0x1.47ae147ae1481p-7 0.01000000000000001 + * # 0x1.47ae147ae1482p-7 0.010000000000000012 + * # 0x1.47ae147ae1483p-7 0.010000000000000014 + * # 0x1.47ae147ae1484p-7 0.010000000000000016 + * # 0x1.47ae147ae1485p-7 0.010000000000000018 + * # 0x1.47ae147ae1486p-7 0.01000000000000002 + * # 0x1.47ae147ae1487p-7 0.010000000000000021 + * # 0x1.47ae147ae1488p-7 0.010000000000000023 + * # 0x1.47ae147ae1489p-7 0.010000000000000024 + * # 0x1.47ae147ae148ap-7 0.010000000000000026 + * # 0x1.47ae147ae148bp-7 0.010000000000000028 + * # 0x1.47ae147ae148cp-7 0.01000000000000003 + * # 0x1.47ae147ae148dp-7 0.010000000000000031 + * # 0x1.47ae147ae148ep-7 0.010000000000000033 + * + * f = 0.0 + * 100.times { f += 0.1 } + * p f #=> 9.99999999999998 # should be 10.0 in the ideal world. + * p 10-f #=> 1.9539925233402755e-14 # the floating-point error. + * p(10.0.next_float-10) #=> 1.7763568394002505e-15 # 1 ulp (units in the last place). + * p((10-f)/(10.0.next_float-10)) #=> 11.0 # the error is 11 ulp. + * p((10-f)/(10*Float::EPSILON)) #=> 8.8 # approximation of the above. + * p "%a" % f #=> "0x1.3fffffffffff5p+3" # the last hex digit is 5. 16 - 5 = 11 ulp. + * + */ +static VALUE +flo_next_float(VALUE vx) +{ + double x, y; + x = NUM2DBL(vx); + y = nextafter(x, INFINITY); + return DBL2NUM(y); +} + +/* + * call-seq: + * float.prev_float -> float + * + * Returns the previous representable floatint-point number. + * + * (-Float::MAX).prev_float and (-Float::INFINITY).prev_float is -Float::INFINITY. + * + * Float::NAN.prev_float is Float::NAN. + * + * For example: + * + * p 0.01.prev_float #=> 0.009999999999999998 + * p 1.0.prev_float #=> 0.9999999999999999 + * p 100.0.prev_float #=> 99.99999999999999 + * + * p 0.01 - 0.01.prev_float #=> 1.734723475976807e-18 + * p 1.0 - 1.0.prev_float #=> 1.1102230246251565e-16 + * p 100.0 - 100.0.prev_float #=> 1.4210854715202004e-14 + * + * f = 0.01; 20.times { printf "%-20a %s\n", f, f.to_s; f = f.prev_float } + * #=> 0x1.47ae147ae147bp-7 0.01 + * # 0x1.47ae147ae147ap-7 0.009999999999999998 + * # 0x1.47ae147ae1479p-7 0.009999999999999997 + * # 0x1.47ae147ae1478p-7 0.009999999999999995 + * # 0x1.47ae147ae1477p-7 0.009999999999999993 + * # 0x1.47ae147ae1476p-7 0.009999999999999992 + * # 0x1.47ae147ae1475p-7 0.00999999999999999 + * # 0x1.47ae147ae1474p-7 0.009999999999999988 + * # 0x1.47ae147ae1473p-7 0.009999999999999986 + * # 0x1.47ae147ae1472p-7 0.009999999999999985 + * # 0x1.47ae147ae1471p-7 0.009999999999999983 + * # 0x1.47ae147ae147p-7 0.009999999999999981 + * # 0x1.47ae147ae146fp-7 0.00999999999999998 + * # 0x1.47ae147ae146ep-7 0.009999999999999978 + * # 0x1.47ae147ae146dp-7 0.009999999999999976 + * # 0x1.47ae147ae146cp-7 0.009999999999999974 + * # 0x1.47ae147ae146bp-7 0.009999999999999972 + * # 0x1.47ae147ae146ap-7 0.00999999999999997 + * # 0x1.47ae147ae1469p-7 0.009999999999999969 + * # 0x1.47ae147ae1468p-7 0.009999999999999967 + * + */ +static VALUE +flo_prev_float(VALUE vx) +{ + double x, y; + x = NUM2DBL(vx); + y = nextafter(x, -INFINITY); + return DBL2NUM(y); +} + +/* + * call-seq: * float.floor -> integer * * Returns the largest integer less than or equal to +float+. @@ -4106,6 +4219,8 @@ Init_Numeric(void) rb_define_method(rb_cFloat, "nan?", flo_is_nan_p, 0); rb_define_method(rb_cFloat, "infinite?", flo_is_infinite_p, 0); rb_define_method(rb_cFloat, "finite?", flo_is_finite_p, 0); + rb_define_method(rb_cFloat, "next_float", flo_next_float, 0); + rb_define_method(rb_cFloat, "prev_float", flo_prev_float, 0); sym_to = ID2SYM(rb_intern("to")); sym_by = ID2SYM(rb_intern("by")); diff --git a/configure.in b/configure.in index ba40553..83ce1c5 100644 --- a/configure.in +++ b/configure.in @@ -1934,6 +1934,7 @@ AC_REPLACE_FUNCS(isinf) AC_REPLACE_FUNCS(isnan) AC_REPLACE_FUNCS(lgamma_r) AC_REPLACE_FUNCS(memmove) +AC_REPLACE_FUNCS(nextafter) AC_REPLACE_FUNCS(setproctitle) AC_REPLACE_FUNCS(strchr) AC_REPLACE_FUNCS(strerror) diff --git a/include/ruby/missing.h b/include/ruby/missing.h index 5dba4ed..53e0833 100644 --- a/include/ruby/missing.h +++ b/include/ruby/missing.h @@ -168,6 +168,10 @@ RUBY_EXTERN int isnan(double); # endif #endif +#ifndef HAVE_NEXTAFTER +RUBY_EXTERN double nextafter(double x, double y); +#endif + /* #ifndef HAVE_MEMCMP RUBY_EXTERN int memcmp(const void *, const void *, size_t); diff --git a/missing/nextafter.c b/missing/nextafter.c new file mode 100644 index 0000000..4e7a09b --- /dev/null +++ b/missing/nextafter.c @@ -0,0 +1,68 @@ +#include +#include + +/* This function doesn't set errno. It should on POSIX, though. */ + +double +nextafter(double x, double y) +{ + double x1, x2, d; + int e; + + if (isnan(x)) + return x; + if (isnan(y)) + return y; + + if (x == y) + return y; + + if (x == 0) { + /* the minimum "subnormal" float */ + x1 = ldexp(0.5, DBL_MIN_EXP - DBL_MANT_DIG + 1); + if (x1 == 0) + x1 = DBL_MIN; /* the minimum "normal" float */ + if (0 < y) + return x1; + else + return -x1; + } + + if (x < 0) { + if (isinf(x)) + return -DBL_MAX; + if (x == -DBL_MAX && y < 0 && isinf(y)) + return y; + } + else { + if (isinf(x)) + return DBL_MAX; + if (x == DBL_MAX && 0 < y && isinf(y)) + return y; + } + + x1 = frexp(x, &e); + + if (x < y) { + d = DBL_EPSILON/2; + if (x1 == -0.5) { + x1 *= 2; + e--; + } + } + else { + d = -DBL_EPSILON/2; + if (x1 == 0.5) { + x1 *= 2; + e--; + } + } + + if (e < DBL_MIN_EXP) { + d = ldexp(d, DBL_MIN_EXP-e); + } + + x2 = x1 + d; + + return ldexp(x2, e); +} diff --git a/test/ruby/test_float.rb b/test/ruby/test_float.rb index 3ee0a23..952fa7c 100644 --- a/test/ruby/test_float.rb +++ b/test/ruby/test_float.rb @@ -619,4 +619,31 @@ class TestFloat < Test::Unit::TestCase assert_in_epsilon(10.0, ("1."+"1"*300000).to_f*9) end; end + + def test_next_float + smallest = 0.0.next_float + assert_equal(-Float::MAX, (-Float::INFINITY).next_float) + assert_operator(-Float::MAX, :<, (-Float::MAX).next_float) + assert_equal(Float::EPSILON/2, (-1.0).next_float + 1.0) + assert_operator(0.0, :<, smallest) + assert_operator([0.0, smallest], :include?, smallest/2) + assert_equal(Float::EPSILON, 1.0.next_float - 1.0) + assert_equal(Float::INFINITY, Float::MAX.next_float) + assert_equal(Float::INFINITY, Float::INFINITY.next_float) + assert(Float::NAN.next_float.nan?) + end + + def test_prev_float + smallest = 0.0.next_float + assert_equal(-Float::INFINITY, (-Float::INFINITY).prev_float) + assert_equal(-Float::INFINITY, (-Float::MAX).prev_float) + assert_equal(-Float::EPSILON, (-1.0).prev_float + 1.0) + assert_equal(-smallest, 0.0.prev_float) + assert_operator([0.0, 0.0.prev_float], :include?, 0.0.prev_float/2) + assert_equal(-Float::EPSILON/2, 1.0.prev_float - 1.0) + assert_operator(Float::MAX, :>, Float::MAX.prev_float) + assert_equal(Float::MAX, Float::INFINITY.prev_float) + assert(Float::NAN.prev_float.nan?) + end + end diff --git a/ext/-test-/float/depend b/ext/-test-/float/depend new file mode 100644 index 0000000..dff1455 --- /dev/null +++ b/ext/-test-/float/depend @@ -0,0 +1,3 @@ +$(OBJS): $(HDRS) $(ruby_headers) + +nextafter.o: nextafter.c $(top_srcdir)/missing/nextafter.c diff --git a/ext/-test-/float/extconf.rb b/ext/-test-/float/extconf.rb new file mode 100644 index 0000000..0a9a299 --- /dev/null +++ b/ext/-test-/float/extconf.rb @@ -0,0 +1,7 @@ +$INCFLAGS << " -I$(topdir) -I$(top_srcdir)" +$srcs = Dir[File.join($srcdir, "*.{#{SRC_EXT.join(%q{,})}}")] +inits = $srcs.map {|s| File.basename(s, ".*")} +inits.delete("init") +inits.map! {|s|"X(#{s})"} +$defs << "-DTEST_INIT_FUNCS(X)=\"#{inits.join(' ')}\"" +create_makefile("-test-/float") diff --git a/ext/-test-/float/init.c b/ext/-test-/float/init.c new file mode 100644 index 0000000..d962108 --- /dev/null +++ b/ext/-test-/float/init.c @@ -0,0 +1,11 @@ +#include "ruby.h" + +#define init(n) {void Init_##n(VALUE klass); Init_##n(klass);} + +void +Init_float(void) +{ + VALUE mBug = rb_define_module("Bug"); + VALUE klass = rb_define_class_under(mBug, "Float", rb_cObject); + TEST_INIT_FUNCS(init); +} diff --git a/ext/-test-/float/nextafter.c b/ext/-test-/float/nextafter.c new file mode 100644 index 0000000..30fb71f --- /dev/null +++ b/ext/-test-/float/nextafter.c @@ -0,0 +1,36 @@ +#include "ruby.h" + +static VALUE +system_nextafter_m(VALUE klass, VALUE vx, VALUE vy) +{ + double x, y, z; + + x = NUM2DBL(vx); + y = NUM2DBL(vy); + z = nextafter(x, y); + + return DBL2NUM(z); +} + +#define nextafter missing_nextafter +#include "../../../missing/nextafter.c" +#undef nextafter + +static VALUE +missing_nextafter_m(VALUE klass, VALUE vx, VALUE vy) +{ + double x, y, z; + + x = NUM2DBL(vx); + y = NUM2DBL(vy); + z = missing_nextafter(x, y); + + return DBL2NUM(z); +} + +void +Init_nextafter(VALUE klass) +{ + rb_define_singleton_method(klass, "system_nextafter", system_nextafter_m, 2); + rb_define_singleton_method(klass, "missing_nextafter", missing_nextafter_m, 2); +} diff --git a/test/-ext-/float/test_nextafter.rb b/test/-ext-/float/test_nextafter.rb new file mode 100644 index 0000000..b952d7e --- /dev/null +++ b/test/-ext-/float/test_nextafter.rb @@ -0,0 +1,41 @@ +require 'test/unit' +require "-test-/float" + +class TestFloatExt < Test::Unit::TestCase + def test_nextafter + nums = [ + -Float::INFINITY, + -Float::MAX, + -100.0, + -1.0, + -Float::EPSILON, + -Float::MIN/2, + -Math.ldexp(0.5, Float::MIN_EXP - Float::MANT_DIG + 1), + 0.0, + Math.ldexp(0.5, Float::MIN_EXP - Float::MANT_DIG + 1), + Float::MIN/2, + Float::MIN, + Float::EPSILON, + 1.0, + 100.0, + Float::MAX, + Float::INFINITY, + Float::NAN + ] + nums.each {|n1| + nums.each {|n2| + v1 = Bug::Float.missing_nextafter(n1, n2) + v2 = Bug::Float.system_nextafter(n1, n2) + assert_kind_of(Float, v1) + assert_kind_of(Float, v2) + if v1.nan? + assert(v2.nan?, "Bug::Float.system_nextafter(#{n1}, #{n2}).nan?") + else + assert_equal(v1, v2, + "Bug::Float.missing_nextafter(#{'%a' % n1}, #{'%a' % n2}) = #{'%a' % v1} != " + + "#{'%a' % v2} = Bug::Float.system_nextafter(#{'%a' % n1}, #{'%a' % n2})") + end + } + } + end +end