Project

General

Profile

Feature #9834 ยป next_float-and-prev_float.patch

akr (Akira Tanaka), 05/13/2014 09:47 AM

View differences:

numeric.c
1480 1480

  
1481 1481
/*
1482 1482
 *  call-seq:
1483
 *     float.next_float  ->  float
1484
 *
1485
 *  Returns the next representable floating-point number.
1486
 *
1487
 *  Float::MAX.next_float and Float::INFINITY.next_float is Float::INFINITY.
1488
 *
1489
 *  Float::NAN.next_float is Float::NAN.
1490
 *
1491
 *  For example:
1492
 *
1493
 *    p 0.01.next_float  #=> 0.010000000000000002
1494
 *    p 1.0.next_float   #=> 1.0000000000000002
1495
 *    p 100.0.next_float #=> 100.00000000000001
1496
 *
1497
 *    p 0.01.next_float - 0.01   #=> 1.734723475976807e-18
1498
 *    p 1.0.next_float - 1.0     #=> 2.220446049250313e-16
1499
 *    p 100.0.next_float - 100.0 #=> 1.4210854715202004e-14
1500
 *
1501
 *    f = 0.01; 20.times { printf "%-20a %s\n", f, f.to_s; f = f.next_float }
1502
 *    #=> 0x1.47ae147ae147bp-7 0.01
1503
 *    #   0x1.47ae147ae147cp-7 0.010000000000000002
1504
 *    #   0x1.47ae147ae147dp-7 0.010000000000000004
1505
 *    #   0x1.47ae147ae147ep-7 0.010000000000000005
1506
 *    #   0x1.47ae147ae147fp-7 0.010000000000000007
1507
 *    #   0x1.47ae147ae148p-7  0.010000000000000009
1508
 *    #   0x1.47ae147ae1481p-7 0.01000000000000001
1509
 *    #   0x1.47ae147ae1482p-7 0.010000000000000012
1510
 *    #   0x1.47ae147ae1483p-7 0.010000000000000014
1511
 *    #   0x1.47ae147ae1484p-7 0.010000000000000016
1512
 *    #   0x1.47ae147ae1485p-7 0.010000000000000018
1513
 *    #   0x1.47ae147ae1486p-7 0.01000000000000002
1514
 *    #   0x1.47ae147ae1487p-7 0.010000000000000021
1515
 *    #   0x1.47ae147ae1488p-7 0.010000000000000023
1516
 *    #   0x1.47ae147ae1489p-7 0.010000000000000024
1517
 *    #   0x1.47ae147ae148ap-7 0.010000000000000026
1518
 *    #   0x1.47ae147ae148bp-7 0.010000000000000028
1519
 *    #   0x1.47ae147ae148cp-7 0.01000000000000003
1520
 *    #   0x1.47ae147ae148dp-7 0.010000000000000031
1521
 *    #   0x1.47ae147ae148ep-7 0.010000000000000033
1522
 *
1523
 *    f = 0.0
1524
 *    100.times { f += 0.1 }
1525
 *    p f                            #=> 9.99999999999998       # should be 10.0 in the ideal world.
1526
 *    p 10-f                         #=> 1.9539925233402755e-14 # the floating-point error.
1527
 *    p(10.0.next_float-10)          #=> 1.7763568394002505e-15 # 1 ulp (units in the last place).
1528
 *    p((10-f)/(10.0.next_float-10)) #=> 11.0                   # the error is 11 ulp.
1529
 *    p((10-f)/(10*Float::EPSILON))  #=> 8.8                    # approximation of the above.
1530
 *    p "%a" % f                     #=> "0x1.3fffffffffff5p+3" # the last hex digit is 5.  16 - 5 = 11 ulp.
1531
 *
1532
 */
