extern crate nalgebra as na;
extern crate num_traits;
use std::cmp::Ordering;
use std::fmt::{Debug, Display, Formatter, LowerExp, Result as FmtResult};
use std::iter::{Product, Sum};
use std::num::FpCategory;
use std::ops::{Add, AddAssign, Deref, DerefMut, Div, DivAssign, Mul, MulAssign, Neg, Rem, RemAssign, Sub, SubAssign};
pub use num_traits::{Float, FloatConst, Num, One, Zero};
mod differentials;
pub use differentials::*;
pub mod linalg;
use num_traits::{FromPrimitive, Inv, MulAdd, MulAddAssign, NumCast, Pow, Signed, ToPrimitive, Unsigned};
use na::{OVector, SVector, Scalar};
pub use na::allocator::Allocator;
pub use na::dimension::*;
pub use na::storage::Owned;
pub use na::{DefaultAllocator, Dim, DimName};
#[derive(Clone, Copy)]
pub struct OHyperdual<T: Copy + Scalar, N: Dim + DimName>(OVector<T, N>)
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy;
pub type Hyperdual<T, const N: usize> = OHyperdual<T, Const<N>>;
impl<T: Copy + Scalar, N: Dim + DimName> OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
pub fn from_slice(v: &[T]) -> Self {
Self(OVector::<T, N>::from_row_slice(v))
}
#[inline]
pub fn from_real(real: T) -> Self
where
T: Zero,
{
let mut dual = OVector::<T, N>::zeros();
dual[0] = real;
Self(dual)
}
#[inline]
pub fn real(&self) -> T {
self[0]
}
#[inline]
pub fn real_ref(&self) -> &T {
&self[0]
}
pub fn real_mut(&mut self) -> &mut T {
&mut self[0]
}
#[inline]
pub fn map_dual<F>(&self, r: T, f: F) -> Self
where
F: Fn(&T) -> T,
{
let mut v = self.map(|x| f(&x));
v[0] = r;
Self(v)
}
#[inline]
pub fn from_fn<F>(mut f: F) -> Self
where
F: FnMut(usize) -> T,
{
Self(OVector::<T, N>::from_fn(|i, _| f(i)))
}
}
impl<T: Copy + Scalar, N: Dim + DimName> Debug for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
fn fmt(&self, f: &mut Formatter) -> FmtResult {
let mut a = f.debug_tuple("Dual");
for x in self.iter() {
a.field(x);
}
a.finish()
}
}
impl<T: Copy + Scalar + Num + Zero, N: Dim + DimName> Default for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn default() -> Self {
Self::zero()
}
}
impl<T: Copy + Scalar + Zero, N: Dim + DimName> From<T> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn from(real: T) -> Self {
Self::from_real(real)
}
}
impl<T: Copy + Scalar, N: Dim + DimName> Deref for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Target = OVector<T, N>;
#[inline]
fn deref(&self) -> &OVector<T, N> {
&self.0
}
}
impl<T: Copy + Scalar, N: Dim + DimName> DerefMut for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn deref_mut(&mut self) -> &mut OVector<T, N> {
&mut self.0
}
}
impl<T: Copy + Scalar, N: Dim + DimName> AsRef<OVector<T, N>> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn as_ref(&self) -> &OVector<T, N> {
&self.0
}
}
impl<T: Copy + Scalar, N: Dim + DimName> AsMut<OVector<T, N>> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn as_mut(&mut self) -> &mut OVector<T, N> {
&mut self.0
}
}
impl<T: Copy + Scalar + Neg<Output = T>, N: Dim + DimName> OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
pub fn conjugate(self) -> Self {
self.map_dual(self.real(), |x| x.neg())
}
}
impl<T: Copy + Scalar + Display, N: Dim + DimName> Display for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
fn fmt(&self, f: &mut Formatter) -> FmtResult {
let precision = f.precision().unwrap_or(4);
write!(f, "{:.p$}", self.real(), p = precision)?;
for (i, x) in self.iter().skip(1).enumerate() {
write!(f, " + {:.p$}\u{03B5}{}", x, i + 1, p = precision)?;
}
Ok(())
}
}
impl<T: Copy + Scalar + LowerExp, N: Dim + DimName> LowerExp for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
fn fmt(&self, f: &mut Formatter) -> FmtResult {
let precision = f.precision().unwrap_or(4);
write!(f, "{:.p$e}", self.real(), p = precision)?;
for (i, x) in self.iter().skip(1).enumerate() {
write!(f, " + {:.p$e}\u{03B5}{}", x, i + 1, p = precision)?;
}
Ok(())
}
}
impl<T: Copy + Scalar + PartialEq, N: Dim + DimName> PartialEq<Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn eq(&self, rhs: &Self) -> bool {
self.0 == rhs.0
}
}
impl<T: Copy + Scalar + PartialOrd, N: Dim + DimName> PartialOrd<Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn partial_cmp(&self, rhs: &Self) -> Option<Ordering> {
PartialOrd::partial_cmp(self.real_ref(), rhs.real_ref())
}
#[inline]
fn lt(&self, rhs: &Self) -> bool {
self.real() < rhs.real()
}
#[inline]
fn gt(&self, rhs: &Self) -> bool {
self.real() > rhs.real()
}
}
impl<T: Copy + Scalar + PartialEq, N: Dim + DimName> PartialEq<T> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn eq(&self, rhs: &T) -> bool {
*self.real_ref() == *rhs
}
}
impl<T: Copy + Scalar + PartialOrd, N: Dim + DimName> PartialOrd<T> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn partial_cmp(&self, rhs: &T) -> Option<Ordering> {
PartialOrd::partial_cmp(self.real_ref(), rhs)
}
#[inline]
fn lt(&self, rhs: &T) -> bool {
self.real() < *rhs
}
#[inline]
fn gt(&self, rhs: &T) -> bool {
self.real() > *rhs
}
}
macro_rules! impl_to_primitive {
($($name:ident, $ty:ty),*) => {
impl<T: Copy + Scalar + ToPrimitive, N: Dim + DimName> ToPrimitive for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
$(
#[inline]
fn $name(&self) -> Option<$ty> {
self.real_ref().$name()
}
)*
}
}
}
macro_rules! impl_from_primitive {
($($name:ident, $ty:ty),*) => {
impl<T: Copy + Scalar + FromPrimitive, N: Dim + DimName> FromPrimitive for OHyperdual<T, N>
where
T: Zero,
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
$(
#[inline]
fn $name(n: $ty) -> Option<OHyperdual<T,N>> {
T::$name(n).map(Self::from_real)
}
)*
}
}
}
macro_rules! impl_primitive_cast {
($($to:ident, $from:ident - $ty:ty),*) => {
impl_to_primitive!($($to, $ty),*);
impl_from_primitive!($($from, $ty),*);
}
}
impl_primitive_cast! {
to_isize, from_isize - isize,
to_i8, from_i8 - i8,
to_i16, from_i16 - i16,
to_i32, from_i32 - i32,
to_i64, from_i64 - i64,
to_usize, from_usize - usize,
to_u8, from_u8 - u8,
to_u16, from_u16 - u16,
to_u32, from_u32 - u32,
to_u64, from_u64 - u64,
to_f32, from_f32 - f32,
to_f64, from_f64 - f64
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Add<T> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = OHyperdual<T, N>;
#[inline]
fn add(self, rhs: T) -> Self {
let mut d = self;
d[0] = d[0] + rhs;
d
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> AddAssign<T> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn add_assign(&mut self, rhs: T) {
*self = (*self) + Self::from_real(rhs)
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Sub<T> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = OHyperdual<T, N>;
#[inline]
fn sub(self, rhs: T) -> Self {
let mut d = self;
d[0] = d[0] - rhs;
d
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> SubAssign<T> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn sub_assign(&mut self, rhs: T) {
*self = (*self) - Self::from_real(rhs)
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Mul<T> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = OHyperdual<T, N>;
#[inline]
fn mul(self, rhs: T) -> Self {
self * Self::from_real(rhs)
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> MulAssign<T> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn mul_assign(&mut self, rhs: T) {
*self = (*self) * Self::from_real(rhs)
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Div<T> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = OHyperdual<T, N>;
#[inline]
fn div(self, rhs: T) -> Self {
self / Self::from_real(rhs)
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> DivAssign<T> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn div_assign(&mut self, rhs: T) {
*self = (*self) / Self::from_real(rhs)
}
}
impl<T: Copy + Scalar + Signed, N: Dim + DimName> Neg for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
#[inline]
fn neg(self) -> Self {
Self(self.map(|x| x.neg()))
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Add<Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
#[inline]
fn add(self, rhs: Self) -> Self {
Self(self.zip_map(&rhs, |x, y| x + y))
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Add<&Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
#[inline]
fn add(self, rhs: &Self) -> Self {
Self(self.zip_map(rhs, |x, y| x + y))
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> AddAssign<Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn add_assign(&mut self, rhs: Self) {
*self = (*self) + rhs
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Sub<Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
#[inline]
fn sub(self, rhs: Self) -> Self {
Self(self.zip_map(&rhs, |x, y| x - y))
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Sub<&Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
#[inline]
fn sub(self, rhs: &Self) -> Self {
Self(self.zip_map(rhs, |x, y| x - y))
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> SubAssign<Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn sub_assign(&mut self, rhs: Self) {
*self = (*self) - rhs
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Mul<Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
#[allow(clippy::suspicious_arithmetic_impl)]
#[inline]
fn mul(self, rhs: Self) -> Self {
let mut v = self.zip_map(&rhs, |x, y| rhs.real() * x + self.real() * y);
v[0] = self.real() * rhs.real();
Self(v)
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Mul<&Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
#[allow(clippy::suspicious_arithmetic_impl)]
#[inline]
fn mul(self, rhs: &Self) -> Self {
let mut v = self.zip_map(rhs, |x, y| rhs.real() * x + self.real() * y);
v[0] = self.real() * rhs.real();
Self(v)
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> MulAssign<Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn mul_assign(&mut self, rhs: Self) {
*self = (*self) * rhs
}
}
macro_rules! impl_mul_add {
($(<$a:ident, $b:ident>),*) => {
$(
impl<T: Copy + Scalar + Num + Mul + Add, N: Dim + DimName> MulAdd<$a, $b> for OHyperdual<T,N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = OHyperdual<T,N>;
#[inline]
fn mul_add(self, a: $a, b: $b) -> Self {
(self * a) + b
}
}
impl<T: Copy + Scalar + Num + Mul + Add, N: Dim + DimName> MulAddAssign<$a, $b> for OHyperdual<T,N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn mul_add_assign(&mut self, a: $a, b: $b) {
*self = (*self * a) + b;
}
}
)*
}
}
impl_mul_add! {
<Self, Self>,
<T, Self>,
<Self, T>,
<T, T>
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Div<Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
#[inline]
#[allow(clippy::suspicious_arithmetic_impl)]
fn div(self, rhs: Self) -> Self {
let d = rhs.real() * rhs.real();
let mut v = self.zip_map(&rhs, |x, y| (rhs.real() * x - self.real() * y) / d);
v[0] = self.real() / rhs.real();
Self(v)
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Div<&Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
#[inline]
#[allow(clippy::suspicious_arithmetic_impl)]
fn div(self, rhs: &Self) -> Self {
let d = rhs.real() * rhs.real();
let mut v = self.zip_map(rhs, |x, y| (rhs.real() * x - self.real() * y) / d);
v[0] = self.real() / rhs.real();
Self(v)
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> DivAssign<Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn div_assign(&mut self, rhs: Self) {
*self = (*self) / rhs
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Rem<Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
fn rem(self, _: Self) -> Self {
unimplemented!()
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Rem<&Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
fn rem(self, _: &Self) -> Self {
unimplemented!()
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> RemAssign<Self> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
fn rem_assign(&mut self, _: Self) {
unimplemented!()
}
}
impl<T: Copy + Scalar, P: Into<OHyperdual<T, N>>, N: Dim + DimName> Pow<P> for OHyperdual<T, N>
where
OHyperdual<T, N>: Float,
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
#[inline]
fn pow(self, rhs: P) -> Self {
self.powf(rhs.into())
}
}
impl<T: Copy + Scalar, N: Dim + DimName> Inv for OHyperdual<T, N>
where
Self: One + Div<Output = Self>,
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type Output = Self;
#[inline]
fn inv(self) -> Self {
Self::one() / self
}
}
impl<T: Copy + Scalar, N: Dim + DimName> Signed for OHyperdual<T, N>
where
T: Signed + PartialOrd,
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn abs(&self) -> Self {
let s = self.real().signum();
Self(self.map(|x| x * s))
}
#[inline]
fn abs_sub(&self, rhs: &Self) -> Self {
if self.real() > rhs.real() {
self.sub(*rhs)
} else {
Self::zero()
}
}
#[inline]
fn signum(&self) -> Self {
Self::from_real(self.real().signum())
}
#[inline]
fn is_positive(&self) -> bool {
self.real().is_positive()
}
#[inline]
fn is_negative(&self) -> bool {
self.real().is_negative()
}
}
impl<T: Copy + Scalar + Unsigned, N: Dim + DimName> Unsigned for OHyperdual<T, N>
where
Self: Num,
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
}
impl<T: Copy + Scalar + Num + Zero, N: Dim + DimName> Zero for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn zero() -> Self {
Self::from_real(T::zero())
}
#[inline]
fn is_zero(&self) -> bool {
self.iter().all(|x| x.is_zero())
}
}
impl<T: Copy + Scalar + Num + One, N: Dim + DimName> One for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn one() -> Self {
Self::from_real(T::one())
}
#[inline]
fn is_one(&self) -> bool
where
Self: PartialEq,
{
self.real().is_one()
}
}
impl<T: Copy + Scalar + Num, N: Dim + DimName> Num for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
type FromStrRadixErr = <T as Num>::FromStrRadixErr;
#[inline]
fn from_str_radix(str: &str, radix: u32) -> Result<OHyperdual<T, N>, Self::FromStrRadixErr> {
<T as Num>::from_str_radix(str, radix).map(Self::from_real)
}
}
impl<T: Copy + Scalar + Float, N: Dim + DimName> NumCast for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
#[inline]
fn from<P: ToPrimitive>(n: P) -> Option<OHyperdual<T, N>> {
<T as NumCast>::from(n).map(Self::from_real)
}
}
macro_rules! impl_float_const {
($($c:ident),*) => {
$(
fn $c()-> Self { Self::from_real(T::$c()) }
)*
}
}
impl<T: Copy + Scalar + FloatConst + Zero, N: Dim + DimName> FloatConst for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
impl_float_const!(
E,
FRAC_1_PI,
FRAC_1_SQRT_2,
FRAC_2_PI,
FRAC_2_SQRT_PI,
FRAC_PI_2,
FRAC_PI_3,
FRAC_PI_4,
FRAC_PI_6,
FRAC_PI_8,
LN_10,
LN_2,
LOG10_E,
LOG2_E,
PI,
SQRT_2
);
}
macro_rules! impl_real_constant {
($($prop:ident),*) => {
$(
fn $prop() -> Self { Self::from_real(<T as Float>::$prop()) }
)*
}
}
macro_rules! impl_single_boolean_op {
($op:ident REAL) => {
fn $op(self) -> bool {
self.real().$op()
}
};
($op:ident OR) => {
fn $op(self) -> bool {
let mut b = self.real().$op();
for x in self.iter().skip(1) {
b |= x.$op();
}
b
}
};
($op:ident AND) => {
fn $op(self) -> bool {
let mut b = self.real().$op();
for x in self.iter().skip(1) {
b &= x.$op();
}
b
}
};
}
macro_rules! impl_boolean_op {
($($op:ident $t:tt),*) => {
$(impl_single_boolean_op!($op $t);)*
};
}
macro_rules! impl_real_op {
($($op:ident),*) => {
$(
fn $op(self) -> Self { Self::from_real(self.real().$op()) }
)*
}
}
impl<T: Copy + Scalar + Num + Zero, N: Dim + DimName> Sum for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
fn sum<I: Iterator<Item = OHyperdual<T, N>>>(iter: I) -> Self {
iter.fold(Self::zero(), |a, b| a + b)
}
}
impl<'a, T: Copy + Scalar + Num + Zero, N: Dim + DimName> Sum<&'a OHyperdual<T, N>> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
fn sum<I: Iterator<Item = &'a OHyperdual<T, N>>>(iter: I) -> Self {
iter.fold(Self::zero(), |a, b| a + *b)
}
}
impl<T: Copy + Scalar + Num + One, N: Dim + DimName> Product for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
fn product<I: Iterator<Item = OHyperdual<T, N>>>(iter: I) -> Self {
iter.fold(Self::one(), |a, b| a * b)
}
}
impl<'a, T: Copy + Scalar + Num + One, N: Dim + DimName> Product<&'a OHyperdual<T, N>> for OHyperdual<T, N>
where
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
fn product<I: Iterator<Item = &'a OHyperdual<T, N>>>(iter: I) -> Self {
iter.fold(Self::one(), |a, b| a * *b)
}
}
impl<T: Copy + Scalar, N: Dim + DimName> Float for OHyperdual<T, N>
where
T: Float + Signed + FloatConst,
DefaultAllocator: Allocator<N>,
Owned<T, N>: Copy,
{
impl_real_constant!(nan, infinity, neg_infinity, neg_zero, min_positive_value, epsilon, min_value, max_value);
impl_boolean_op!(
is_nan OR,
is_infinite OR,
is_finite AND,
is_normal AND,
is_sign_positive REAL,
is_sign_negative REAL
);
#[inline]
fn classify(self) -> FpCategory {
self.real().classify()
}
impl_real_op!(floor, ceil, round, trunc);
#[inline]
fn fract(self) -> Self {
let mut v = self;
v[0] = self.real().fract();
v
}
#[inline]
fn signum(self) -> Self {
Self::from_real(self.real().signum())
}
#[inline]
fn abs(self) -> Self {
let s = self.real().signum();
Self(self.map(|x| x * s))
}
#[inline]
fn max(self, other: Self) -> Self {
if self.real() > other.real() {
self
} else {
other
}
}
#[inline]
fn min(self, other: Self) -> Self {
if self.real() < other.real() {
self
} else {
other
}
}
#[inline]
fn abs_sub(self, rhs: Self) -> Self {
if self.real() > rhs.real() {
self.sub(rhs)
} else {
Self::zero()
}
}
#[inline]
fn mul_add(self, a: Self, b: Self) -> Self {
let mut dual = Self::from_real(self.real().mul_add(a.real(), b.real()));
for x in 1..self.len() {
dual[x] = self[x] * a.real() + self.real() * a[x] + b[x];
}
dual
}
#[inline]
fn recip(self) -> Self {
Self::one() / self
}
#[inline]
fn powi(self, n: i32) -> Self {
let r = self.real().powi(n - 1);
let nf = <T as NumCast>::from(n).expect("Invalid value") * r;
self.map_dual(self.real() * r, |x| nf * *x)
}
#[inline]
fn powf(self, n: Self) -> Self {
let real_part = self.real().powf(n.real());
let a = n.real() * self.real().powf(n.real() - T::one());
let b = real_part * self.real().ln();
let mut v = self.zip_map(&n, |x, y| a * x + b * y);
v[0] = real_part;
Self(v)
}
#[inline]
fn exp(self) -> Self {
let real = self.real().exp();
self.map_dual(real, |x| real * *x)
}
#[inline]
fn exp2(self) -> Self {
let real = self.real().exp2();
self.map_dual(real, |x| *x * T::LN_2() * real)
}
#[inline]
fn ln(self) -> Self {
self.map_dual(self.real().ln(), |x| *x / self.real())
}
#[inline]
fn log(self, base: Self) -> Self {
self.ln() / base.ln()
}
#[inline]
fn log2(self) -> Self {
self.map_dual(self.real().log2(), |x| *x / (self.real() * T::LN_2()))
}
#[inline]
fn log10(self) -> Self {
self.map_dual(self.real().log10(), |x| *x / (self.real() * T::LN_10()))
}
#[inline]
fn sqrt(self) -> Self {
let real = self.real().sqrt();
let d = T::from(1).unwrap() / (T::from(2).unwrap() * real);
self.map_dual(real, |x| *x * d)
}
#[inline]
fn cbrt(self) -> Self {
let real = self.real().cbrt();
self.map_dual(real, |x| *x / (T::from(3).unwrap() * real))
}
#[inline]
fn hypot(self, other: Self) -> Self {
let c = self.real().hypot(other.real());
let mut v = self.zip_map(&other, |x, y| (self.real() * x + other.real() * y) / c);
v[0] = c;
Self(v)
}
#[inline]
fn sin(self) -> Self {
let c = self.real().cos();
self.map_dual(self.real().sin(), |x| *x * c)
}
#[inline]
fn cos(self) -> Self {
let c = self.real().sin();
self.map_dual(self.real().cos(), |x| x.neg() * c)
}
#[inline]
fn tan(self) -> Self {
let t = self.real().tan();
let c = t * t + T::one();
self.map_dual(t, |x| *x * c)
}
#[inline]
fn asin(self) -> Self {
let c = (T::one() - self.real().powi(2)).sqrt();
self.map_dual(self.real().asin(), |x| *x / c)
}
#[inline]
fn acos(self) -> Self {
let c = (T::one() - self.real().powi(2)).sqrt();
self.map_dual(self.real().acos(), |x| x.neg() / c)
}
#[inline]
fn atan(self) -> Self {
let c = T::one() + self.real().powi(2);
self.map_dual(self.real().atan(), |x| *x / c)
}
#[inline]
fn atan2(self, other: Self) -> Self {
let c = self.real().powi(2) + other.real().powi(2);
let mut v = self.zip_map(&other, |x, y| (other.real() * x - self.real() * y) / c);
v[0] = self.real().atan2(other.real());
Self(v)
}
#[inline]
fn sin_cos(self) -> (Self, Self) {
let (s, c) = self.real().sin_cos();
let sn = self.map_dual(s, |x| *x * c);
let cn = self.map_dual(c, |x| x.neg() * s);
(sn, cn)
}
#[inline]
fn exp_m1(self) -> Self {
let c = self.real().exp();
self.map_dual(self.real().exp_m1(), |x| *x * c)
}
#[inline]
fn ln_1p(self) -> Self {
let c = self.real() + T::one();
self.map_dual(self.real().ln_1p(), |x| *x / c)
}
#[inline]
fn sinh(self) -> Self {
let c = self.real().cosh();
self.map_dual(self.real().sinh(), |x| *x * c)
}
#[inline]
fn cosh(self) -> Self {
let c = self.real().sinh();
self.map_dual(self.real().cosh(), |x| *x * c)
}
#[inline]
fn tanh(self) -> Self {
let real = self.real().tanh();
let c = T::one() - real.powi(2);
self.map_dual(real, |x| *x * c)
}
#[inline]
fn asinh(self) -> Self {
let c = (self.real().powi(2) + T::one()).sqrt();
self.map_dual(self.real().asinh(), |x| *x / c)
}
#[inline]
fn acosh(self) -> Self {
let c = (self.real() + T::one()).sqrt() * (self.real() - T::one()).sqrt();
self.map_dual(self.real().acosh(), |x| *x / c)
}
#[inline]
fn atanh(self) -> Self {
let c = T::one() - self.real().powi(2);
self.map_dual(self.real().atanh(), |x| *x / c)
}
#[inline]
fn integer_decode(self) -> (u64, i16, i8) {
self.real().integer_decode()
}
#[inline]
fn to_degrees(self) -> Self {
self.map_dual(self.real().to_degrees(), |x| x.to_degrees())
}
#[inline]
fn to_radians(self) -> Self {
self.map_dual(self.real().to_radians(), |x| x.to_radians())
}
}
pub type Dual<T> = Hyperdual<T, 2>;
impl<T: Copy + Scalar> Dual<T> {
#[inline]
pub fn new(real: T, dual: T) -> Dual<T> {
Dual::from_slice(&[real, dual])
}
#[inline]
pub fn dual(&self) -> T {
self[1]
}
#[inline]
pub fn dual_ref(&self) -> &T {
&self[1]
}
#[inline]
pub fn dual_mut(&mut self) -> &mut T {
&mut self[1]
}
}
pub type DualN<T, const N: usize> = Hyperdual<T, N>;
pub fn hyperspace_from_vector<T: Copy + Scalar + Num + Zero, const D: usize, const N: usize>(v: &SVector<T, N>) -> SVector<Hyperdual<T, D>, N> {
let mut space_slice = vec![Hyperdual::<T, D>::zero(); N];
for i in 0..N {
space_slice[i] = Hyperdual::<T, D>::from_fn(|j| {
if j == 0 {
v[i]
} else if i + 1 == j {
T::one()
} else {
T::zero()
}
});
}
SVector::<Hyperdual<T, D>, N>::from_row_slice(&space_slice)
}
pub fn vector_from_hyperspace<T: Scalar + Zero + Float, const DIM_VECTOR: usize, const DIM_HYPER: usize>(
x_dual: &SVector<DualN<T, DIM_HYPER>, { DIM_VECTOR }>,
) -> SVector<T, DIM_VECTOR> {
x_dual.map(|x| x.real())
}
pub fn hyperspace_from_ovector<T: Copy + Scalar + Num + Zero, D: Dim + DimName, N: Dim + DimName>(v: &OVector<T, N>) -> OVector<OHyperdual<T, D>, N>
where
DefaultAllocator: Allocator<D> + Allocator<N>,
Owned<T, D>: Copy,
{
let mut space_slice = vec![OHyperdual::<T, D>::zero(); N::dim()];
for i in 0..N::dim() {
space_slice[i] = OHyperdual::<T, D>::from_fn(|j| {
if j == 0 {
v[i]
} else if i + 1 == j {
T::one()
} else {
T::zero()
}
});
}
OVector::<OHyperdual<T, D>, N>::from_row_slice(&space_slice)
}
pub fn ovector_from_hyperspace<T: Scalar + Zero + Float, DimVector: Dim + DimName, DimHyper: Dim + DimName>(
x_dual: &OVector<OHyperdual<T, DimHyper>, DimVector>,
) -> OVector<T, DimVector>
where
DefaultAllocator: Allocator<DimVector> + Allocator<DimVector> + Allocator<DimHyper>,
<DefaultAllocator as Allocator<DimHyper>>::Buffer<T>: Copy,
{
x_dual.map(|x| x.real())
}