From 221f469c3fdc619b81379b064fdc226d9de5ac3c Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Sun, 2 Oct 2011 23:36:25 +0000
Subject: [PATCH] core: implement asin() and acos() for real numbers and add
 unit tests for these functions.

---
 src/real.cpp       | 54 ++++++++++++++++++++++++++++++++++++++++++++++
 src/real.h         |  3 +++
 test/unit/real.cpp | 20 +++++++++++++++++
 3 files changed, 77 insertions(+)

diff --git a/src/real.cpp b/src/real.cpp
index e6178885..c2b93497 100644
--- a/src/real.cpp
+++ b/src/real.cpp
@@ -96,6 +96,11 @@ real::operator double() const
     return u.d;
 }
 
+real real::operator +() const
+{
+    return *this;
+}
+
 real real::operator -() const
 {
     real ret = *this;
@@ -668,6 +673,55 @@ real cos(real const &x)
     return ret;
 }
 
+static real asinacos(real const &x, bool is_asin, bool is_negative)
+{
+    /* Strategy for asin(): in [-0.5..0.5], use a Taylor series around
+     * zero. In [0.5..1], use asin(x) = π/2 - 2*asin(sqrt((1-x)/2)), and
+     * in [-1..-0.5] just revert the sign.
+     * Strategy for acos(): use acos(x) = π/2 - asin(x) and try not to
+     * lose the precision around x=1. */
+    real absx = fabs(x);
+    bool around_zero = (absx < (real::R_1 >> 1));
+
+    if (!around_zero)
+        absx = sqrt((real::R_1 - absx) >> 1);
+
+    real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1;
+    for (int i = 1; i < 280; i++)
+    {
+        xn *= x2;
+        ret += (fact1 * xn / ((real)(2 * i + 1) * fact2)) >> (i * 2);
+        fact1 *= (real)((2 * i + 1) * (2 * i + 2));
+        fact2 *= (real)((i + 1) * (i + 1));
+    }
+
+    if (is_negative)
+        ret = -ret;
+
+    if (around_zero)
+        ret = is_asin ? ret : real::R_PI_2 - ret;
+    else
+    {
+        real adjust = is_negative ? real::R_PI : real::R_0;
+        if (is_asin)
+            ret = real::R_PI_2 - adjust - (ret << 1);
+        else
+            ret = adjust + (ret << 1);
+    }
+
+    return ret;
+}
+
+real asin(real const &x)
+{
+    return asinacos(x, true, x.m_signexp >> 31);
+}
+
+real acos(real const &x)
+{
+    return asinacos(x, false, x.m_signexp >> 31);
+}
+
 real atan(real const &x)
 {
     /* Computing atan(x): we choose a different Taylor series depending on
diff --git a/src/real.h b/src/real.h
index 5d3838cb..0b2b8a05 100644
--- a/src/real.h
+++ b/src/real.h
@@ -36,6 +36,7 @@ public:
     operator int() const;
     operator unsigned int() const;
 
+    real operator +() const;
     real operator -() const;
     real operator +(real const &x) const;
     real operator -(real const &x) const;
@@ -70,6 +71,8 @@ public:
 
     friend real sin(real const &x);
     friend real cos(real const &x);
+    friend real asin(real const &x);
+    friend real acos(real const &x);
     friend real atan(real const &x);
 
     void print(int ndigits = 150) const;
diff --git a/test/unit/real.cpp b/test/unit/real.cpp
index 8753e7fc..ace5b4b9 100644
--- a/test/unit/real.cpp
+++ b/test/unit/real.cpp
@@ -236,6 +236,26 @@ LOLUNIT_FIXTURE(RealTest)
         LOLUNIT_ASSERT(a);
         LOLUNIT_ASSERT(!!a);
     }
+
+    LOLUNIT_TEST(AsinAcos)
+    {
+        double tests[14] =
+        {
+            -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++)
+        {
+            double a = tests[n] / 1024;
+            double b = sin(asin((real)a));
+            double c = cos(acos((real)a));
+
+            LOLUNIT_SET_CONTEXT(a);
+            LOLUNIT_ASSERT_DOUBLES_EQUAL(b, a, 1e-100);
+            LOLUNIT_ASSERT_DOUBLES_EQUAL(c, a, 1e-100);
+        }
+    }
 };
 
 } /* namespace lol */