1533
static VALUE
1534
flo_next_float(VALUE vx)
1535
{
1536
    double x, y;
1537
    x = NUM2DBL(vx);
1538
    y = nextafter(x, INFINITY);
1539
    return DBL2NUM(y);
1540
}
1541

  
1542
/*
1543
 *  call-seq:
1544
 *     float.prev_float  ->  float
1545
 *
1546
 *  Returns the previous representable floatint-point number.
1547
 *
1548
 *  (-Float::MAX).prev_float and (-Float::INFINITY).prev_float is -Float::INFINITY.
1549
 *
1550
 *  Float::NAN.prev_float is Float::NAN.
1551
 *
1552
 *  For example:
1553
 *
1554
 *    p 0.01.prev_float  #=> 0.009999999999999998
1555
 *    p 1.0.prev_float   #=> 0.9999999999999999
1556
 *    p 100.0.prev_float #=> 99.99999999999999
1557
 *
1558
 *    p 0.01 - 0.01.prev_float   #=> 1.734723475976807e-18
1559
 *    p 1.0 - 1.0.prev_float     #=> 1.1102230246251565e-16
1560
 *    p 100.0 - 100.0.prev_float #=> 1.4210854715202004e-14
1561
 *
1562
 *    f = 0.01; 20.times { printf "%-20a %s\n", f, f.to_s; f = f.prev_float }
1563
 *    #=> 0x1.47ae147ae147bp-7 0.01
1564
 *    #   0x1.47ae147ae147ap-7 0.009999999999999998
1565
 *    #   0x1.47ae147ae1479p-7 0.009999999999999997
1566
 *    #   0x1.47ae147ae1478p-7 0.009999999999999995
1567
 *    #   0x1.47ae147ae1477p-7 0.009999999999999993
1568
 *    #   0x1.47ae147ae1476p-7 0.009999999999999992
1569
 *    #   0x1.47ae147ae1475p-7 0.00999999999999999
1570
 *    #   0x1.47ae147ae1474p-7 0.009999999999999988
1571
 *    #   0x1.47ae147ae1473p-7 0.009999999999999986
1572
 *    #   0x1.47ae147ae1472p-7 0.009999999999999985
1573
 *    #   0x1.47ae147ae1471p-7 0.009999999999999983
1574
 *    #   0x1.47ae147ae147p-7  0.009999999999999981
1575
 *    #   0x1.47ae147ae146fp-7 0.00999999999999998
1576
 *    #   0x1.47ae147ae146ep-7 0.009999999999999978
1577
 *    #   0x1.47ae147ae146dp-7 0.009999999999999976
1578
 *    #   0x1.47ae147ae146cp-7 0.009999999999999974
1579
 *    #   0x1.47ae147ae146bp-7 0.009999999999999972
1580
 *    #   0x1.47ae147ae146ap-7 0.00999999999999997
1581
 *    #   0x1.47ae147ae1469p-7 0.009999999999999969
1582
 *    #   0x1.47ae147ae1468p-7 0.009999999999999967
1583
 *
1584
 */
