From e96ef7af323009ef3690c650ba6aae168e42d277 Mon Sep 17 00:00:00 2001
From: Danny Robson <danny@nerdcruft.net>
Date: Mon, 17 Oct 2016 22:14:57 +1100
Subject: [PATCH] quaternion: move out of coord infrastructure

---
 Makefile.am         |   1 +
 coord/fwd.hpp       |   1 -
 coord/ops.hpp       |   3 --
 quaternion.cpp      | 127 +++++++++++++++++++++-----------------------
 quaternion.hpp      | 110 +++++++++++++++++++++++++++-----------
 quaternion.ipp      |  99 ++++++++++++++++++++++++++++++++++
 test/quaternion.cpp |  20 +++----
 7 files changed, 251 insertions(+), 110 deletions(-)
 create mode 100644 quaternion.ipp

diff --git a/Makefile.am b/Makefile.am
index a9157056..a443ef24 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -206,6 +206,7 @@ UTIL_FILES   =			\
     preprocessor.hpp		\
     quaternion.cpp		\
     quaternion.hpp		\
+    quaternion.ipp		\
     raii.hpp			\
     rand/lcg.cpp		\
     rand/lcg.hpp		\
diff --git a/coord/fwd.hpp b/coord/fwd.hpp
index 3d144148..e1c3347e 100644
--- a/coord/fwd.hpp
+++ b/coord/fwd.hpp
@@ -21,7 +21,6 @@ namespace util {
     template <size_t,typename> struct colour;
     template <size_t,typename> struct extent;
     template <size_t,typename> struct point;
-    template <size_t,typename> struct quaternion;
     template <size_t,typename> struct vector;
 }
 
