lisyarus/psemek
1
0
Fork 1
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Add n-dimensional dual numbers implementation

Browse source
This commit is contained in:
Nikita Lisitsa 2021-12-01 21:23:13 +03:00
parent ec86b0c0aa
commit 93f4147941
1 changed files with 160 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

160
libs/geom/include/psemek/geom/dual.hpp Normal file
View file

@ -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;
}
}