1585
static VALUE
1586
flo_prev_float(VALUE vx)
1587
{
1588
    double x, y;
1589
    x = NUM2DBL(vx);
1590
    y = nextafter(x, -INFINITY);
1591
    return DBL2NUM(y);
1592
}
1593

  
1594
/*
1595
 *  call-seq:
1483 1596
 *     float.floor  ->  integer
1484 1597
 *
1485 1598
 *  Returns the largest integer less than or equal to +float+.
......
4106 4219
    rb_define_method(rb_cFloat, "nan?",      flo_is_nan_p, 0);
4107 4220
    rb_define_method(rb_cFloat, "infinite?", flo_is_infinite_p, 0);
4108 4221
    rb_define_method(rb_cFloat, "finite?",   flo_is_finite_p, 0);
4222
    rb_define_method(rb_cFloat, "next_float", flo_next_float, 0);
4223
    rb_define_method(rb_cFloat, "prev_float", flo_prev_float, 0);
4109 4224

  
4110 4225
    sym_to = ID2SYM(rb_intern("to"));
4111 4226
    sym_by = ID2SYM(rb_intern("by"));
configure.in
1934 1934
AC_REPLACE_FUNCS(isnan)
1935 1935
AC_REPLACE_FUNCS(lgamma_r)
1936 1936
AC_REPLACE_FUNCS(memmove)
1937
AC_REPLACE_FUNCS(nextafter)
1937 1938
AC_REPLACE_FUNCS(setproctitle)
1938 1939
AC_REPLACE_FUNCS(strchr)
1939 1940
AC_REPLACE_FUNCS(strerror)
include/ruby/missing.h
168 168
# endif
169 169
#endif
170 170

  
171
#ifndef HAVE_NEXTAFTER
172
RUBY_EXTERN double nextafter(double x, double y);
173
#endif
174

  
171 175
/*
172 176
#ifndef HAVE_MEMCMP
173 177
RUBY_EXTERN int memcmp(const void *, const void *, size_t);
missing/nextafter.c
1
#include <math.h>
2
#include <float.h>
3

  
4
/* This function doesn't set errno.  It should on POSIX, though. */
5

  
6
double
7
nextafter(double x, double y)
8
{
9
    double x1, x2, d;
10
    int e;
11

  
12
    if (isnan(x))
13
        return x;
14
    if (isnan(y))
15
        return y;
16

  
17
    if (x == y)
18
        return y;
19

  
20
    if (x == 0) {
21
        /* the minimum "subnormal" float */
22
        x1 = ldexp(0.5, DBL_MIN_EXP - DBL_MANT_DIG + 1);
23
        if (x1 == 0)
24
            x1 = DBL_MIN; /* the minimum "normal" float */
25
        if (0 < y)
26
            return x1;
27
        else
28
            return -x1;
29
    }
30

  
31
    if (x < 0) {
32
        if (isinf(x))
33
            return -DBL_MAX;
34
        if (x == -DBL_MAX && y < 0 && isinf(y))
35
            return y;
36
    }
37
    else {
38
        if (isinf(x))
39
            return DBL_MAX;
40
        if (x == DBL_MAX && 0 < y && isinf(y))
41
            return y;
42
    }
43

  
44
    x1 = frexp(x, &e);
45

  
46
    if (x < y) {
47
        d = DBL_EPSILON/2;
48
        if (x1 == -0.5) {
49
            x1 *= 2;
50
            e--;
51
        }
52
    }
53
    else {
54
        d = -DBL_EPSILON/2;
55
        if (x1 == 0.5) {
56
            x1 *= 2;
57
            e--;
58
        }
59
    }
60

  
61
    if (e < DBL_MIN_EXP) {
62
        d = ldexp(d, DBL_MIN_EXP-e);
63
    }
64

  
65
    x2 = x1 + d;
66

  
67
    return ldexp(x2, e);
68
}
test/ruby/test_float.rb
619 619
    assert_in_epsilon(10.0, ("1."+"1"*300000).to_f*9)
620 620
    end;
621 621
  end
622

  
623
  def test_next_float
624
    smallest = 0.0.next_float
625
    assert_equal(-Float::MAX, (-Float::INFINITY).next_float)
626
    assert_operator(-Float::MAX, :<, (-Float::MAX).next_float)
627
    assert_equal(Float::EPSILON/2, (-1.0).next_float + 1.0)
628
    assert_operator(0.0, :<, smallest)
629
    assert_operator([0.0, smallest], :include?, smallest/2)
630
    assert_equal(Float::EPSILON, 1.0.next_float - 1.0)
631
    assert_equal(Float::INFINITY, Float::MAX.next_float)
632
    assert_equal(Float::INFINITY, Float::INFINITY.next_float)
633
    assert(Float::NAN.next_float.nan?)
634
  end
635

  
636
  def test_prev_float
637
    smallest = 0.0.next_float
638
    assert_equal(-Float::INFINITY, (-Float::INFINITY).prev_float)
639
    assert_equal(-Float::INFINITY, (-Float::MAX).prev_float)
640
    assert_equal(-Float::EPSILON, (-1.0).prev_float + 1.0)
641
    assert_equal(-smallest, 0.0.prev_float)
642
    assert_operator([0.0, 0.0.prev_float], :include?, 0.0.prev_float/2)
643
    assert_equal(-Float::EPSILON/2, 1.0.prev_float - 1.0)
644
    assert_operator(Float::MAX, :>, Float::MAX.prev_float)
645
    assert_equal(Float::MAX, Float::INFINITY.prev_float)
