From bed5dea0c0d266ae4d7ae92c05ab56be9faacd0a Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Wed, 11 Jan 2012 23:36:29 +0000
Subject: [PATCH] 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 |  2 ++
 src/real.cpp        | 20 ++++++++++++++++++++
 test/unit/real.cpp  |  8 ++++++++
 3 files changed, 30 insertions(+)

diff --git a/src/lol/math/real.h b/src/lol/math/real.h
index eaba718d..e37bcb26 100644
--- a/src/lol/math/real.h
+++ b/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);
diff --git a/src/real.cpp b/src/real.cpp
index c2dfcaba..a96cb112 100644
--- a/src/real.cpp
+++ b/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;
diff --git a/test/unit/real.cpp b/test/unit/real.cpp
index b3536646..664d68d6 100644
--- a/test/unit/real.cpp
+++ b/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;