import math
import numpy as np
[docs]class Dual:
"""
A class to represent dual numbers, enabling basic operations
and automatic differentiation.
"""
# Initialising the 'Dual number' class
[docs] def __init__(self, real: float, dual: float):
"""
Initialize a dual number with its real and dual parts.
Parameters:
real (float): The real part of the dual number.
dual (float): The dual part of the dual number.
"""
# Preforming input validation using seperate class function to allow reuse
self.real = self._validate_input(real, "Real part")
self.dual = self._validate_input(dual, "Dual part")
# Use of '__repr__' to set a dual numbers output format 'Dual(real = ..., dual = ...)'
def __repr__(self):
"""
Output the dual number in a readable format.
Returns:
str: A string describing the dual number.
"""
return f"Dual(real={self.real}, dual={self.dual})"
# Overwrites addition operator for use on dual numbers
# Written to support addition of dual numbers together and dual numbers with scalar
[docs] def __add__(self, other):
"""
Define addition for Dual numbers and scalars.
This method supports addition in both standard and reverse cases:
- Dual + Dual
- Dual + scalar
- scalar + Dual (via __radd__)
Parameters:
other (Dual, int, or float): The value to add to the current Dual number.
Returns:
Dual: A new Dual number representing the result of the addition.
Raises:
TypeError: If 'other' is not a Dual, int, or float.
Examples:
1. Dual + Dual:
>>> x = Dual(2, 1)
>>> y = Dual(3, 2)
>>> x + y
Dual(real=5, dual=3)
2. Dual + scalar:
>>> x = Dual(2, 1)
>>> x + 3
Dual(real=5, dual=1)
3. Scalar + Dual:
>>> x = Dual(2, 1)
>>> 3 + x
Dual(real=5, dual=1)
4. Invalid input:
>>> x = Dual(2, 1)
>>> x + "string"
Traceback (most recent call last):
...
TypeError: Addition is only supported with Dual or scalar values.
"""
# If both operands are dual numbers add real and dual parts respectively
# (a+bε) + (c+dε) = (a+c) + (b+d)ε
if isinstance(other, Dual):
real_part = self.real + other.real
dual_part = self.dual + other.dual
return Dual(real_part, dual_part)
# If dual + int/float
# (a+bε) + c = (a+c) + bε
elif isinstance(other, (int, float)):
return Dual(self.real + other, self.dual)
# If not addition of dual and scalar - throw 'Type Error'
else:
raise TypeError("Addition failed as only supported with Dual or scalar values.")
# This overwrites the reverse addition opertor to support int/float + dual
# a + (c+dε) = (a+c) + dε
# Can be used as commutative
__radd__ = __add__
# Overwrites subtraction operator for use on dual numbers
# Written to support subtraction of dual numbers together and dual numbers with scalar
[docs] def __sub__(self, other):
"""
Define subtraction for Dual numbers and scalars.
This method supports subtraction in both standard and reverse cases:
- Dual - Dual
- Dual - scalar
- scalar - Dual (via __rsub__)
Parameters:
other (Dual, int, or float): The value to subtract from the current Dual number.
Returns:
Dual: A new Dual number representing the result of the subtraction.
Raises:
TypeError: If 'other' is not a Dual, int, or float.
Examples:
1. Dual - Dual:
>>> x = Dual(5, 3)
>>> y = Dual(2, 1)
>>> x - y
Dual(real=3, dual=2)
2. Dual - scalar:
>>> x = Dual(5, 3)
>>> x - 2
Dual(real=3, dual=3)
3. Scalar - Dual:
>>> x = Dual(5, 3)
>>> 10 - x
Dual(real=5, dual=-3)
4. Invalid input:
>>> x = Dual(5, 3)
>>> x - "string"
Traceback (most recent call last):
...
TypeError: Subtraction is only supported with Dual or scalar values.
"""
# If both operands are dual numbers subtract real and dual parts respectively
# (a+bε) - (c+dε) = (a-c) + (b-d)ε
if isinstance(other, Dual):
real_part = self.real - other.real
dual_part = self.dual - other.dual
return Dual(real_part, dual_part)
# If dual - int/float
# (a+bε) - c = (a-c) + bε
elif isinstance(other, (int, float)):
return Dual(self.real - other, self.dual)
# If not subtraction of dual and scalar - throw 'Type Error'
else:
raise TypeError("Subtraction failed as only supported with Dual or scalar values.")
# This overwrites the reverse subtraction opertor to support int/float - dual
# a - (c+dε) = (a-c) - dε
def __rsub__(self, other):
if isinstance(other, (int, float)):
return Dual(other - self.real, -self.dual)
else:
raise TypeError("Subtraction failed as only supported with Dual or scalar values.")
# Overwrites multiplication operator for use on dual numbers
# Written to support multiplication of dual numbers together and dual numbers with scalar
[docs] def __mul__(self, other):
"""
Define multiplication for Dual numbers and scalars.
This method supports multiplication in both standard and reverse cases:
- Dual * Dual
- Dual * scalar
- scalar * Dual (via __rmul__)
Parameters:
other (Dual, int, or float): The value to multiply with the current Dual number.
Returns:
Dual: A new Dual number representing the result of the multiplication.
Raises:
TypeError: If 'other' is not a Dual, int, or float.
Examples:
1. Dual * Dual:
>>> x = Dual(5, 3)
>>> y = Dual(2, 1)
>>> x * y
Dual(real=10, dual=11)
2. Dual * scalar:
>>> x = Dual(5, 3)
>>> x * 2
Dual(real=10, dual=6)
3. Scalar * Dual:
>>> x = Dual(5, 3)
>>> 2 * x
Dual(real=10, dual=6)
4. Invalid input:
>>> x = Dual(5, 3)
>>> x * "string"
Traceback (most recent call last):
...
TypeError: Multiplication is only supported with Dual or scalar values.
"""
# If both operands are dual numbers:
# (a + bε) * (c + dε) = ac + (ad + bc)ε
if isinstance(other, Dual):
real_part = self.real * other.real
dual_part = self.real * other.dual + self.dual * other.real
return Dual(real_part, dual_part)
# If dual * int/float
# (a + bε) * c = a*c + (b*c)ε
elif isinstance(other, (int, float)):
return Dual(self.real * other, self.dual * other)
# If not multiplication of dual and scalar - throw 'TypeError'
else:
raise TypeError("Multiplication failed as only supported with Dual or scalar values.")
# This overwrites the reverse multiplication opertor to support int/float * dual
# a - (c+dε) = (a-c) - dε
# Can be used as commutative
__rmul__ = __mul__
[docs] def __truediv__(self, other):
"""
Define division for Dual numbers and scalars.
This method supports division in both standard and reverse cases:
- Dual / Dual
- Dual / scalar
- scalar / Dual (via __rtruediv__)
Parameters:
other (Dual, int, or float): The value to divide by.
Returns:
Dual: A new Dual number representing the result of the division.
Raises:
TypeError: If 'other' is not a Dual, int, or float.
ZeroDivisionError: If dividing by zero (in the real part for Dual, or scalar zero).
Examples:
1. Dual / Dual:
>>> x = Dual(6, 4)
>>> y = Dual(2, 1)
>>> x / y
Dual(real=3.0, dual=0.5)
2. Dual / scalar:
>>> x = Dual(6, 4)
>>> x / 2
Dual(real=3.0, dual=2.0)
3. Scalar / Dual:
>>> x = Dual(6, 4)
>>> 12 / x
Dual(real=2.0, dual=-1.3333...)
4. Invalid input:
>>> x = Dual(6, 4)
>>> y = Dual(0, 1)
>>> x / y
Traceback (most recent call last):
...
ZeroDivisionError: Cannot divide by a Dual number with no real part.
"""
# If both operands are Dual numbers:
# (a + bε) / (c + dε) = (a/c) + ((b*c - a*d) / c^2)ε
if isinstance(other, Dual):
# If division by 0 (undefined) - throws ZeroDivisionError
if other.real == 0:
raise ZeroDivisionError("Cannot divide by a Dual number with no real part.")
real_part = self.real / other.real
dual_part = (self.dual * other.real - self.real * other.dual) / (other.real ** 2)
return Dual(real_part, dual_part)
# If dividing a Dual number by a scalar (int or float):
# (a + bε) / c = (a/c) + (b/c)ε
elif isinstance(other, (int, float)):
# If division by 0 (undefined) - throws ZeroDivisionError
if other == 0:
raise ZeroDivisionError("Cannot divide by zero.")
return Dual(self.real / other, self.dual / other)
# If not division of dual and scalar - throw 'TypeError'
else:
raise TypeError("Division is only supported with Dual or scalar values.")
# This overwrites the reverse multiplication opertor to support int/float / dual
# a / (c+dε) = (a / c) - (a * d / c^2)ε
def __rtruediv__(self, other):
if isinstance(other, (int, float)):
if self.real == 0:
raise ZeroDivisionError("Cannot divide by a Dual number with zero as its real part.")
real_part = other / self.real
dual_part = (-self.dual * other) / (self.real ** 2)
return Dual(real_part, dual_part)
else:
raise TypeError("Division is only supported with Dual or scalar values.")
# Overwrites the power operator for Dual numbers
def __pow__(self, n):
"""
Compute the power of the Dual number.
Parameters:
n (float or int): The exponent.
Returns:
Dual: The Dual number raised to the power `n`.
Raises:
TypeError: If `n` is not a float or int.
Examples:
>>> x = Dual(2, 1)
>>> x**2
Dual(real=4.0, dual=4.0)
"""
# If power is not an int or float - throw 'TypeError'
if not isinstance(n, (int, float)):
raise TypeError(f"Power must be an int or float.")
real_part = self.real**n
dual_part = n * self.real**(n - 1) * self.dual
return Dual(real_part, dual_part)
# Overwrites the equality operator for Dual numbers by comparing real and dual parts
# Allows tolerance of 1e-9 for floating point comparison
def __eq__(self, other):
"""
Check equality between two Dual objects.
Parameters:
other (Dual): The Dual number to compare against.
Returns:
bool: True if both the real and dual parts are equal, False otherwise.
Raises:
TypeError: If `other` is not an instance of Dual.
"""
# If not comparing Dual numbers - throw 'TypeError'
if not isinstance(other, Dual):
raise TypeError(f"Equality comparison is only supported between Dual objects, got {type(other).__name__}.")
return math.isclose(self.real, other.real, rel_tol=1e-9) and math.isclose(self.dual, other.dual, rel_tol=1e-9)
# Overwrites the inequality operator for Dual numbers by comparing real and dual parts
def __ne__(self, other):
"""
Check inequality between two Dual objects.
Parameters:
other (Dual): The Dual number to compare against.
Returns:
bool: True if either the real or dual parts are not equal, False otherwise.
Raises:
TypeError: If `other` is not an instance of Dual.
"""
# If not comparing Dual numbers - throw 'TypeError'
if not isinstance(other, Dual):
raise TypeError(f"Inequality comparison is only supported between Dual objects, got {type(other).__name__}.")
return not self.__eq__(other)
# Overwrites the 'less than' operator for Dual numbers by comparing real and dual parts
def __lt__(self, other):
"""
Check if this Dual object is less than another Dual object (based on real part).
Parameters:
other (Dual): The Dual number to compare against.
Returns:
bool: True if this Dual object is less than the other, False otherwise.
Raises:
TypeError: If `other` is not an instance of Dual.
"""
# If not comparing Dual numbers - throw 'TypeError'
if not isinstance(other, Dual):
raise TypeError(f"Less than comparison is only supported between Dual objects, got {type(other).__name__}.")
return self.real < other.real or (self.real == other.real and self.dual < other.dual)
# Overwrites the 'less than or equal to' operator for Dual numbers by comparing real and dual parts
def __le__(self, other):
"""
Check if this Dual object is less than or equal to another Dual object (based on real part).
Parameters:
other (Dual): The Dual number to compare against.
Returns:
bool: True if this Dual object is less than or equal to the other, False otherwise.
Raises:
TypeError: If `other` is not an instance of Dual.
"""
# If not comparing Dual numbers - throw 'TypeError'
if not isinstance(other, Dual):
raise TypeError(f"Less than or equal comparison is only supported between Dual objects, got {type(other).__name__}.")
return self.real < other.real or (self.real == other.real and self.dual <= other.dual)
# Overwrites the 'greater than' operator for Dual numbers by comparing real and dual parts
def __gt__(self, other):
"""
Check if this Dual object is greater than another Dual object (based on real part).
Parameters:
other (Dual): The Dual number to compare against.
Returns:
bool: True if this Dual object is greater than the other, False otherwise.
Raises:
TypeError: If `other` is not an instance of Dual.
"""
# If not comparing Dual numbers - throw 'TypeError'
if not isinstance(other, Dual):
raise TypeError(f"Greater than comparison is only supported between Dual objects, got {type(other).__name__}.")
return self.real > other.real or (self.real == other.real and self.dual > other.dual)
# Overwrites the 'greater than or equal to' operator for Dual numbers by comparing real and dual parts
def __ge__(self, other):
"""
Check if this Dual object is greater than or equal to another Dual object (based on real part).
Parameters:
other (Dual): The Dual number to compare against.
Returns:
bool: True if this Dual object is greater than or equal to the other, False otherwise.
Raises:
TypeError: If `other` is not an instance of Dual.
"""
# If not comparing Dual numbers - throw 'TypeError'
if not isinstance(other, Dual):
raise TypeError(f"Greater than or equal comparison is only supported between Dual objects, got {type(other).__name__}.")
return self.real > other.real or (self.real == other.real and self.dual >= other.dual)
# Defines a generic function implementing the preperty of dual numbers
# f(a+bε) = f(a) = f'(a)*bε
# Used internally by specific functions which input function and its derivative
[docs] def _dual_function(self, func, func_deriv):
"""
Internal method to apply a function to a dual number.
Supported Functions:
- Sine: .sin()
- Cosine: .cos()
- Tangent: .tan()
- Arcsin: .arcsin()
- Arccos: .arccos()
- Arctan: .arctan()
- Sinh: .sinh()
- Cosh: .cosh()
- Tanh: .tanh()
- Exponential: .exp()
- Logarithm: .log()
- Powers: .pow(n)
- Square Root: .sqrt()
Parameters:
func (callable): The mathematical function applied to dual number.
func_deriv (callable): The derivative of the function.
Returns:
Dual: The out the function on a dual number
"""
real_part = func(self.real)
dual_part = func_deriv(self.real) * self.dual
return Dual(real_part, dual_part)
# Defines Sin() operator using _dual_function
# Derivative: Cos()
[docs] def sin(self):
"""
Compute the sine of the Dual number.
Returns:
Dual: Sine of original dual number
Examples:
>>> x = Dual(2, 1)
>>> x.sin()
Dual(real=0.9092..., dual=-0.4161...)
"""
return self._dual_function(math.sin, math.cos)
# Defines Cos() operator using _dual_function
# Derivative: -Sin()
[docs] def cos(self):
"""
Compute the cosine of the Dual number.
Returns:
Dual: Cosine of original dual number
Examples:
>>> x = Dual(2, 1)
>>> x.cos()
Dual(real=-0.4161..., dual=-0.9093...)
"""
return self._dual_function(math.cos, lambda x: -math.sin(x))
# Defines tan() operator using _dual_function - log(x)
# Derivative: 1/cos^2(x)
[docs] def tan(self):
"""
Compute the tangent of the Dual number.
Returns:
Dual: Tangent of original dual number
Raises:
ValueError: The tangent function is undefined as cosine of real part equals 0
Examples:
>>> x = Dual(2, 1)
>>> x.tan()
Dual(real=-2.1850..., dual=5.7744...)
>>> Dual(math.pi / 2, 1).tan()
Traceback (most recent call last):
...
ValueError: Tangent undefined when cosine of real part equals 0.
"""
# Checks for regions in which tan is undefined - eg. pi/2
# Same regions in which cos(x) = 0 - uses tolerance of 1e-9 around value
# If case: throws 'ValueError'
if math.isclose(math.cos(self.real), 0, abs_tol=1e-9):
raise ValueError("Tangent undefined when cosine of real part equals 0.")
return self._dual_function(math.tan, lambda x: 1 / (math.cos(x) ** 2))
# Defines exp() operator using _dual_function - e^x
# Derivative: exp()
[docs] def exp(self):
"""
Compute the exponential of the Dual number.
Returns:
Dual: Expontential of original dual number
Examples:
>>> x = Dual(2, 1)
>>> x.exp()
Dual(real=7.3891..., dual=7.3891...)
"""
return self._dual_function(math.exp, math.exp)
# Defines log() operator using _dual_function - log(x)
# Derivative: 1/x
[docs] def log(self):
"""
Compute the natural logarithm of the Dual number.
Returns:
Dual: Natural logirithim of original dual number
Raises:
ValueError: If the real part of dual number is non-positive
Examples:
>>> x = Dual(2, 1)
>>> x.log()
Dual(real=0.6931..., dual=0.5)
>>> x = Dual(0, 1)
>>> x.log()
Traceback (most recent call last):
...
ValueError: Logarithm is undefined for dual numbers with non-positive real parts.
"""
# Checks for Logirithm of non positive number for which undefined
# If so throws 'Value Error'
if self.real <= 0:
raise ValueError("Logarithm is undefined for dual numbers with non-positive real parts.")
return self._dual_function(math.log, lambda x: 1 / x)
# Defines arcsin() operator using _dual_function
# Derivative: 1/sqrt(1-x^2)
[docs] def arcsin(self):
"""
Compute the arcsine (inverse sine) of the Dual number.
Returns:
Dual: Arcsine of the original Dual number.
Raises:
ValueError: If the real part is outside the range [-1, 1].
Examples:
>>> x = Dual(0.5, 1)
>>> x.arcsin()
Dual(real=0.5236..., dual=1.1547...)
>>> x = Dual(1.5, 1)
>>> x.arcsin()
Traceback (most recent call last):
...
ValueError: Arcsine is only defined for real parts in the range [-1, 1].
"""
# Checks for defined inputs where arcsin is defined: [-1,1]
# If so throws 'Value Error'
if not -1 <= self.real <= 1:
raise ValueError("Arcsine is only defined for real parts in the range [-1, 1].")
return self._dual_function(math.asin, lambda x: 1 / math.sqrt(1 - x**2))
# Defines arccos() operator using _dual_function
# Derivative: -1/sqrt(1-x^2)
[docs] def arccos(self):
"""
Compute the arccosine (inverse cosine) of the Dual number.
Returns:
Dual: Arccosine of the original Dual number.
Raises:
ValueError: If the real part is outside the range [-1, 1].
Examples:
>>> x = Dual(0.5, 1)
>>> x.arccos()
Dual(real=1.0472..., dual=-1.1547...)
>>> x = Dual(-1.5, 1)
>>> x.arccos()
Traceback (most recent call last):
...
ValueError: Arccosine is only defined for real parts in the range [-1, 1].
"""
# Checks for defined inputs where arcccos is defined: [-1,1]
# If so throws 'Value Error'
if not -1 <= self.real <= 1:
raise ValueError("Arccosine is only defined for real parts in the range [-1, 1].")
return self._dual_function(math.acos, lambda x: -1 / math.sqrt(1 - x**2))
# Defines arctan() operator using _dual_function
# Derivative: 1/sqrt(1+x^2)
[docs] def arctan(self):
"""
Compute the arctangent (inverse tangent) of the Dual number.
Returns:
Dual: Arctangent of the original Dual number.
Examples:
>>> x = Dual(1, 1)
>>> x.arctan()
Dual(real=0.7854..., dual=0.5)
"""
return self._dual_function(math.atan, lambda x: 1 / (1 + x**2))
# Defines hyperbolic operator: sinh() using _dual_function
# Derivative: hyperbolic operator: cosh()
[docs] def sinh(self):
"""
Compute the hyperbolic sine of the Dual number.
Returns:
Dual: Hyperbolic sine of the original Dual number.
Examples:
>>> x = Dual(1, 1)
>>> x.sinh()
Dual(real=1.1752..., dual=1.5431...)
"""
return self._dual_function(math.sinh, math.cosh)
# Defines hyperbolic operator: cosh() using _dual_function
# Derivative: hyperbolic operator: sinh()
[docs] def cosh(self):
"""
Compute the hyperbolic cosine of the Dual number.
Returns:
Dual: Hyperbolic cosine of the original Dual number.
Examples:
>>> x = Dual(1, 1)
>>> x.cosh()
Dual(real=1.5431..., dual=1.1752...)
"""
return self._dual_function(math.cosh, math.sinh)
# Defines hyperbolic operator: tanh() using _dual_function
# Derivative: hyperbolic operator: 1-tanh^2()
[docs] def tanh(self):
"""
Compute the hyperbolic tangent of the Dual number.
Returns:
Dual: Hyperbolic tangent of the original Dual number.
Examples:
>>> x = Dual(1, 1)
>>> x.tanh()
Dual(real=0.7616..., dual=0.4200...)
"""
return self._dual_function(math.tanh, lambda x: 1 - math.tanh(x)**2)
# Defines square root operator - sqrt() using _dual_function
# Derivative: (1/2) x^(-1/2)
[docs] def sqrt(self):
"""
Compute the square root of the Dual number.
Returns:
Dual: Square root of the original dual number.
Raises:
ValueError: If the real part of dual number is negative.
Examples:
>>> x = Dual(4, 1)
>>> x.sqrt()
Dual(real=2.0, dual=0.25)
>>> x = Dual(0, 1)
>>> x.sqrt()
Traceback (most recent call last):
...
ValueError: Square root is undefined when real parts of dual number are negetive.
"""
if self.real < 0:
raise ValueError("Square root is undefined when real parts of dual number are negetive.")
return self._dual_function(math.sqrt, lambda x: 0.5 / math.sqrt(x))
# Defines 'to the power of' operator - pow(n) using _dual_function
# Derivative: n * x^(n-1)
[docs] def pow(self, n):
"""
Compute the Dual number raised to a power n.
Parameters:
n (float or int): The power to which the Dual number is raised.
Returns:
Dual: The original dual number raised to the given power.
Raises:
TypeError: If n is not an int or float.
Examples:
>>> x = Dual(2, 1)
>>> x.pow(3)
Dual(real=8.0, dual=12.0)
"""
# Validate the input type of n
if not isinstance(n, (int, float)):
raise TypeError(f"Power must be an int or float, got {type(n).__name__}.")
return self._dual_function(lambda x: x**n, lambda x: n * x**(n-1))
def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
"""
Handle NumPy ufuncs (universal functions) for Dual objects.
Supported unfuncs:
- Addition: np.add()
- Subtraction: np.subtract()
- Multiplication: np.multiply()
- Division: np.divide()
- Sine: np.sin()
- Cosine: np.cos()
- Tangent: np.tan()
- Arcsin: np.arcsin()
- Arccos: np.arccos()
- Arctan: np.arctan()
- Sinh: np.sinh()
- Cosh: np.cosh()
- Tanh: np.tanh()
- Exponential: np.exp()
- Logarithm: np.log()
- Powers: np.pow(n)
- Square Root: np.sqrt()
- Equal: np.equal()
- Not Equal: np.not_equal()
- Less: np.less()
- Less Equal: np.less_equal()
- Greater: np.greater()
- Greater Equal: np.greater_equal()
Returns:
Dual or NotImplemented: The result of the ufunc operation, or NotImplemented if the operation is not supported.
Example:
>>> x = Dual(2, 1)
>>> y = Dual(3, 2)
>>> np.add(x, y)
Dual(real=5.0, dual=3.0)
>>> np.multiply(x, y)
Dual(real=6.0, dual=7.0)
"""
# Ensure the method is '__call__' (standard operation)
print(f"Entered __array_ufunc__ for ufunc: {ufunc}")
if method != '__call__':
return NotImplemented
# Extract real and dual parts from inputs
real_parts = []
dual_parts = []
for x in inputs:
if isinstance(x, Dual):
real_parts.append(x.real)
dual_parts.append(x.dual)
else: # Treat scalars or standard arrays as real numbers with dual=0
real_parts.append(x)
dual_parts.append(0)
# Convert lists to arrays for broadcasting
real_parts = np.array(real_parts)
dual_parts = np.array(dual_parts)
# Now define the operations for each option of ufunc
# add and sub are element-wise operations and therefore can keep unfunc as operation to support vectorization
if ufunc in {np.add, np.subtract}:
real_result = ufunc(*real_parts, **kwargs)
dual_result = ufunc(*dual_parts, **kwargs)
# Multiplication
elif ufunc == np.multiply:
# Real part is the product of the real parts
result_real = np.prod(real_parts, axis=0)
# Calculates the dual parts (ie first order coefficients) of the product by summing the products of the real parts excluding the current index
# non vectorized version below for reference
# result_dual = 0 # Initialize the result for the dual part
# # Loop through each dual component
# for i, dual in enumerate(duals):
# # Compute the product of the real parts excluding the current index
# product_of_reals = 1 # Start with a multiplicative identity
# for j, r in enumerate(reals):
# if j != i: # Exclude the current index
# product_of_reals *= r
# # Add the contribution of the current dual to the total dual result
# result_dual += dual * product_of_reals
# Vectorized version exploiting numpy's efficeint broadcasting
result_dual = sum( (dual * np.prod([r for j, r in enumerate(real_parts) if j != i],axis=0)) for i, dual in enumerate(dual_parts))
# Division
# As np.divide does not accept variafic arguments, only two inputs should be supported and vectorisation is not required
elif ufunc == np.divide: # Division
# Ensures only two inputs for division
if len(inputs) != 2:
raise ValueError("np.divide only supports two inputs (numerator and denominator).")
# Real parts divided by respectively
result_real = real_parts[0] / real_parts[1]
# Dual parts calculated as before
result_dual = (dual_parts[0] * real_parts[1] - real_parts[0] * dual_parts[1]) / (real_parts[1] ** 2)
# Trigonometric functions
# Done part wise and using numpy's vectorization
# Use same logic as _dual_function
#sin
if ufunc == np.sin:
real_result = np.sin(real_parts)
dual_result = np.cos(real_parts) * dual_parts
#cos
elif ufunc == np.cos:
real_result = np.cos(real_parts)
dual_result = -np.sin(real_parts) * dual_parts
#tan
elif ufunc == np.tan:
real_result = np.tan(real_parts)
dual_result = (1 / np.cos(real_parts) ** 2) * dual_parts
# Inverse trigonometric functions
# Done part wise and using numpy's vectorization
# Use same logic as _dual_function
#arcsin
elif ufunc == np.arcsin:
real_result = np.arcsin(real_parts)
dual_result = dual_parts / np.sqrt(1 - real_parts ** 2)
#arccos
elif ufunc == np.arccos:
real_result = np.arccos(real_parts)
dual_result = -dual_parts / np.sqrt(1 - real_parts ** 2)
#arctan
elif ufunc == np.arctan:
real_result = np.arctan(real_parts)
dual_result = dual_parts / (1 + real_parts ** 2)
# Hyperbolic functions
# Done part wise and using numpy's vectorization
# Use same logic as _dual_function
# sinh
elif ufunc == np.sinh:
real_result = np.sinh(real_parts)
dual_result = np.cosh(real_parts) * dual_parts
# cosh
elif ufunc == np.cosh:
real_result = np.cosh(real_parts)
dual_result = np.sinh(real_parts) * dual_parts
# tanh
elif ufunc == np.tanh:
real_result = np.tanh(real_parts)
dual_result = (1 - np.tanh(real_parts) ** 2) * dual_parts
# Exponential, Logarithmic and Square root functions
# Done part wise and using numpy's vectorization
# Use same logic as _dual_function
# exp
elif ufunc == np.exp:
real_result = np.exp(real_parts)
dual_result = np.exp(real_parts) * dual_parts
# log
elif ufunc == np.log:
real_result = np.log(real_parts)
dual_result = dual_parts / real_parts
# Square Root
elif ufunc == np.sqrt:
real_result = np.sqrt(real_parts)
dual_result = 0.5 * dual_parts / np.sqrt(real_parts)
# Power
elif ufunc == np.power: # Power
if len(inputs) != 2:
raise ValueError("np.power only supports two inputs (base and exponent).")
if not isinstance(inputs[1], (int, float)):
raise ValueError(f"Exponent must be an int/float for the use of np.power with Dual class, got type {type(inputs[1]).__name__}.")
real_result = np.power(real_parts[0], real_parts[1])
dual_result = real_parts[1] * np.power(real_parts[0], real_parts[1] - 1) * dual_parts[0]
# Comparison operators
# Apart from (not) equal to are done on real parts only
elif ufunc in {np.less, np.less_equal, np.greater, np.greater_equal}:
return ufunc(*real_parts, **kwargs)
# not equal to and equal to are done on both real and dual parts
elif ufunc in {np.equal, np.not_equal}:
return ufunc(*real_parts, **kwargs) and ufunc(*dual_parts, **kwargs)
else:
raise NotImplementedError(f"Ufunc {ufunc} not implemented for Dual numbers")
# Return final result of all operations
#
print(real_result)
print(dual_result)
return np.array([Dual(r, d) for r, d in zip(real_result, dual_result)])