Math: dual number implementation.

src/Math/CMakeLists.txt

@@ -1,6 +1,7 @@
 set(MagnumMath_HEADERS
     BoolVector.h
     Constants.h
+    Dual.h
     Functions.h
     Math.h
     MathTypeTraits.h

src/Math/Dual.h

@@ -0,0 +1,167 @@
+#ifndef Magnum_Math_Dual_h
+#define Magnum_Math_Dual_h
+/*
+    Copyright © 2010, 2011, 2012 Vladimír Vondruš <mosra@centrum.cz>
+
+    This file is part of Magnum.
+
+    Magnum is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser General Public License version 3
+    only, as published by the Free Software Foundation.
+
+    Magnum is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU Lesser General Public License version 3 for more details.
+*/
+
+/** @file
+ * @brief Class Magnum::Math::Dual
+ */
+
+#include <Utility/Debug.h>
+
+#include "Math/MathTypeTraits.h"
+
+namespace Magnum { namespace Math {
+
+/** @brief %Dual number */
+template<class T> class Dual {
+    public:
+        /**
+         * @brief Default constructor
+         *
+         * Both parts are default-constructed.
+         */
+        inline constexpr /*implicit*/ Dual(): _real(), _dual() {}
+
+        /**
+         * @brief Construct dual number from real and dual part
+         *
+         * @f[
+         *      \hat a = a_0 + \epsilon a_\epsilon
+         * @f]
+         */
+        inline constexpr /*implicit*/ Dual(const T& real, const T& dual = T()): _real(real), _dual(dual) {}
+
+        /** @brief Equality comparison */
+        inline bool operator==(const Dual<T>& other) const {
+            return MathTypeTraits<T>::equals(_real, other._real) &&
+                   MathTypeTraits<T>::equals(_dual, other._dual);
+        }
+
+        /** @brief Non-equality comparison */
+        inline bool operator!=(const Dual<T>& other) const {
+            return !operator==(other);
+        }
+
+        /** @brief Real part */
+        inline constexpr T real() const { return _real; }
+
+        /** @brief %Dual part */
+        inline constexpr T dual() const { return _dual; }
+
+        /**
+         * @brief Add and assign dual number
+         *
+         * The computation is done in-place. @f[
+         *      \hat a + \hat b = a_0 + b_0 + \epsilon (a_\epsilon + b_\epsilon)
+         * @f]
+         */
+        inline Dual<T>& operator+=(const Dual<T>& other) {
+            _real += other._real;
+            _dual += other._dual;
+            return *this;
+        }
+
+        /**
+         * @brief Add dual number
+         *
+         * @see operator+=()
+         */
+        inline Dual<T> operator+(const Dual<T>& other) const {
+            return Dual<T>(*this)+=other;
+        }
+
+        /**
+         * @brief Negated dual number
+         *
+         * @f[
+         *      -\hat a = -a_0 - \epsilon a_\epsilon
+         * @f]
+         */
+        inline Dual<T> operator-() const {
+            return {-_real, -_dual};
+        }
+
+        /**
+         * @brief Subtract and assign dual number
+         *
+         * The computation is done in-place. @f[
+         *      \hat a - \hat b = a_0 - b_0 + \epsilon (a_\epsilon - b_\epsilon)
+         * @f]
+         */
+        inline Dual<T>& operator-=(const Dual<T>& other) {
+            _real -= other._real;
+            _dual -= other._dual;
+            return *this;
+        }
+
+        /**
+         * @brief Subtract dual number
+         *
+         * @see operator-=()
+         */
+        inline Dual<T> operator-(const Dual<T>& other) const {
+            return Dual<T>(*this)-=other;
+        }
+
+        /**
+         * @brief Multiply by dual number
+         *
+         * @f[
+         *      \hat a \hat b = a_0 b_0 + \epsilon (a_0 b_\epsilon + a_\epsilon b_0)
+         * @f]
+         */
+        inline Dual<T> operator*(const Dual<T>& other) const {
+            return {_real*other._real, _real*other._dual + _dual*other._real};
+        }
+
+        /**
+         * @brief Divide by dual number
+         *
+         * @f[
+         *      \frac{\hat a}{\hat b} = \frac{a_0}{b_0} + \epsilon \frac{a_\epsilon b_0 - a_0 b_\epsilon}{b_0^2}
+         * @f]
+         */
+        inline Dual<T> operator/(const Dual<T>& other) const {
+            return {_real/other._real, (_dual*other._real - _real*other._dual)/(other._real*other._real)};
+        }
+
+        /**
+         * @brief Conjugated dual number
+         *
+         * @f[
+         *      \overline{\hat a} = a_0 - \epsilon a_\epsilon
+         * @f]
+         */
+        inline Dual<T> conjugated() const {
+            return {_real, -_dual};
+        }
+
+    private:
+        T _real, _dual;
+};
+
+/** @debugoperator{Magnum::Math::Dual} */
+template<class T> Corrade::Utility::Debug operator<<(Corrade::Utility::Debug debug, const Dual<T>& value) {
+    debug << "Dual(";
+    debug.setFlag(Corrade::Utility::Debug::SpaceAfterEachValue, false);
+    debug << value.real() << ", " << value.dual() << ")";
+    debug.setFlag(Corrade::Utility::Debug::SpaceAfterEachValue, true);
+    return debug;
+}
+
+}}
+
+#endif

src/Math/Math.h

@@ -27,6 +27,8 @@ namespace Magnum { namespace Math {
 #ifndef DOXYGEN_GENERATING_OUTPUT
 /* Class Constants used only statically */
 
+template<class T> class Dual;
+
 template<std::size_t, class> class Matrix;
 template<class> class Matrix3;
 template<class> class Matrix4;

src/Math/Test/CMakeLists.txt

@@ -17,6 +17,7 @@ corrade_add_test(MathMatrix4Test Matrix4Test.cpp LIBRARIES MagnumMathTestLib)
 
 corrade_add_test(MathSwizzleTest SwizzleTest.cpp LIBRARIES MagnumMathTestLib)
 
+corrade_add_test(MathDualTest DualTest.cpp)
 corrade_add_test(MathQuaternionTest QuaternionTest.cpp LIBRARIES MagnumMathTestLib)
 
 set_target_properties(MathVectorTest MathMatrix3Test MathMatrix4Test MathQuaternionTest PROPERTIES COMPILE_FLAGS -DCORRADE_GRACEFUL_ASSERT)

src/Math/Test/DualTest.cpp

@@ -0,0 +1,116 @@
+/*
+    Copyright © 2010, 2011, 2012 Vladimír Vondruš <mosra@centrum.cz>
+
+    This file is part of Magnum.
+
+    Magnum is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser General Public License version 3
+    only, as published by the Free Software Foundation.
+
+    Magnum is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU Lesser General Public License version 3 for more details.
+*/
+
+#include <sstream>
+#include <TestSuite/Tester.h>
+
+#include "Math/Dual.h"
+
+namespace Magnum { namespace Math { namespace Test {
+
+class DualTest: public Corrade::TestSuite::Tester {
+    public:
+        explicit DualTest();
+
+        void construct();
+        void constructDefault();
+        void compare();
+
+        void addSubtract();
+        void negated();
+        void multiplyDivide();
+
+        void conjugated();
+
+        void debug();
+};
+
+typedef Math::Dual<float> Dual;
+
+DualTest::DualTest() {
+    addTests(&DualTest::construct,
+             &DualTest::constructDefault,
+             &DualTest::compare,
+
+             &DualTest::addSubtract,
+             &DualTest::negated,
+             &DualTest::multiplyDivide,
+
+             &DualTest::conjugated,
+
+             &DualTest::debug);
+}
+
+void DualTest::construct() {
+    Dual a(2.0f, -7.5f);
+    CORRADE_COMPARE(a.real(), 2.0f);
+    CORRADE_COMPARE(a.dual(), -7.5f);
+
+    Dual b(3.0f);
+    CORRADE_COMPARE(b.real(), 3.0f);
+    CORRADE_COMPARE(b.dual(), 0.0f);
+}
+
+void DualTest::constructDefault() {
+    CORRADE_COMPARE(Dual(), Dual(0.0f, 0.0f));
+}
+
+void DualTest::compare() {
+    CORRADE_VERIFY(Dual(1.0f, 1.0f+MathTypeTraits<float>::epsilon()/2) == Dual(1.0f, 1.0f));
+    CORRADE_VERIFY(Dual(1.0f, 1.0f+MathTypeTraits<float>::epsilon()*2) != Dual(1.0f, 1.0f));
+    CORRADE_VERIFY(Dual(1.0f+MathTypeTraits<float>::epsilon()/2, 1.0f) == Dual(1.0f, 1.0f));
+    CORRADE_VERIFY(Dual(1.0f+MathTypeTraits<float>::epsilon()*2, 1.0f) != Dual(1.0f, 1.0f));
+
+    /* Compare to real part only */
+    CORRADE_VERIFY(Dual(1.0f, 0.0f) == 1.0f);
+    CORRADE_VERIFY(Dual(1.0f, 3.0f) != 1.0f);
+}
+
+void DualTest::addSubtract() {
+    Dual a(2.0f, -7.5f);
+    Dual b(-3.3f, 0.2f);
+    Dual c(-1.3f, -7.3f);
+
+    CORRADE_COMPARE(a + b, c);
+    CORRADE_COMPARE(c - b, a);
+}
+
+void DualTest::negated() {
+    CORRADE_COMPARE(-Dual(1.0f, -6.5f), Dual(-1.0f, 6.5f));
+}
+
+void DualTest::multiplyDivide() {
+    Dual a(1.5f, -4.0f);
+    Dual b(-2.0f, 0.5f);
+    Dual c(-3.0f, 8.75f);
+
+    CORRADE_COMPARE(a*b, c);
+    CORRADE_COMPARE(c/b, a);
+}
+
+void DualTest::conjugated() {
+    CORRADE_COMPARE(Dual(1.0f, -6.5f).conjugated(), Dual(1.0f, 6.5f));
+}
+
+void DualTest::debug() {
+    std::ostringstream o;
+
+    Debug(&o) << Dual(2.5f, -0.3f);
+    CORRADE_COMPARE(o.str(), "Dual(2.5, -0.3)\n");
+}
+
+}}}
+
+CORRADE_TEST_MAIN(Magnum::Math::Test::DualTest)