math: implement ulp() for reals, which returns the smallest real y > 0 such

that x + y != x, and nextafter() which behaves like the C function.

src/lol/math/real.h

@@ -88,6 +88,8 @@ public:
     friend real frexp(real const &x, int *exp);
     friend real ldexp(real const &x, int exp);
     friend real modf(real const &x, real *iptr);
+    friend real ulp(real const &x);
+    friend real nextafter(real const &x, real const &y);
 
     /* Power functions */
     friend real re(real const &x);

src/real.cpp

@@ -938,6 +938,26 @@ real modf(real const &x, real *iptr)
     return copysign(absx - tmp, x);
 }
 
+real ulp(real const &x)
+{
+    real ret = real::R_1;
+    if (x)
+        ret.m_signexp = x.m_signexp + 1 - real::BIGITS * real::BIGIT_BITS;
+    else
+        ret.m_signexp = 0;
+    return ret;
+}
+
+real nextafter(real const &x, real const &y)
+{
+    if (x == y)
+        return x;
+    else if (x < y)
+        return x + ulp(x);
+    else
+        return x - ulp(x);
+}
+
 real copysign(real const &x, real const &y)
 {
     real ret = x;

test/unit/real.cpp

@@ -242,6 +242,14 @@ LOLUNIT_FIXTURE(RealTest)
         LOLUNIT_ASSERT_EQUAL((double)ldexp(a3, -7), 0.0);
     }
 
+    LOLUNIT_TEST(Ulp)
+    {
+        real a1 = real::R_PI;
+
+        LOLUNIT_ASSERT_NOT_EQUAL((double)(a1 + ulp(a1) - a1), 0.0);
+        LOLUNIT_ASSERT_EQUAL((double)(a1 + ulp(a1) / 2 - a1), 0.0);
+    }
+
     LOLUNIT_TEST(Bool)
     {
         real a = 0.0;