From 064856b443fff5217c895cc3c43e688e942e33db Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Sun, 2 Oct 2011 23:36:28 +0000
Subject: [PATCH] core: implement floor() and ceil() for real numbers.

---
 src/real.cpp       | 46 ++++++++++++++++++++++++++++++++++++++++++++++
 src/real.h         |  3 +++
 test/unit/real.cpp | 43 +++++++++++++++++++++++++++++++++++++++++--
 3 files changed, 90 insertions(+), 2 deletions(-)

diff --git a/src/real.cpp b/src/real.cpp
index c2b93497..ff8579e5 100644
--- a/src/real.cpp
+++ b/src/real.cpp
@@ -639,6 +639,52 @@ real exp(real const &x)
     return x1;
 }
 
+real floor(real const &x)
+{
+    /* Strategy for floor(x):
+     *  - if negative, return -ceil(-x)
+     *  - if zero or negative zero, return x
+     *  - if less than one, return zero
+     *  - otherwise, if e is the exponent, clear all bits except the
+     *    first e. */
+    if (x < -real::R_0)
+        return -ceil(-x);
+    if (!x)
+        return x;
+    if (x < real::R_1)
+        return real::R_0;
+
+    real ret = x;
+    int exponent = x.m_signexp - (1 << 30) + 1;
+
+    for (int i = 0; i < real::BIGITS; i++)
+    {
+        if (exponent <= 0)
+            ret.m_mantissa[i] = 0;
+        else if (exponent < 8 * (int)sizeof(uint16_t))
+            ret.m_mantissa[i] &= ~((1 << (16 - exponent)) - 1);
+
+        exponent -= 8 * sizeof(uint16_t);
+    }
+
+    return ret;
+}
+
+real ceil(real const &x)
+{
+    /* Strategy for ceil(x):
+     *  - if negative, return -floor(-x)
+     *  - if x == floor(x), return x
+     *  - otherwise, return floor(x) + 1 */
+    if (x < -real::R_0)
+        return -floor(-x);
+    real ret = floor(x);
+    if (x == ret)
+        return ret;
+    else
+        return ret + real::R_1;
+}
+
 real sin(real const &x)
 {
     real ret = 0.0, fact = 1.0, xn = x, x2 = x * x;
diff --git a/src/real.h b/src/real.h
index 0b2b8a05..e88707a5 100644
--- a/src/real.h
+++ b/src/real.h
@@ -69,6 +69,9 @@ public:
     friend real log(real const &x);
     friend real exp(real const &x);
 
+    friend real floor(real const &x);
+    friend real ceil(real const &x);
+
     friend real sin(real const &x);
     friend real cos(real const &x);
     friend real asin(real const &x);
diff --git a/test/unit/real.cpp b/test/unit/real.cpp
index ace5b4b9..cf482789 100644
--- a/test/unit/real.cpp
+++ b/test/unit/real.cpp
@@ -239,13 +239,13 @@ LOLUNIT_FIXTURE(RealTest)
 
     LOLUNIT_TEST(AsinAcos)
     {
-        double tests[14] =
+        double tests[] =
         {
             -1024.0, -1023.0, -513.0, -512.0, -511.0, -1.0, -0.0,
             0.0, 1.0, 511.0, 512.0, 513.0, 1023.0, 1024.0
         };
 
-        for (int n = 0; n < 14; n++)
+        for (unsigned int n = 0; n < sizeof(tests) / sizeof(*tests); n++)
         {
             double a = tests[n] / 1024;
             double b = sin(asin((real)a));
@@ -256,6 +256,45 @@ LOLUNIT_FIXTURE(RealTest)
             LOLUNIT_ASSERT_DOUBLES_EQUAL(c, a, 1e-100);
         }
     }
+
+    LOLUNIT_TEST(FloorCeil)
+    {
+        double tests[] =
+        {
+            -2.0,  -2.0, -2.0,
+            -1.5,  -2.0, -1.0,
+            -1.0,  -1.0, -1.0,
+            -0.0,  -0.0, -0.0,
+             0.0,   0.0,  0.0,
+             0.25,  0.0,  1.0,
+             0.375, 0.0,  1.0,
+             0.5,   0.0,  1.0,
+             1.0,   1.0,  1.0,
+             1.5,   1.0,  2.0,
+             2.0,   2.0,  2.0,
+             2.5,   2.0,  3.0,
+             3.0,   3.0,  3.0,
+            8192.0,     8192.0, 8192.0,
+            8192.03125, 8192.0, 8193.0,
+            8192.5,     8192.0, 8193.0,
+            8193.0,     8193.0, 8193.0,
+            549755813888.0,     549755813888.0, 549755813888.0,
+            549755813888.03125, 549755813888.0, 549755813889.0,
+            549755813888.5,     549755813888.0, 549755813889.0,
+            549755813889.0,     549755813889.0, 549755813889.0,
+        };
+
+        for (unsigned int n = 0; n < sizeof(tests) / sizeof(*tests); n += 3)
+        {
+            double a0 = floor((real)tests[n]);
+            double b0 = tests[n + 1];
+            double a1 = ceil((real)tests[n]);
+            double b1 = tests[n + 2];
+
+            LOLUNIT_ASSERT_EQUAL(b0, a0);
+            LOLUNIT_ASSERT_EQUAL(b1, a1);
+        }
+    }
 };
 
 } /* namespace lol */