646
    assert(Float::NAN.prev_float.nan?)
647
  end
648

  
622 649
end
ext/-test-/float/depend
1
$(OBJS): $(HDRS) $(ruby_headers)
2

  
3
nextafter.o: nextafter.c $(top_srcdir)/missing/nextafter.c
ext/-test-/float/extconf.rb
1
$INCFLAGS << " -I$(topdir) -I$(top_srcdir)"
2
$srcs = Dir[File.join($srcdir, "*.{#{SRC_EXT.join(%q{,})}}")]
3
inits = $srcs.map {|s| File.basename(s, ".*")}
4
inits.delete("init")
5
inits.map! {|s|"X(#{s})"}
6
$defs << "-DTEST_INIT_FUNCS(X)=\"#{inits.join(' ')}\""
7
create_makefile("-test-/float")
ext/-test-/float/init.c
1
#include "ruby.h"
2

  
3
#define init(n) {void Init_##n(VALUE klass); Init_##n(klass);}
4

  
5
void
6
Init_float(void)
7
{
8
    VALUE mBug = rb_define_module("Bug");
9
    VALUE klass = rb_define_class_under(mBug, "Float", rb_cObject);
10
    TEST_INIT_FUNCS(init);
11
}
ext/-test-/float/nextafter.c
1
#include "ruby.h"
2

  
3
static VALUE
4
system_nextafter_m(VALUE klass, VALUE vx, VALUE vy)
5
{
6
    double x, y, z;
7

  
8
    x = NUM2DBL(vx);
9
    y = NUM2DBL(vy);
10
    z = nextafter(x, y);
11

  
12
    return DBL2NUM(z);
13
}
14

  
15
#define nextafter missing_nextafter
16
#include "../../../missing/nextafter.c"
17
#undef nextafter
18

  
19
static VALUE
20
missing_nextafter_m(VALUE klass, VALUE vx, VALUE vy)
21
{
22
    double x, y, z;
23

  
24
    x = NUM2DBL(vx);
25
    y = NUM2DBL(vy);
26
    z = missing_nextafter(x, y);
27

  
28
    return DBL2NUM(z);
29
}
30

  
31
void
32
Init_nextafter(VALUE klass)
33
{
34
    rb_define_singleton_method(klass, "system_nextafter", system_nextafter_m, 2);
35
    rb_define_singleton_method(klass, "missing_nextafter", missing_nextafter_m, 2);
36
}
test/-ext-/float/test_nextafter.rb
1
require 'test/unit'
2
require "-test-/float"
3

  
4
class TestFloatExt < Test::Unit::TestCase
5
  def test_nextafter
6
    nums = [
7
      -Float::INFINITY,
8
      -Float::MAX,
9
      -100.0,
10
      -1.0,
11
      -Float::EPSILON,
12
      -Float::MIN/2,
13
      -Math.ldexp(0.5, Float::MIN_EXP - Float::MANT_DIG + 1),
14
      0.0,
15
      Math.ldexp(0.5, Float::MIN_EXP - Float::MANT_DIG + 1),
16
      Float::MIN/2,
17
      Float::MIN,
18
      Float::EPSILON,
19
      1.0,
20
      100.0,
21
      Float::MAX,
22
      Float::INFINITY,
23
      Float::NAN
24
    ]
25
    nums.each {|n1|
26
      nums.each {|n2|
27
        v1 = Bug::Float.missing_nextafter(n1, n2)
28
        v2 = Bug::Float.system_nextafter(n1, n2)
29
        assert_kind_of(Float, v1)
30
        assert_kind_of(Float, v2)
31
        if v1.nan?
32
          assert(v2.nan?, "Bug::Float.system_nextafter(#{n1}, #{n2}).nan?")
33
        else
34
          assert_equal(v1, v2,
35
            "Bug::Float.missing_nextafter(#{'%a' % n1}, #{'%a' % n2}) = #{'%a' % v1} != " +
36
            "#{'%a' % v2} = Bug::Float.system_nextafter(#{'%a' % n1}, #{'%a' % n2})")
37
        end
38
      }
39
    }
40
  end
41
end