Dual.hpp

Go to the documentation of this file.

#pragma once
 
#include <Math/math/Functions.hpp>
#include <cassert>
#include <concepts>
#include <type_traits>
 
 
namespace Math
{
 
template<class T>
class BasicDual
{
public:
    using value_type = T; 
 
    constexpr BasicDual() = default; 
    
    constexpr BasicDual(const T &r) : real(r) { }
 
    constexpr BasicDual(const T &r, const T &d) : real(r), dual(d) { }
 
    constexpr static BasicDual<T> identity() { return { T{1}, T{0} }; }
    constexpr static BasicDual<T> zero() { return { }; }
 
    constexpr BasicDual<T> conjugate() const
    {
        if constexpr (std::is_floating_point_v<T>)
            return { real, -dual };
        else
            return { real, conjugate(dual) };
    }
 
    constexpr T magnitude() const { return real; }
 
    constexpr static BasicDual<T> make_pure(const T &input)
    {
        return { T{}, input };
    }
 
    bool isNaN() const
    {
        if constexpr ( std::is_floating_point_v<T> )
            return std::isnan(real) || std::isnan(dual);
        else
            return real.isNaN() || dual.isNaN();
    }
 
    bool isInf() const
    {
        if constexpr ( std::is_floating_point_v<T> )
            return std::isinf(real) || std::isinf(dual);
        else
            return real.isInf() || dual.isInf();
    }
 
    T real{};
    T dual{};
 
private:
 
    friend constexpr bool operator ==(const BasicDual<T> &left, const BasicDual<T> &right)
    {
        return approximately_equal_to(left, right);
    }
 
    template <std::floating_point OT = T>
        requires std::is_floating_point_v<T>
    friend constexpr bool approximately_equal_to(const BasicDual<T> &value_to_test,
                                                 const BasicDual<T> &value_it_should_be,
                                                       OT            tolerance = OT{0.0002})
    {
        return approximately_equal_to(value_to_test.real, value_it_should_be.real, tolerance) &&
               approximately_equal_to(value_to_test.dual, value_it_should_be.dual, tolerance);
    }
 
    friend constexpr BasicDual<T> operator *(const BasicDual<T> &left, const BasicDual<T> &right)
    {
        return { left.real * right.real,
                 left.real * right.dual + right.real * left.dual };
    }
 
    template <std::floating_point OT = double>
    friend constexpr BasicDual<T> operator *(OT scalar, const BasicDual<T> &d)
    {
        return BasicDual<T>( T(scalar) ) * d;
    }
 
    template <std::floating_point OT = double>
    friend constexpr BasicDual<T> operator *(const BasicDual<T> &d, OT scalar)
    {
        return d * BasicDual<T>( T(scalar) );
    }
 
    friend constexpr BasicDual<T> operator /(const BasicDual<T> &left, const BasicDual<T> &right)
    {
        return { left.real * right.real / (right.real * right.real),
                (left.dual * right.real - left.real * right.dual) / (right.real * right.real) };
    }
 
    template <std::floating_point OT = double>
    friend constexpr BasicDual<T> operator /(OT scalar, const BasicDual<T> &d)
    {
        return BasicDual<T>( T(scalar) ) / d;
    }
 
    template <std::floating_point OT = double>
    friend constexpr BasicDual<T> operator /(const BasicDual<T> &left, OT scalar)
    {
        return left / BasicDual<T>( T(scalar) );
    }
 
    friend constexpr BasicDual<T> operator +(const BasicDual<T> &left, const BasicDual<T> &right)
    {
        return { left.real + right.real, left.dual + right.dual };
    }
 
    template <std::floating_point OT = double>
    friend constexpr BasicDual<T> operator +(OT scalar, const BasicDual<T> &d)
    {
        return BasicDual<T>( T(scalar) ) + d;
    }
 
    template <std::floating_point OT = double>
    friend constexpr BasicDual<T> operator +(const BasicDual<T> &left, OT scalar)
    {
        return left + BasicDual<T>( T(scalar) );
    }
 
