## \file arithmetic.py
# \brief This file contains several functions and classes that are used
# for numeric computations in the other modules of this library.
# \author <a href="http://thomas.reidemeister.org/" target="_blank">
# Thomas Reidemeister</a>
## \namespace scuq::arithmetic
# \brief This namespace contains several functions and classes that are used
# for numeric computations within this class library.
## \defgroup arithmetic The Arithmetic Module
#
# This module contains several functions, classes, and constants
# that are used for numeric computations in the other modules of
# this library.
# \author <a href="http://thomas.reidemeister.org/" target="_blank">
# Thomas Reidemeister</a>
# \addtogroup arithmetic
# @{
# standard module
import operator
import numbers
import numpy
import warnings
def _require_int(value, name):
if not isinstance(value, int):
raise TypeError("%s must be an int" % name)
def _require_number(value, name="value"):
if not isinstance(value, (RationalNumber, numbers.Number)):
raise TypeError("%s must be a RationalNumber or number" % name)
def _deprecated_python2_method(name, replacement):
warnings.warn(
"%s is deprecated and kept only for Python 2 compatibility; use %s instead"
% (name, replacement),
DeprecationWarning,
stacklevel=2,
)
[docs]
def Coerce(x, y):
try:
t = type(x + y)
return (t(x), t(y))
except TypeError:
raise NotImplementedError
[docs]
def gcd(m, n):
"""Calculate the greatest common divisor.
:param n: First integer value (greater or equal to zero).
:param m: Second value (greater or equal to zero).
:returns: The greatest common divisor of the inputs.
"""
_require_int(m, "m")
_require_int(n, "n")
if n < 0 or m < 0:
raise ValueError("gcd arguments must be greater than or equal to zero")
m = int(m)
n = int(n)
if n == 0:
return m
else:
return gcd(n, m % n)
[docs]
def rational(n, d):
"""This function provides an interface for rational numbers
creation, as suggested in PEP 239.
:param d: The denominator (must be an interger type).
:param n: The nominator (must be an integer type).
:returns: An instance of RationalNumber
"""
return RationalNumber(n, d)
[docs]
def complex_to_matrix(c):
"""Convert a complex number to its 2x2 NumPy array representation."""
val = complex(c)
return numpy.array([[val.real, -val.imag], [val.imag, val.real]])
[docs]
class RationalNumber:
"""Represent a rational number with integer dividend and divisor.
The arithmetic operators emulate rational arithmetic where possible.
For unsupported mixed-type operations, methods may fall back to
floating-point behavior.
"""
[docs]
def __init__(self, dividend, divisor=1):
"""Default constructor.
This initializes the rational number.
:param dividend: An integer representing the dividend of this
rational number.
:param divisor: An integer representing the divisor of this
rational number. If this parameter is obmitted it is
initialized to 1.
"""
_require_int(dividend, "dividend")
_require_int(divisor, "divisor")
if divisor == 0:
raise ArithmeticError("Divide by zero")
self._divisor = int(divisor)
self._dividend = int(dividend)
self.normalize()
[docs]
def __str__(self):
"""This method returns a string representing this rational number.
:returns: A string representing this rational number.
"""
if self._divisor == 1:
return str(self._dividend)
else:
return "(" + str(self._dividend) + "/" + str(self._divisor) + ")"
[docs]
def normalize(self):
"""This method maintains the canonical form of this rational
number and avoids negative divisors.
"""
if self._divisor < 0:
self._dividend = -self._dividend
self._divisor = -self._divisor
mygcd = gcd(abs(self._dividend), self._divisor)
self._dividend = self._dividend // mygcd
self._divisor = self._divisor // mygcd
### The following methods are used to emulate the
### numeric behaviour.
[docs]
def __complex__(self):
"""Cast this rational number to a complex number.
:returns: A complex number, having a zero imaginary part.
"""
return complex(self.__float__())
[docs]
def __int__(self):
"""Cast this rational number to an integer.
integer overflow.
:returns: An integer.
"""
return int(self._dividend) // int(self._divisor)
[docs]
def __float__(self):
"""Cast this rational number to a floating point number.
:returns: An integer.
"""
return float(operator.truediv(self._dividend, self._divisor))
[docs]
def __add__(self, value):
"""Add a number and return the result.
:param value: The number to add.
:returns: The sum of this instance and the argument.
"""
if not isinstance(value, (RationalNumber, numbers.Number)):
return value.__radd__(self)
if isinstance(value, RationalNumber):
selfDividend = self._dividend * value._divisor
otherDividend = value._dividend * self._divisor
newDivisor = self._divisor * value._divisor
return RationalNumber(selfDividend + otherDividend, newDivisor)
# elif isinstance(value, ucomponents.UncertainComponent):
# return value.__add__(self)
elif isinstance(value, int):
return RationalNumber(self._dividend + value * self._divisor, self._divisor)
else:
return float(self) + value
[docs]
def __sub__(self, value):
"""Substract a number and return the result.
:param value: The number to substract.
:returns: The difference of this instance and the argument.
"""
if not isinstance(value, (RationalNumber, numbers.Number)):
return value.__rsub__(self)
if isinstance(value, RationalNumber):
selfDividend = self._dividend * value._divisor
otherDividend = value._dividend * self._divisor
newDivisor = self._divisor * value._divisor
return RationalNumber(selfDividend - otherDividend, newDivisor)
# elif isinstance(value, ucomponents.UncertainComponent): #self-value = -value+self
# v=-value
# return v.__add__(self)
elif isinstance(value, int):
return RationalNumber(self._dividend - value * self._divisor, self._divisor)
else:
return float(self) - value
[docs]
def __mul__(self, value):
"""Multiply a number and return the result.
:param value: The number to multiply.
:returns: The product of this instance and the argument.
"""
if not isinstance(value, (RationalNumber, numbers.Number)):
return value.__rmul__(self)
if isinstance(value, RationalNumber):
newDividend = self._dividend * value._dividend
newDivisor = self._divisor * value._divisor
return RationalNumber(newDividend, newDivisor)
elif isinstance(value, int):
return RationalNumber(self._dividend * value, self._divisor)
else:
return float(self) * value
[docs]
def __pow__(self, value):
"""Raise this rational number to the given power and return
the result.
a floating point number will be returned. Note that
this may result in a loss of accuracy.
:param value: A numeric value representing the power.
:returns: A new rational number representing power of this instance.
"""
if not isinstance(value, (RationalNumber, numbers.Number)):
return value.__rpow__(self)
if not isinstance(value, float) and (
isinstance(value, int) or value.is_integer()
):
if value < 0:
return pow(~self, -value)
else:
value = int(value)
return RationalNumber(self._dividend**value, self._divisor**value)
## something else (maybe float)
return float(self) ** value
[docs]
def __rpow__(self, value):
"""Raise another value to the power of this rational number.
the result.
:param value: A value to be raised to the power.
:returns: A new rational number representing power of this instance.
"""
_require_number(value)
if self.is_integer():
return value ** int(self)
else:
return value ** float(self)
[docs]
def __truediv__(self, value):
"""Divide by another number and return the result.
:param value: A number.
:returns: The fraction of this instance and the number.
"""
if not isinstance(value, (RationalNumber, numbers.Number)):
return value.__rtruediv__(self)
if isinstance(value, RationalNumber):
if value._dividend == 0:
raise ArithmeticError("Divide by zero")
newDividend = self._dividend * value._divisor
newDivisor = self._divisor * value._dividend
return RationalNumber(newDividend, newDivisor)
elif isinstance(value, int):
if value == 0:
raise ArithmeticError("Divide by zero")
newDividend = self._dividend
newDivisor = self._divisor * value
return RationalNumber(newDividend, newDivisor)
else:
return float(self) / value
[docs]
def __neg__(self):
"""This method returns the negative of this instance.
:returns: A new rational number.
"""
return RationalNumber(-self._dividend, self._divisor)
[docs]
def __pos__(self):
"""This method returns a copy of this instance.
:returns: A new rational number.
"""
return RationalNumber(self._dividend, self._divisor)
[docs]
def __abs__(self):
"""This method returns the absolute value of this instance.
:returns: A new rational number.
"""
if self._dividend < 0:
return self.__neg__()
return self.__pos__()
[docs]
def __invert__(self):
"""This method returns a new rational number that
swapped dividend and divsor of this instance.
:returns: A new rational number.
"""
return RationalNumber(self._divisor, self._dividend)
[docs]
def get_dividend(self):
"""Returns the dividend of this instance.
:returns: The dividend of this instance.
"""
return self._dividend
[docs]
def get_divisor(self):
"""Returns the divisor of this instance.
:returns: The divisor of this instance.
"""
return self._divisor
[docs]
def __eq__(self, value):
"""Checks if this instance is equal to a number.
:param value: The value to compare to.
:returns: If this rational number is equal to the argument.
"""
if isinstance(value, RationalNumber):
return self._divisor == value._divisor and self._dividend == value._dividend
elif isinstance(value, int):
return self._divisor == 1 and self._dividend == value
elif isinstance(value, numbers.Number):
return float(self) == value
else:
return False
[docs]
def __lt__(self, value):
"""Checks if this instance is less than another number.
:param value: The value to compare to.
:returns: True, if this rational number is less than the argument.
"""
_require_number(value)
if isinstance(value, int):
return self._dividend < self._divisor * int(value)
if isinstance(value, RationalNumber):
return self._dividend * value._divisor < value._dividend * self._divisor
# something else
return float(self) < value
[docs]
def __ne__(self, value):
"""Checks if this instance unequal to another number.
:param value: The value to compare to.
:returns: True, if this rational number unequal to the argument.
"""
if value == None:
return True
_require_number(value)
if isinstance(value, int):
return self._divisor != 1 or self._dividend != value
if isinstance(value, RationalNumber):
return self._divisor != value._divisor or self._dividend != value._dividend
# something else
return float(self) != value
[docs]
def __gt__(self, value):
"""Checks if this instance is greater than another number.
:param value: The value to compare to.
:returns: True, if this rational number is greater than the argument.
"""
_require_number(value)
if isinstance(value, int):
return self._dividend > self._divisor * value
if isinstance(value, RationalNumber):
return self._dividend * value._divisor > value._dividend * self._divisor
# something else
return float(self) > value
[docs]
def __ge__(self, value):
"""Checks if this instance is greater or equal to another number.
:param value: The value to compare to.
:returns: True, if this rational number is greater or equal to the argument.
"""
return self.__gt__(value) or self.__eq__(value)
[docs]
def __le__(self, value):
"""Checks if this instance is less or equal to another number.
:param value: The value to compare to.
:returns: True, if this rational number is less or equal to the argument.
"""
return self.__lt__(value) or self.__eq__(value)
[docs]
def __cmp__(self, value):
"""Compares this instance to another number.
:param value: The value to compare to.
:returns: -1: if this instance is less...; +1: if this is greater than the argument; 0 otherwise
"""
_deprecated_python2_method("__cmp__", "the rich comparison operators")
if not isinstance(value, numbers.Number):
raise TypeError("value must be a number")
if self.__lt__(value):
return -1
if self.__gt__(value):
return +1
return 0
[docs]
def __bool__(self):
"""Check if this instance is nonzero.
:returns: True, if the dividend is nonzero.
"""
return self._dividend != 0
### The same arithmetic Operations again, now for
### left arguments.
### Attention: only +,-,*,/ are defined in this case.
[docs]
def __radd__(self, value):
"""Left addition of a numeric value.
:param value: A value to left from this instance.
"""
# since this operation is symmetric
if isinstance(value, (RationalNumber, int)):
return self.__add__(value)
else:
return float(self) + value
[docs]
def __rsub__(self, value):
"""Right substraction of a numeric value.
:param value: A value to left from this instance.
"""
if isinstance(value, (RationalNumber, int)):
return (-self).__add__(value)
else:
return value - float(self)
[docs]
def __rmul__(self, value):
"""Right multiplication of a numeric value.
:param value: A value to left from this instance.
"""
if isinstance(value, (RationalNumber, int)):
return self.__mul__(value)
else:
return float(self) * value
[docs]
def __rtruediv__(self, value):
"""Right division of a numeric value.
:param value: A value to left from this instance.
"""
return (~self).__mul__(value)
[docs]
def __getstate__(self):
"""Serialization using pickle.
:returns: A string that represents the serialized form of this instance.
"""
return (self._dividend, self._divisor)
[docs]
def __setstate__(self, state):
"""Deserialization using pickle.
:param state: The state of the object.
"""
self._dividend, self._divisor = state
[docs]
def value_of(number):
"""Return a ``RationalNumber`` for an integer-like input.
:raises TypeError: If ``number`` is neither integer-like nor a
``RationalNumber``.
"""
# no conversion needed
if isinstance(number, RationalNumber):
return number
if isinstance(number, int) or (
isinstance(number, float) and number.is_integer()
):
return RationalNumber(int(number))
raise TypeError("Illegal Argument")
value_of = staticmethod(value_of)
[docs]
def is_integer(self):
"""Check wether this instance could be
be casted to long accurately.
:returns: True, if the divisor is equal to one.
"""
return self._divisor == 1
### The definition of numpy ufuncts
### All of these methods cast the rational numbers
### to float
[docs]
def arccos(self):
"""This method provides the broadcast interface for
numpy.arccos.
:returns: The inverse Cosine of this number.
"""
return numpy.arccos(float(self))
[docs]
def arccosh(self):
"""This method provides the broadcast interface for
numpy.arccosh.
:returns: The inverse hyperbolic Cosine of this number.
"""
return numpy.arccosh(float(self))
[docs]
def arcsin(self):
"""This method provides the broadcast interface for
numpy.arcsin.
:returns: The inverse Sine of this number.
"""
return numpy.arcsin(float(self))
[docs]
def arcsinh(self):
"""This method provides the broadcast interface for
numpy.arcsinh.
:returns: The inverse hyperbolic Sine of this number.
"""
return numpy.arcsinh(float(self))
[docs]
def arctan(self):
"""This method provides the broadcast interface for
numpy.arctan.
:returns: The inverse Tangent of this number.
"""
return numpy.arctan(float(self))
[docs]
def arctanh(self):
"""This method provides the broadcast interface for
numpy.arctanh.
:returns: The inverse hyperbolic Tangent of this number.
"""
return numpy.arctanh(float(self))
[docs]
def cos(self):
"""This method provides the broadcast interface for
numpy.cos.
:returns: The Cosine of this number.
"""
return numpy.cos(float(self))
[docs]
def cosh(self):
"""This method provides the broadcast interface for
numpy.cosh.
:returns: The hyperbolic Cosine of this number.
"""
return numpy.cosh(float(self))
[docs]
def tan(self):
"""This method provides the broadcast interface for
numpy.tan.
:returns: The Tangent of this number.
"""
return numpy.tan(float(self))
[docs]
def tanh(self):
"""This method provides the broadcast interface for
numpy.tanh.
:returns: The hyperbolic Tangent of this number.
"""
return numpy.tanh(float(self))
[docs]
def log10(self):
"""This method provides the broadcast interface for
numpy.log10.
:returns: The decadic Logarithm of this number.
"""
return numpy.log10(float(self))
[docs]
def sin(self):
"""This method provides the broadcast interface for
numpy.sin.
:returns: The Sine of this number.
"""
return numpy.sin(float(self))
[docs]
def sinh(self):
"""This method provides the broadcast interface for
numpy.sinh.
:returns: The hyperbolic Sine of this number.
"""
return numpy.sinh(float(self))
[docs]
def sqrt(self):
"""This method provides the broadcast interface for
numpy.sqrt.
:returns: The Square Root of this number.
"""
return numpy.sqrt(float(self))
[docs]
def fabs(self):
"""This method provides the broadcast interface for
numpy.fabs.
:returns: The absolute value of this number.
"""
return numpy.fabs(float(self))
[docs]
def absolute(self):
"""This method provides the broadcast interface for
numpy.absolute.
:returns: The absolute value of this number.
"""
return numpy.absolute(float(self))
[docs]
def floor(self):
"""This method provides the broadcast interface for
numpy.floor.
:returns: The largest integer less than or equal to this number.
"""
return numpy.floor(float(self))
[docs]
def ceil(self):
"""This method provides the broadcast interface for
numpy.ceil.
:returns: The largest integer greater than or equal to this number.
"""
return numpy.ceil(float(self))
[docs]
def exp(self):
"""This method provides the broadcast interface for
numpy.exp.
:returns: The Exponential of this number.
"""
return numpy.exp(float(self))
[docs]
def log(self):
"""This method provides the broadcast interface for
numpy.log.
:returns: The Natural Logarithm of this number.
"""
return numpy.log(float(self))
[docs]
def log2(self):
"""This method provides the broadcast interface for
numpy.log2.
:returns: The binary logarithm of this number.
"""
return numpy.log2(float(self))
[docs]
def square(self):
"""This method provides the broadcast interface for
numpy.square.
:returns: The binary logarithm of this number.
"""
return self * self
[docs]
def arctan2(self, other):
"""This method provides the broadcast interface for
numpy.arctan2.
:param other: Another rational number.
:returns: The binary logarithm of this number.
"""
numpy.arctan2(float(self), float(other))
[docs]
def hypot(self, other):
"""This method provides the broadcast interface for
numpy.hypot.
:param other: Another rational number.
:returns: The binary logarithm of this number.
"""
if not isinstance(other, RationalNumber):
tmp, other = self.coerce(other)
return numpy.hypot(tmp, other)
return numpy.sqrt(self * self + other * other)
[docs]
def fmod(self, other):
"""Return ``numpy.fmod(float(self), other)``.
Note:
This operation is intentionally one-way and may fail when
called via reverse dispatch from other numeric types.
"""
return numpy.fmod(float(self), other)
[docs]
def conjugate(self):
"""This method provides the broadcast interface for
numpy.conjugate.
:returns: This number.
"""
return float(self)
[docs]
def coerce(self, other):
"""Implementation of coercion rules.
See also the "Coercion" documentation page.
"""
# A x A
if isinstance(other, RationalNumber):
return (self, other)
elif isinstance(other, int) or other.is_integer():
return (self, RationalNumber.value_of(int(other)))
# fall back to float
else:
ret = float(self)
return Coerce(ret, other)
## \brief Global constant for infinity that is used in combination with the
# degrees of freedom evaluation.
INFINITY = "inf"
## @}