diff --git a/coord/ops.hpp b/coord/ops.hpp
index 5ce9b01c..2484ba92 100644
--- a/coord/ops.hpp
+++ b/coord/ops.hpp
@@ -58,7 +58,6 @@ namespace util {
         struct has_norm : public std::false_type { };
 
         template <> struct has_norm<vector>     : public std::true_type { };
-        template <> struct has_norm<quaternion> : public std::true_type { };
 
         template <template <size_t,typename> class K>
         constexpr auto has_norm_v = has_norm<K>::value;
@@ -71,7 +70,6 @@ namespace util {
         template <> struct has_scalar_op<colour>     : public std::true_type { };
         template <> struct has_scalar_op<extent>     : public std::true_type { };
         template <> struct has_scalar_op<point>      : public std::true_type { };
-        template <> struct has_scalar_op<quaternion> : public std::true_type { };
         template <> struct has_scalar_op<vector>     : public std::true_type { };
 
         template <template <size_t,typename> class K>
@@ -84,7 +82,6 @@ namespace util {
     template <size_t S, typename T> struct is_coord<extent<S,T>>     : std::true_type { };
     template <size_t S, typename T> struct is_coord<vector<S,T>>     : std::true_type { };
     template <size_t S, typename T> struct is_coord<colour<S,T>>     : std::true_type { };
-    template <size_t S, typename T> struct is_coord<quaternion<S,T>> : std::true_type { };
 
     template <class K>
     constexpr bool
diff --git a/quaternion.cpp b/quaternion.cpp
index a5268073..7f956afe 100644
--- a/quaternion.cpp
+++ b/quaternion.cpp
@@ -11,7 +11,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  *
- * Copyright 2011 Danny Robson <danny@nerdcruft.net>
+ * Copyright 2011-2016 Danny Robson <danny@nerdcruft.net>
  */
 
 
@@ -22,19 +22,15 @@
 
 #include <cmath>
 
-//-----------------------------------------------------------------------------
+
+///////////////////////////////////////////////////////////////////////////////
 using util::quaternion;
 
 
-//-----------------------------------------------------------------------------
-template<> const quaternion<4,  float> quaternion<4,  float>::IDENTITY = { 1, 0, 0, 0 };
-template<> const quaternion<4, double> quaternion<4, double>::IDENTITY = { 1, 0, 0, 0 };
-
-
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-quaternion<S,T>
-quaternion<S,T>::angle_axis (const T radians, const vector<3,T> axis)
+template <typename T>
+quaternion<T>
+quaternion<T>::angle_axis (const T radians, const vector<3,T> axis)
 {
     CHECK (is_normalised (axis));
 
@@ -48,9 +44,9 @@ quaternion<S,T>::angle_axis (const T radians, const vector<3,T> axis)
 
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-quaternion<S,T>
-quaternion<S,T>::from_euler (vector<3,T> angles)
+template <typename T>
+quaternion<T>
+quaternion<T>::from_euler (vector<3,T> angles)
 {
     auto half = angles / 2;
 
@@ -69,9 +65,9 @@ quaternion<S,T>::from_euler (vector<3,T> angles)
 ///////////////////////////////////////////////////////////////////////////////
 // vector-to-vector rotation algorithm from:
 //     http://lolengine.net/blog/2014/02/24/quaternion-from-two-vectors-final
-template <size_t S, typename T>
-quaternion<S,T>
-quaternion<S,T>::from_to (const vector<3,T> u, const vector<3,T> v)
+template <typename T>
+quaternion<T>
+quaternion<T>::from_to (const vector<3,T> u, const vector<3,T> v)
 {
     CHECK (is_normalised (u));
     CHECK (is_normalised (v));
@@ -103,28 +99,24 @@ quaternion<S,T>::from_to (const vector<3,T> u, const vector<3,T> v)
         w = cross(u, v);
     }
 
-    return normalised (util::quaternion<4,T> (real_part, w.x, w.y, w.z));
+    return normalised (util::quaternion<T> {real_part, w.x, w.y, w.z});
 #endif
 }
 
 
 ///////////////////////////////////////////////////////////////////////////////
 template <typename T>
-quaternion<4,T>
-util::conjugate (quaternion<4,T> q)
+quaternion<T>
+util::conjugate (quaternion<T> q)
 {
     return { q.w, -q.x, -q.y, -q.z };
 }
 
 
-//-----------------------------------------------------------------------------
-template quaternion<4,float> util::conjugate (quaternion<4,float>);
-
-
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
-quaternion<S,T>
-util::operator* (const quaternion<S,T> a, const quaternion<S,T> b)
+template <typename T>
+quaternion<T>
+util::operator* (const quaternion<T> a, const quaternion<T> b)
 {
     return {
         a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
@@ -134,18 +126,16 @@ util::operator* (const quaternion<S,T> a, const quaternion<S,T> b)
     };
 }
 
-template quaternion<4,float> util::operator* (quaternion<4,float>, quaternion<4,float>);
-
 
 //-----------------------------------------------------------------------------
-template <size_t S, typename T>
-quaternion<S,T>
-util::operator/ (const quaternion<S,T> a, const quaternion<S,T> b)
+template <typename T>
+quaternion<T>
+util::operator/ (const quaternion<T> a, const quaternion<T> b)
 {
     CHECK (is_normalised (a));
     CHECK (is_normalised (b));
 
-    return quaternion<S,T> {
+    return quaternion<T> {
           a.w * b.w + a.x * b.x + a.y * b.y + a.z * b.z,
         - a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
         - a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
@@ -153,13 +143,11 @@ util::operator/ (const quaternion<S,T> a, const quaternion<S,T> b)
     };
 }
 
-template quaternion<4,float> util::operator/ (quaternion<4,float>, quaternion<4,float>);
-
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
+template <typename T>
 util::matrix4<T>
-quaternion<S, T>::as_matrix (void) const
+quaternion<T>::as_matrix (void) const
 {
     CHECK (is_normalised (*this));
 
@@ -180,13 +168,13 @@ quaternion<S, T>::as_matrix (void) const
 // https://gamedev.stackexchange.com/questions/28395/rotating-vector3-by-a-quaternion
 template <typename T>
 util::vector3<T>
-util::rotate (vector3<T> v, quaternion<4,T> q)
+util::rotate (vector3<T> v, quaternion<T> q)
 {
     CHECK (is_normalised (v));
 
 #if 0
     // Naive:
-    quaternion<4,T> p { 0, v.x, v.y, v.z };
+    quaternion<T> p { 0, v.x, v.y, v.z };
     auto p_ = q * p * conjugate (q);
     return { p_.x, p_.y, p_.z };
 #else
@@ -196,17 +184,12 @@ util::rotate (vector3<T> v, quaternion<4,T> q)
 }
 
 
-//-----------------------------------------------------------------------------
-template util::vector3f util::rotate (util::vector3f, util::quaternionf);
-template util::vector3d util::rotate (util::vector3d, util::quaterniond);
-
-
 ///////////////////////////////////////////////////////////////////////////////
 // based on the implementation at:
 // http://www.opengl-tutorial.org/intermediate-tutorials/tutorial-17-quaternions/
-template <size_t S, typename T>
-quaternion<S,T>
-quaternion<S,T>::look (vector<3,T> fwd, vector<3,T> up)
+template <typename T>
+quaternion<T>
+quaternion<T>::look (vector<3,T> fwd, vector<3,T> up)
 {
     CHECK (is_normalised (fwd));
     CHECK (is_normalised (up));
@@ -232,36 +215,50 @@ quaternion<S,T>::look (vector<3,T> fwd, vector<3,T> up)
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template <size_t S, typename T>
+template <typename T>
+bool
+util::almost_equal (quaternion<T> a, quaternion<T> b)
+{
+    return almost_equal (a.w, b.w) &&
+           almost_equal (a.x, b.x) &&
+           almost_equal (a.y, b.y) &&
+           almost_equal (a.z, b.z);
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+template <typename T>
 std::ostream&
-util::operator<< (std::ostream &os, const quaternion<S,T> q)
+util::operator<< (std::ostream &os, const quaternion<T> q)
 {
     return os << "[" << q.w << ", " << q.x << ", " << q.y << ", " << q.z << "]";
 }
 
 
-//-----------------------------------------------------------------------------
-template std::ostream& util::operator<< (std::ostream&, quaternion<4,float>);
-template std::ostream& util::operator<< (std::ostream&, quaternion<4,double>);
-
-
 ///////////////////////////////////////////////////////////////////////////////
-namespace util { namespace debug {
-    template <size_t S, typename T>
-    struct validator<quaternion<S,T>> {
-        static bool is_valid (const quaternion<S,T> &q)
+namespace util::debug {
+    template <typename T>
+    struct validator<quaternion<T>> {
+        static constexpr
+        bool
+        is_valid (const quaternion<T> &q)
         {
             return is_normalised (q);
         }
     };
-} }
-
-
-//-----------------------------------------------------------------------------
-template bool util::debug::is_valid(const quaternion<4,float>&);
-template bool util::debug::is_valid(const quaternion<4,double>&);
+}
 
 
 ///////////////////////////////////////////////////////////////////////////////
-template struct util::quaternion<4,float>;
-template struct util::quaternion<4,double>;
+#define INSTANTIATE(T)                                                          \
+template util::vector3<T> util::rotate (util::vector3<T>, util::quaternion<T>);   \
+template quaternion<T> util::conjugate (quaternion<T>);                         \
+template quaternion<T> util::operator* (quaternion<T>, quaternion<T>);          \
+template quaternion<T> util::operator/ (quaternion<T>, quaternion<T>);          \
+template bool util::almost_equal (util::quaternion<T>, util::quaternion<T>);    \
+template std::ostream& util::operator<< (std::ostream&, quaternion<T>);         \
+template bool util::debug::is_valid(const quaternion<T>&);                      \
+template struct util::quaternion<T>;
+
+INSTANTIATE(float)
+INSTANTIATE(double)
diff --git a/quaternion.hpp b/quaternion.hpp
index 58f8193d..362313a3 100644
--- a/quaternion.hpp
+++ b/quaternion.hpp
@@ -11,7 +11,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  *
- * Copyright 2011 Danny Robson <danny@nerdcruft.net>
+ * Copyright 2011-2016 Danny Robson <danny@nerdcruft.net>
  */
 
 #ifndef __UTIL_QUATERNION_HPP
@@ -26,20 +26,17 @@
 
 
 namespace util {
-    // quaternions must be 4 elements, but we include a size parameter so it
-    // fits with the generic coord infrastructure more easily.
+    // quaternion's are _just_ different enough to other coord types that we
+    // special case as a distinct POD type and provide many of the same
+    // functions as distinct declarations.
     //
-    // specifically:
-    // large regions of base code require a template template parameter with
-    // size and type arguments, which is annoying to work around for this one
-    // case.
-    //
-    // we protect against invalid instantiations through static_assert
-    template <size_t S, typename T>
-    struct quaternion : public coord::base<4,T,quaternion,coord::wxyz,coord::abcd> {
-        static_assert (S == 4, "quaternions must be 4 elements");
-
-        using coord::base<S,T,::util::quaternion,::util::coord::wxyz,::util::coord::abcd>::base;
+    // issues include:
+    // * strictly 4 dimensions
+    // * scalar operations sometimes don't make sense on the w component
+    // * objects must be normalised to make sense
+    template <typename T>
+    struct quaternion {
+        T w, x, y, z;
 
         static quaternion angle_axis (T radians, vector<3,T> axis);
         static quaternion from_euler (vector<3,T>);
@@ -49,32 +46,81 @@ namespace util {
 
         matrix4<T> as_matrix (void) const;
 
-        static const quaternion IDENTITY;
+        static constexpr quaternion<T> identity (void);
     };
 
-    template <typename T>
-    quaternion<4,T>
-    conjugate (quaternion<4,T>);
-
-    template <size_t S, typename T>
-    quaternion<S,T>
-    operator* (const quaternion<S,T>, const quaternion<S,T>);
-
-    template <size_t S, typename T>
-    quaternion<S,T>
-    operator/ (const quaternion<S,T>, const quaternion<S,T>);
-
-    typedef quaternion<4,float>  quaternionf;
-    typedef quaternion<4,double> quaterniond;
 
+    ///////////////////////////////////////////////////////////////////////////
     template <typename T>
     vector3<T>
-    rotate (vector3<T>, quaternion<4,T>);
+    rotate (vector3<T>, quaternion<T>);
 
-    template <size_t S, typename T>
+
+    ///////////////////////////////////////////////////////////////////////////
+    template <typename T>
+    constexpr
+    T
+    norm2 (quaternion<T>);
+
+    template <typename T>
+    constexpr
+    T
+    norm (quaternion<T>);
+
+    template <typename T>
+    constexpr
+    bool
+    is_normalised (quaternion<T>);
+
+    template <typename T>
+    constexpr
+    quaternion<T>
+    normalised (quaternion<T>);
+
+
+    ///////////////////////////////////////////////////////////////////////////
+    template <typename T>
+    quaternion<T>
+    conjugate (quaternion<T>);
+
+
+    ///////////////////////////////////////////////////////////////////////////
+    template <typename T>
+    quaternion<T>
+    operator* (quaternion<T>, quaternion<T>);
+
+    //-------------------------------------------------------------------------
+    template <typename T>
+    quaternion<T>
+    operator/ (quaternion<T>, quaternion<T>);
+
+    //-------------------------------------------------------------------------
+    template <typename T>
+    constexpr
+    quaternion<T>
+    operator/ (quaternion<T>, T);
+
+
+    ///////////////////////////////////////////////////////////////////////////
+    template <typename T>
+    constexpr
+    bool operator== (quaternion<T>, quaternion<T>);
+
+    //-------------------------------------------------------------------------
+    template <typename T>
+    bool almost_equal (quaternion<T>, quaternion<T>);
+
+    ///////////////////////////////////////////////////////////////////////////
+    typedef quaternion<float>  quaternionf;
+    typedef quaternion<double> quaterniond;
+
+
+    ///////////////////////////////////////////////////////////////////////////
+    template <typename T>
     std::ostream&
-    operator<< (std::ostream&, quaternion<S,T>);
+    operator<< (std::ostream&, quaternion<T>);
 }
 
+#include "./quaternion.ipp"
 
 #endif
diff --git a/quaternion.ipp b/quaternion.ipp
new file mode 100644
index 00000000..93e4e11f
--- /dev/null
+++ b/quaternion.ipp
@@ -0,0 +1,99 @@
+/*
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ * Copyright 2011-2016 Danny Robson <danny@nerdcruft.net>
+ */
+
+#if defined(CRUFT_UTIL_QUATERNION_IPP)
+#error
+#endif
+
+#define CRUFT_UTIL_QUATERNION_IPP
+
+#include <cmath>
+
+
+///////////////////////////////////////////////////////////////////////////////
+template <typename T>
+constexpr
+util::quaternion<T>
+util::quaternion<T>::identity (void)
+{
+    return { 1, 0, 0, 0 };
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+template <typename T>
+constexpr
+T
+util::norm2 (quaternion<T> q)
+{
+    return q.w * q.w +
+           q.x * q.x +
+           q.y * q.y +
+           q.z * q.z;
+}
+
+
+//-----------------------------------------------------------------------------
+template <typename T>
+constexpr
+T
+util::norm (quaternion<T> q)
+{
+    return std::sqrt (norm2 (q));
+}
+
+
+//-----------------------------------------------------------------------------
+template <typename T>
+constexpr
+bool
+util::is_normalised (quaternion<T> q)
+{
+    return almost_equal (T{1}, norm2 (q));
+}
+
+
+//-----------------------------------------------------------------------------
+template <typename T>
+constexpr
+util::quaternion<T>
+util::normalised (quaternion<T> q)
+{
+    return q / norm (q);
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+template <typename T>
+constexpr
+util::quaternion<T>
+util::operator/ (quaternion<T> q, T t)
+{
+    return { q.w / t, q.x / t, q.y / t, q.z / t };
+}
+
+
+///////////////////////////////////////////////////////////////////////////////
+template <typename T>
+constexpr
+bool
+util::operator== (quaternion<T> a, quaternion<T> b)
+{
+    return exactly_equal (a.w, b.w) &&
+           exactly_equal (a.x, b.x) &&
+           exactly_equal (a.y, b.y) &&
+           exactly_equal (a.z, b.z);
+}
diff --git a/test/quaternion.cpp b/test/quaternion.cpp
index 8cd5d677..eba9c9a0 100644
--- a/test/quaternion.cpp
+++ b/test/quaternion.cpp
@@ -17,22 +17,22 @@ main (void)
 
     // identity relations
     tap.expect_eq (
-        norm (quaternionf::IDENTITY), 1.f,
+        norm (quaternionf::identity ()), 1.f,
         "identity magnitude is unit"
     );
 
     tap.expect_eq (
-        quaternionf::IDENTITY * quaternionf::IDENTITY,
-        quaternionf::IDENTITY,
+        quaternionf::identity () * quaternionf::identity (),
+        quaternionf::identity (),
         "identity multiplication with identity"
     );
 
     // normalisation
     {
-        auto val = normalised (quaternionf (2, 3, 4, 7));
+        auto val = normalised (quaternionf {2, 3, 4, 7});
 
         tap.expect_eq (
-            val * quaternionf::IDENTITY,
+            val * quaternionf::identity (),
             val,
             "identity multiplication with quaternion constant"
         );
@@ -47,9 +47,11 @@ main (void)
         // towards rotations than general maths).
         util::vector4f a_v { 2, -11,  5, -17};
         util::vector4f b_v { 3,  13, -7, -19};
+        a_v = normalised (a_v);
+        b_v = normalised (b_v);
 
-        auto a = normalised (a_v).as<quaternion> ();
-        auto b = normalised (b_v).as<quaternion> ();
+        auto a = quaternionf { a_v[0], a_v[1], a_v[2], a_v[3] };
+        auto b = quaternionf { b_v[0], b_v[1], b_v[2], b_v[3] };
         auto c = quaternionf {
             -0.27358657116960006f,
             -0.43498209092420004f,
@@ -62,7 +64,7 @@ main (void)
     }
 
     tap.expect_eq (
-        quaternionf::IDENTITY.as_matrix (),
+        quaternionf::identity ().as_matrix (),
         util::matrix4f::IDENTITY,
         "identity quaternion to matrix"
     );
@@ -90,7 +92,7 @@ main (void)
             tap.expect_lt (util::sum (diff), 1e-6f, "single basis rotation %zu", i);
         }
 
-        auto q = quaternionf::IDENTITY;
+        auto q = quaternionf::identity ();
         auto m = util::matrix4f::IDENTITY;
 
         for (auto r: ROTATIONS) {