    friend constexpr BasicDual<T> operator -(const BasicDual<T> &left, const BasicDual<T> &right)
    {
        return { left.real - right.real, left.dual - right.dual };
    }
 
    template <std::floating_point OT = double>
    friend constexpr BasicDual<T> operator -(OT scalar, const BasicDual<T> &d)
    {
        return BasicDual<T>( T(scalar) ) - d;
    }
 
    template <std::floating_point OT = double>
    friend constexpr BasicDual<T> operator -(const BasicDual<T> &left, OT scalar)
    {
        return left - BasicDual<T>( T(scalar) );
    }
 
    template <std::floating_point OT = T>
        requires std::is_floating_point_v<T>
    friend bool check_if_equal(const BasicDual<T> &input,
                               const BasicDual<T> &near_to,
                                     OT            tolerance = OT{0.0002})
    {
        if (!approximately_equal_to(input, near_to, tolerance))
        {
            auto diff{ near_to - input };
 
            std::cout << std::format("input: {} is not equal to near_to: {} within tolerance: {}.  Difference is {} .",
                                     format(input),
                                     format(near_to),
                                     tolerance,
                                     format(near_to - input))
            << std::endl;
            return  false;
        }
        return true;
    }
 
    template <std::floating_point OT = T>
        requires std::is_floating_point_v<T>
    friend bool check_if_not_equal(const BasicDual<T> &input,
                                   const BasicDual<T> &near_to,
                                         OT            tolerance = OT{0.0002})
    {
        if (approximately_equal_to(input, near_to, tolerance))
        {
            auto diff{ near_to - input };
 
            std::cout << std::format("input: {} is equal to near_to: {} within tolerance: {}.  Difference is {} .",
                                     format(input),
                                     format(near_to),
                                     tolerance,
                                     format(near_to - input))
            << std::endl;
            return  false;
        }
        return true;
    }
 
    template <std::floating_point OT = T>
        requires std::is_floating_point_v<T>
    friend void CHECK_IF_EQUAL(const BasicDual<T> &input,
                               const BasicDual<T> &near_to,
                                     OT            tolerance = OT{0.0002})
    {
        assert( check_if_equal(input, near_to, tolerance) );
    }
 
    template <std::floating_point OT = T>
        requires std::is_floating_point_v<T>
    friend void CHECK_IF_NOT_EQUAL(const BasicDual<T> &input,
                                   const BasicDual<T> &near_to,
                                         OT            tolerance = OT{0.0002})
    {
        assert( check_if_not_equal(input, near_to, tolerance) );
    }
 
    template <std::floating_point OT = T>
        requires std::is_floating_point_v<T>
    friend void CHECK_IF_ZERO(const BasicDual<T> &input, OT tolerance = OT{0.0002})
    {
        assert( check_if_equal(input, BasicDual<T>::zero(), tolerance));
    }
 
    friend constexpr T dot(const BasicDual<T> &left, const BasicDual<T> &right)
    {
        return left.real * right.real +
               left.dual * right.dual;
    }
 
    friend constexpr BasicDual<T> dualscalar_sqrt(const BasicDual<T> &input)
    {
        // Expect that T is a scalar type (float, double, int, etc.)
        T root = std::sqrt(input.real);
 
        return { root, input.dual / (T{2} * root) };
    }
 
    friend constexpr T dualscalar_normsquared(const BasicDual<T> &d)
    {
        BasicDual<T> result = d * d.conjugate();
 
        assert( approximately_equal_to(result.dual, T{0}) );
 
        return result.real;
    }
 
    friend constexpr T accumulate(const BasicDual<T> &input)
    {
        return T{input.real + input.dual};
    }
 
    friend constexpr BasicDual<T> conjugate(const BasicDual<T> &input)
    {
        return input.conjugate();
    }
 
    friend std::string format(const BasicDual<T> &input)
    {
        return std::format("[real: {}, dual: {}]", input.real, input.dual);
    }
};
 
 
using Dualf = BasicDual<float>;
using Duald = BasicDual<double>;
using Dual  = BasicDual<double>;
using Duall = BasicDual<long double>;
 
}