From 93f414794166703088773bc39f1171ecd3601e6d Mon Sep 17 00:00:00 2001
From: lisyarus <lisyarus@gmail.com>
Date: Wed, 1 Dec 2021 21:23:13 +0300
Subject: [PATCH] Add n-dimensional dual numbers implementation

---
 libs/geom/include/psemek/geom/dual.hpp | 160 +++++++++++++++++++++++++
 1 file changed, 160 insertions(+)
 create mode 100644 libs/geom/include/psemek/geom/dual.hpp

diff --git a/libs/geom/include/psemek/geom/dual.hpp b/libs/geom/include/psemek/geom/dual.hpp
new file mode 100644
index 00000000..d1b190aa
--- /dev/null
+++ b/libs/geom/include/psemek/geom/dual.hpp
@@ -0,0 +1,160 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+
+#include <iostream>
+
+namespace psemek::geom
+{
+
+	template <typename T, std::size_t N>
+	struct dual
+	{
+		T scalar = T{0};
+		vector<T, N> delta = vector<T, N>::zero();
+
+		dual() = default;
+
+		dual(T scalar)
+			: scalar{std::move(scalar)}
+		{}
+
+		dual(T scalar, vector<T, N> delta)
+			: scalar{std::move(scalar)}
+			, delta{std::move(delta)}
+		{}
+
+		dual & operator += (dual const & d);
+		dual & operator -= (dual const & d);
+		dual & operator *= (dual const & d);
+		dual & operator /= (dual const & d);
+	};
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator + (dual<T, N> const & d)
+	{
+		return d;
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator - (dual<T, N> const & d)
+	{
+		return {-d.scalar, -d.delta};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator + (dual<T, N> const & d1, dual<T, N> const & d2)
+	{
+		return {d1.scalar + d2.scalar, d1.delta + d2.delta};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator - (dual<T, N> const & d1, dual<T, N> const & d2)
+	{
+		return {d1.scalar - d2.scalar, d1.delta - d2.delta};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator * (dual<T, N> const & d1, dual<T, N> const & d2)
+	{
+		return {d1.scalar * d2.scalar, d1.scalar * d2.delta + d2.scalar * d1.delta};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> inverse(dual<T, N> const & d)
+	{
+		return {T{1} / d.scalar, - d.delta / d.scalar / d.scalar};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator / (dual<T, N> const & d1, dual<T, N> const & d2)
+	{
+		return d1 * inverse(d2);
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> & dual<T, N>::operator += (dual<T, N> const & d)
+	{
+		*this = *this + d;
+		return *this;
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> & dual<T, N>::operator -= (dual const & d)
+	{
+		*this = *this - d;
+		return *this;
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> & dual<T, N>::operator *= (dual const & d)
+	{
+		*this = *this * d;
+		return *this;
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> & dual<T, N>::operator /= (dual const & d)
+	{
+		*this = *this / d;
+		return *this;
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator + (dual<T, N> const & d, T const & v)
+	{
+		return {d.scalar + v, d.delta};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator + (T const & v, dual<T, N> const & d)
+	{
+		return {d.scalar + v, d.delta};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator - (dual<T, N> const & d, T const & v)
+	{
+		return {d.scalar - v, d.delta};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator - (T const & v, dual<T, N> const & d)
+	{
+		return {v - d.scalar, -d.delta};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator * (dual<T, N> const & d, T const & v)
+	{
+		return {d.scalar * v, d.delta * v};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator * (T const & v, dual<T, N> const & d)
+	{
+		return {d.scalar * v, d.delta * v};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator / (dual<T, N> const & d, T const & v)
+	{
+		return {d.scalar / v, d.delta / v};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> operator / (T const & v, dual<T, N> const & d)
+	{
+		return v * inverse(d);
+	}
+
+	template <typename T, std::size_t N>
+	std::ostream & operator << (std::ostream & os, dual<T, N> const & d)
+	{
+		os << d.scalar;
+		for (std::size_t i = 0; i < N; ++i)
+			os << " + " << d.delta[i] << " e_" << i;
+		return os;
+	}
+
+}