jni/ruby/missing/nextafter.c


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

#include "ruby/missing.h"

#include <math.h>
#include <float.h>

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

    if (x2 == 0.0) {
        if (x1 < 0)
            return -0.0;
        else
            return +0.0;
    }

    return ldexp(x2, e);
}