Add files using upload-large-folder tool

lib/python3.1/site-packages/sympy/algebras/__init__.py

@@ -0,0 +1,3 @@
+from .quaternion import Quaternion
+
+__all__ = ["Quaternion",]

lib/python3.1/site-packages/sympy/algebras/quaternion.py

@@ -0,0 +1,1667 @@
+from sympy.core.numbers import Rational
+from sympy.core.singleton import S
+from sympy.core.relational import is_eq
+from sympy.functions.elementary.complexes import (conjugate, im, re, sign)
+from sympy.functions.elementary.exponential import (exp, log as ln)
+from sympy.functions.elementary.miscellaneous import sqrt
+from sympy.functions.elementary.trigonometric import (acos, asin, atan2)
+from sympy.functions.elementary.trigonometric import (cos, sin)
+from sympy.simplify.trigsimp import trigsimp
+from sympy.integrals.integrals import integrate
+from sympy.matrices.dense import MutableDenseMatrix as Matrix
+from sympy.core.sympify import sympify, _sympify
+from sympy.core.expr import Expr
+from sympy.core.logic import fuzzy_not, fuzzy_or
+from sympy.utilities.misc import as_int
+
+from mpmath.libmp.libmpf import prec_to_dps
+
+
+def _check_norm(elements, norm):
+    """validate if input norm is consistent"""
+    if norm is not None and norm.is_number:
+        if norm.is_positive is False:
+            raise ValueError("Input norm must be positive.")
+
+        numerical = all(i.is_number and i.is_real is True for i in elements)
+        if numerical and is_eq(norm**2, sum(i**2 for i in elements)) is False:
+            raise ValueError("Incompatible value for norm.")
+
+
+def _is_extrinsic(seq):
+    """validate seq and return True if seq is lowercase and False if uppercase"""
+    if type(seq) != str:
+        raise ValueError('Expected seq to be a string.')
+    if len(seq) != 3:
+        raise ValueError("Expected 3 axes, got `{}`.".format(seq))
+
+    intrinsic = seq.isupper()
+    extrinsic = seq.islower()
+    if not (intrinsic or extrinsic):
+        raise ValueError("seq must either be fully uppercase (for extrinsic "
+                         "rotations), or fully lowercase, for intrinsic "
+                         "rotations).")
+
+    i, j, k = seq.lower()
+    if (i == j) or (j == k):
+        raise ValueError("Consecutive axes must be different")
+
+    bad = set(seq) - set('xyzXYZ')
+    if bad:
+        raise ValueError("Expected axes from `seq` to be from "
+                         "['x', 'y', 'z'] or ['X', 'Y', 'Z'], "
+                         "got {}".format(''.join(bad)))
+
+    return extrinsic
+
+
+class Quaternion(Expr):
+    """Provides basic quaternion operations.
+    Quaternion objects can be instantiated as ``Quaternion(a, b, c, d)``
+    as in $q = a + bi + cj + dk$.
+
+    Parameters
+    ==========
+
+    norm : None or number
+        Pre-defined quaternion norm. If a value is given, Quaternion.norm
+        returns this pre-defined value instead of calculating the norm
+
+    Examples
+    ========
+
+    >>> from sympy import Quaternion
+    >>> q = Quaternion(1, 2, 3, 4)
+    >>> q
+    1 + 2*i + 3*j + 4*k
+
+    Quaternions over complex fields can be defined as:
+
+    >>> from sympy import Quaternion
+    >>> from sympy import symbols, I
+    >>> x = symbols('x')
+    >>> q1 = Quaternion(x, x**3, x, x**2, real_field = False)
+    >>> q2 = Quaternion(3 + 4*I, 2 + 5*I, 0, 7 + 8*I, real_field = False)
+    >>> q1
+    x + x**3*i + x*j + x**2*k
+    >>> q2
+    (3 + 4*I) + (2 + 5*I)*i + 0*j + (7 + 8*I)*k
+
+    Defining symbolic unit quaternions:
+
+    >>> from sympy import Quaternion
+    >>> from sympy.abc import w, x, y, z
+    >>> q = Quaternion(w, x, y, z, norm=1)
+    >>> q
+    w + x*i + y*j + z*k
+    >>> q.norm()
+    1
+
+    References
+    ==========
+
+    .. [1] https://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/
+    .. [2] https://en.wikipedia.org/wiki/Quaternion
+
+    """
+    _op_priority = 11.0
+
+    is_commutative = False
+
+    def __new__(cls, a=0, b=0, c=0, d=0, real_field=True, norm=None):
+        a, b, c, d = map(sympify, (a, b, c, d))
+
+        if any(i.is_commutative is False for i in [a, b, c, d]):
+            raise ValueError("arguments have to be commutative")
+        obj = super().__new__(cls, a, b, c, d)
+        obj._real_field = real_field
+        obj.set_norm(norm)
+        return obj
+
+    def set_norm(self, norm):
+        """Sets norm of an already instantiated quaternion.
+
+        Parameters
+        ==========
+
+        norm : None or number
+            Pre-defined quaternion norm. If a value is given, Quaternion.norm
+            returns this pre-defined value instead of calculating the norm
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy.abc import a, b, c, d
+        >>> q = Quaternion(a, b, c, d)
+        >>> q.norm()
+        sqrt(a**2 + b**2 + c**2 + d**2)
+
+        Setting the norm:
+
+        >>> q.set_norm(1)
+        >>> q.norm()
+        1
+
+        Removing set norm:
+
+        >>> q.set_norm(None)
+        >>> q.norm()
+        sqrt(a**2 + b**2 + c**2 + d**2)
+
+        """
+        norm = sympify(norm)
+        _check_norm(self.args, norm)
+        self._norm = norm
+
+    @property
+    def a(self):
+        return self.args[0]
+
+    @property
+    def b(self):
+        return self.args[1]
+
+    @property
+    def c(self):
+        return self.args[2]
+
+    @property
+    def d(self):
+        return self.args[3]
+
+    @property
+    def real_field(self):
+        return self._real_field
+
+    @property
+    def product_matrix_left(self):
+        r"""Returns 4 x 4 Matrix equivalent to a Hamilton product from the
+        left. This can be useful when treating quaternion elements as column
+        vectors. Given a quaternion $q = a + bi + cj + dk$ where a, b, c and d
+        are real numbers, the product matrix from the left is:
+
+        .. math::
+
+            M  =  \begin{bmatrix} a  &-b  &-c  &-d \\
+                                  b  & a  &-d  & c \\
+                                  c  & d  & a  &-b \\
+                                  d  &-c  & b  & a \end{bmatrix}
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy.abc import a, b, c, d
+        >>> q1 = Quaternion(1, 0, 0, 1)
+        >>> q2 = Quaternion(a, b, c, d)
+        >>> q1.product_matrix_left
+        Matrix([
+        [1, 0,  0, -1],
+        [0, 1, -1,  0],
+        [0, 1,  1,  0],
+        [1, 0,  0,  1]])
+
+        >>> q1.product_matrix_left * q2.to_Matrix()
+        Matrix([
+        [a - d],
+        [b - c],
+        [b + c],
+        [a + d]])
+
+        This is equivalent to:
+
+        >>> (q1 * q2).to_Matrix()
+        Matrix([
+        [a - d],
+        [b - c],
+        [b + c],
+        [a + d]])
+        """
+        return Matrix([
+                [self.a, -self.b, -self.c, -self.d],
+                [self.b, self.a, -self.d, self.c],
+                [self.c, self.d, self.a, -self.b],
+                [self.d, -self.c, self.b, self.a]])
+
+    @property
+    def product_matrix_right(self):
+        r"""Returns 4 x 4 Matrix equivalent to a Hamilton product from the
+        right. This can be useful when treating quaternion elements as column
+        vectors. Given a quaternion $q = a + bi + cj + dk$ where a, b, c and d
+        are real numbers, the product matrix from the left is:
+
+        .. math::
+
+            M  =  \begin{bmatrix} a  &-b  &-c  &-d \\
+                                  b  & a  & d  &-c \\
+                                  c  &-d  & a  & b \\
+                                  d  & c  &-b  & a \end{bmatrix}
+
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy.abc import a, b, c, d
+        >>> q1 = Quaternion(a, b, c, d)
+        >>> q2 = Quaternion(1, 0, 0, 1)
+        >>> q2.product_matrix_right
+        Matrix([
+        [1, 0, 0, -1],
+        [0, 1, 1, 0],
+        [0, -1, 1, 0],
+        [1, 0, 0, 1]])
+
+        Note the switched arguments: the matrix represents the quaternion on
+        the right, but is still considered as a matrix multiplication from the
+        left.
+
+        >>> q2.product_matrix_right * q1.to_Matrix()
+        Matrix([
+        [ a - d],
+        [ b + c],
+        [-b + c],
+        [ a + d]])
+
+        This is equivalent to:
+
+        >>> (q1 * q2).to_Matrix()
+        Matrix([
+        [ a - d],
+        [ b + c],
+        [-b + c],
+        [ a + d]])
+        """
+        return Matrix([
+                [self.a, -self.b, -self.c, -self.d],
+                [self.b, self.a, self.d, -self.c],
+                [self.c, -self.d, self.a, self.b],
+                [self.d, self.c, -self.b, self.a]])
+
+    def to_Matrix(self, vector_only=False):
+        """Returns elements of quaternion as a column vector.
+        By default, a ``Matrix`` of length 4 is returned, with the real part as the
+        first element.
+        If ``vector_only`` is ``True``, returns only imaginary part as a Matrix of
+        length 3.
+
+        Parameters
+        ==========
+
+        vector_only : bool
+            If True, only imaginary part is returned.
+            Default value: False
+
+        Returns
+        =======
+
+        Matrix
+            A column vector constructed by the elements of the quaternion.
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy.abc import a, b, c, d
+        >>> q = Quaternion(a, b, c, d)
+        >>> q
+        a + b*i + c*j + d*k
+
+        >>> q.to_Matrix()
+        Matrix([
+        [a],
+        [b],
+        [c],
+        [d]])
+
+
+        >>> q.to_Matrix(vector_only=True)
+        Matrix([
+        [b],
+        [c],
+        [d]])
+
+        """
+        if vector_only:
+            return Matrix(self.args[1:])
+        else:
+            return Matrix(self.args)
+
+    @classmethod
+    def from_Matrix(cls, elements):
+        """Returns quaternion from elements of a column vector`.
+        If vector_only is True, returns only imaginary part as a Matrix of
+        length 3.
+
+        Parameters
+        ==========
+
+        elements : Matrix, list or tuple of length 3 or 4. If length is 3,
+            assume real part is zero.
+            Default value: False
+
+        Returns
+        =======
+
+        Quaternion
+            A quaternion created from the input elements.
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy.abc import a, b, c, d
+        >>> q = Quaternion.from_Matrix([a, b, c, d])
+        >>> q
+        a + b*i + c*j + d*k
+
+        >>> q = Quaternion.from_Matrix([b, c, d])
+        >>> q
+        0 + b*i + c*j + d*k
+
+        """
+        length = len(elements)
+        if length != 3 and length != 4:
+            raise ValueError("Input elements must have length 3 or 4, got {} "
+                             "elements".format(length))
+
+        if length == 3:
+            return Quaternion(0, *elements)
+        else:
+            return Quaternion(*elements)
+
+    @classmethod
+    def from_euler(cls, angles, seq):
+        """Returns quaternion equivalent to rotation represented by the Euler
+        angles, in the sequence defined by ``seq``.
+
+        Parameters
+        ==========
+
+        angles : list, tuple or Matrix of 3 numbers
+            The Euler angles (in radians).
+        seq : string of length 3
+            Represents the sequence of rotations.
+            For extrinsic rotations, seq must be all lowercase and its elements
+            must be from the set ``{'x', 'y', 'z'}``
+            For intrinsic rotations, seq must be all uppercase and its elements
+            must be from the set ``{'X', 'Y', 'Z'}``
+
+        Returns
+        =======
+
+        Quaternion
+            The normalized rotation quaternion calculated from the Euler angles
+            in the given sequence.
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy import pi
+        >>> q = Quaternion.from_euler([pi/2, 0, 0], 'xyz')
+        >>> q
+        sqrt(2)/2 + sqrt(2)/2*i + 0*j + 0*k
+
+        >>> q = Quaternion.from_euler([0, pi/2, pi] , 'zyz')
+        >>> q
+        0 + (-sqrt(2)/2)*i + 0*j + sqrt(2)/2*k
+
+        >>> q = Quaternion.from_euler([0, pi/2, pi] , 'ZYZ')
+        >>> q
+        0 + sqrt(2)/2*i + 0*j + sqrt(2)/2*k
+
+        """
+
+        if len(angles) != 3:
+            raise ValueError("3 angles must be given.")
+
+        extrinsic = _is_extrinsic(seq)
+        i, j, k = seq.lower()
+
+        # get elementary basis vectors
+        ei = [1 if n == i else 0 for n in 'xyz']
+        ej = [1 if n == j else 0 for n in 'xyz']
+        ek = [1 if n == k else 0 for n in 'xyz']
+
+        # calculate distinct quaternions
+        qi = cls.from_axis_angle(ei, angles[0])
+        qj = cls.from_axis_angle(ej, angles[1])
+        qk = cls.from_axis_angle(ek, angles[2])
+
+        if extrinsic:
+            return trigsimp(qk * qj * qi)
+        else:
+            return trigsimp(qi * qj * qk)
+
+    def to_euler(self, seq, angle_addition=True, avoid_square_root=False):
+        r"""Returns Euler angles representing same rotation as the quaternion,
+        in the sequence given by ``seq``. This implements the method described
+        in [1]_.
+
+        For degenerate cases (gymbal lock cases), the third angle is
+        set to zero.
+
+        Parameters
+        ==========
+
+        seq : string of length 3
+            Represents the sequence of rotations.
+            For extrinsic rotations, seq must be all lowercase and its elements
+            must be from the set ``{'x', 'y', 'z'}``
+            For intrinsic rotations, seq must be all uppercase and its elements
+            must be from the set ``{'X', 'Y', 'Z'}``
+
+        angle_addition : bool
+            When True, first and third angles are given as an addition and
+            subtraction of two simpler ``atan2`` expressions. When False, the
+            first and third angles are each given by a single more complicated
+            ``atan2`` expression. This equivalent expression is given by:
+
+            .. math::
+
+                \operatorname{atan_2} (b,a) \pm \operatorname{atan_2} (d,c) =
+                \operatorname{atan_2} (bc\pm ad, ac\mp bd)
+
+            Default value: True
+
+        avoid_square_root : bool
+            When True, the second angle is calculated with an expression based
+            on ``acos``, which is slightly more complicated but avoids a square
+            root. When False, second angle is calculated with ``atan2``, which
+            is simpler and can be better for numerical reasons (some
+            numerical implementations of ``acos`` have problems near zero).
+            Default value: False
+
+
+        Returns
+        =======
+
+        Tuple
+            The Euler angles calculated from the quaternion
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy.abc import a, b, c, d
+        >>> euler = Quaternion(a, b, c, d).to_euler('zyz')
+        >>> euler
+        (-atan2(-b, c) + atan2(d, a),
+         2*atan2(sqrt(b**2 + c**2), sqrt(a**2 + d**2)),
+         atan2(-b, c) + atan2(d, a))
+
+
+        References
+        ==========
+
+        .. [1] https://doi.org/10.1371/journal.pone.0276302
+
+        """
+        if self.is_zero_quaternion():
+            raise ValueError('Cannot convert a quaternion with norm 0.')
+
+        angles = [0, 0, 0]
+
+        extrinsic = _is_extrinsic(seq)
+        i, j, k = seq.lower()
+
+        # get index corresponding to elementary basis vectors
+        i = 'xyz'.index(i) + 1
+        j = 'xyz'.index(j) + 1
+        k = 'xyz'.index(k) + 1
+
+        if not extrinsic:
+            i, k = k, i
+
+        # check if sequence is symmetric
+        symmetric = i == k
+        if symmetric:
+            k = 6 - i - j
+
+        # parity of the permutation
+        sign = (i - j) * (j - k) * (k - i) // 2
+
+        # permutate elements
+        elements = [self.a, self.b, self.c, self.d]
+        a = elements[0]
+        b = elements[i]
+        c = elements[j]
+        d = elements[k] * sign
+
+        if not symmetric:
+            a, b, c, d = a - c, b + d, c + a, d - b
+
+        if avoid_square_root:
+            if symmetric:
+                n2 = self.norm()**2
+                angles[1] = acos((a * a + b * b - c * c - d * d) / n2)
+            else:
+                n2 = 2 * self.norm()**2
+                angles[1] = asin((c * c + d * d - a * a - b * b) / n2)
+        else:
+            angles[1] = 2 * atan2(sqrt(c * c + d * d), sqrt(a * a + b * b))
+            if not symmetric:
+                angles[1] -= S.Pi / 2
+
+        # Check for singularities in numerical cases
+        case = 0
+        if is_eq(c, S.Zero) and is_eq(d, S.Zero):
+            case = 1
+        if is_eq(a, S.Zero) and is_eq(b, S.Zero):
+            case = 2
+
+        if case == 0:
+            if angle_addition:
+                angles[0] = atan2(b, a) + atan2(d, c)
+                angles[2] = atan2(b, a) - atan2(d, c)
+            else:
+                angles[0] = atan2(b*c + a*d, a*c - b*d)
+                angles[2] = atan2(b*c - a*d, a*c + b*d)
+
+        else:  # any degenerate case
+            angles[2 * (not extrinsic)] = S.Zero
+            if case == 1:
+                angles[2 * extrinsic] = 2 * atan2(b, a)
+            else:
+                angles[2 * extrinsic] = 2 * atan2(d, c)
+                angles[2 * extrinsic] *= (-1 if extrinsic else 1)
+
+        # for Tait-Bryan angles
+        if not symmetric:
+            angles[0] *= sign
+
+        if extrinsic:
+            return tuple(angles[::-1])
+        else:
+            return tuple(angles)
+
+    @classmethod
+    def from_axis_angle(cls, vector, angle):
+        """Returns a rotation quaternion given the axis and the angle of rotation.
+
+        Parameters
+        ==========
+
+        vector : tuple of three numbers
+            The vector representation of the given axis.
+        angle : number
+            The angle by which axis is rotated (in radians).
+
+        Returns
+        =======
+
+        Quaternion
+            The normalized rotation quaternion calculated from the given axis and the angle of rotation.
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy import pi, sqrt
+        >>> q = Quaternion.from_axis_angle((sqrt(3)/3, sqrt(3)/3, sqrt(3)/3), 2*pi/3)
+        >>> q
+        1/2 + 1/2*i + 1/2*j + 1/2*k
+
+        """
+        (x, y, z) = vector
+        norm = sqrt(x**2 + y**2 + z**2)
+        (x, y, z) = (x / norm, y / norm, z / norm)
+        s = sin(angle * S.Half)
+        a = cos(angle * S.Half)
+        b = x * s
+        c = y * s
+        d = z * s
+
+        # note that this quaternion is already normalized by construction:
+        # c^2 + (s*x)^2 + (s*y)^2 + (s*z)^2 = c^2 + s^2*(x^2 + y^2 + z^2) = c^2 + s^2 * 1 = c^2 + s^2 = 1
+        # so, what we return is a normalized quaternion
+
+        return cls(a, b, c, d)
+
+    @classmethod
+    def from_rotation_matrix(cls, M):
+        """Returns the equivalent quaternion of a matrix. The quaternion will be normalized
+        only if the matrix is special orthogonal (orthogonal and det(M) = 1).
+
+        Parameters
+        ==========
+
+        M : Matrix
+            Input matrix to be converted to equivalent quaternion. M must be special
+            orthogonal (orthogonal and det(M) = 1) for the quaternion to be normalized.
+
+        Returns
+        =======
+
+        Quaternion
+            The quaternion equivalent to given matrix.
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy import Matrix, symbols, cos, sin, trigsimp
+        >>> x = symbols('x')
+        >>> M = Matrix([[cos(x), -sin(x), 0], [sin(x), cos(x), 0], [0, 0, 1]])
+        >>> q = trigsimp(Quaternion.from_rotation_matrix(M))
+        >>> q
+        sqrt(2)*sqrt(cos(x) + 1)/2 + 0*i + 0*j + sqrt(2 - 2*cos(x))*sign(sin(x))/2*k
+
+        """
+
+        absQ = M.det()**Rational(1, 3)
+
+        a = sqrt(absQ + M[0, 0] + M[1, 1] + M[2, 2]) / 2
+        b = sqrt(absQ + M[0, 0] - M[1, 1] - M[2, 2]) / 2
+        c = sqrt(absQ - M[0, 0] + M[1, 1] - M[2, 2]) / 2
+        d = sqrt(absQ - M[0, 0] - M[1, 1] + M[2, 2]) / 2
+
+        b = b * sign(M[2, 1] - M[1, 2])
+        c = c * sign(M[0, 2] - M[2, 0])
+        d = d * sign(M[1, 0] - M[0, 1])
+
+        return Quaternion(a, b, c, d)
+
+    def __add__(self, other):
+        return self.add(other)
+
+    def __radd__(self, other):
+        return self.add(other)
+
+    def __sub__(self, other):
+        return self.add(other*-1)
+
+    def __mul__(self, other):
+        return self._generic_mul(self, _sympify(other))
+
+    def __rmul__(self, other):
+        return self._generic_mul(_sympify(other), self)
+
+    def __pow__(self, p):
+        return self.pow(p)
+
+    def __neg__(self):
+        return Quaternion(-self.a, -self.b, -self.c, -self.d)
+
+    def __truediv__(self, other):
+        return self * sympify(other)**-1
+
+    def __rtruediv__(self, other):
+        return sympify(other) * self**-1
+
+    def _eval_Integral(self, *args):
+        return self.integrate(*args)
+
+    def diff(self, *symbols, **kwargs):
+        kwargs.setdefault('evaluate', True)
+        return self.func(*[a.diff(*symbols, **kwargs) for a  in self.args])
+
+    def add(self, other):
+        """Adds quaternions.
+
+        Parameters
+        ==========
+
+        other : Quaternion
+            The quaternion to add to current (self) quaternion.
+
+        Returns
+        =======
+
+        Quaternion
+            The resultant quaternion after adding self to other
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy import symbols
+        >>> q1 = Quaternion(1, 2, 3, 4)
+        >>> q2 = Quaternion(5, 6, 7, 8)
+        >>> q1.add(q2)
+        6 + 8*i + 10*j + 12*k
+        >>> q1 + 5
+        6 + 2*i + 3*j + 4*k
+        >>> x = symbols('x', real = True)
+        >>> q1.add(x)
+        (x + 1) + 2*i + 3*j + 4*k
+
+        Quaternions over complex fields :
+
+        >>> from sympy import Quaternion
+        >>> from sympy import I
+        >>> q3 = Quaternion(3 + 4*I, 2 + 5*I, 0, 7 + 8*I, real_field = False)
+        >>> q3.add(2 + 3*I)
+        (5 + 7*I) + (2 + 5*I)*i + 0*j + (7 + 8*I)*k
+
+        """
+        q1 = self
+        q2 = sympify(other)
+
+        # If q2 is a number or a SymPy expression instead of a quaternion
+        if not isinstance(q2, Quaternion):
+            if q1.real_field and q2.is_complex:
+                return Quaternion(re(q2) + q1.a, im(q2) + q1.b, q1.c, q1.d)
+            elif q2.is_commutative:
+                return Quaternion(q1.a + q2, q1.b, q1.c, q1.d)
+            else:
+                raise ValueError("Only commutative expressions can be added with a Quaternion.")
+
+        return Quaternion(q1.a + q2.a, q1.b + q2.b, q1.c + q2.c, q1.d
+                          + q2.d)
+
+    def mul(self, other):
+        """Multiplies quaternions.
+
+        Parameters
+        ==========
+
+        other : Quaternion or symbol
+            The quaternion to multiply to current (self) quaternion.
+
+        Returns
+        =======
+
+        Quaternion
+            The resultant quaternion after multiplying self with other
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy import symbols
+        >>> q1 = Quaternion(1, 2, 3, 4)
+        >>> q2 = Quaternion(5, 6, 7, 8)
+        >>> q1.mul(q2)
+        (-60) + 12*i + 30*j + 24*k
+        >>> q1.mul(2)
+        2 + 4*i + 6*j + 8*k
+        >>> x = symbols('x', real = True)
+        >>> q1.mul(x)
+        x + 2*x*i + 3*x*j + 4*x*k
+
+        Quaternions over complex fields :
+
+        >>> from sympy import Quaternion
+        >>> from sympy import I
+        >>> q3 = Quaternion(3 + 4*I, 2 + 5*I, 0, 7 + 8*I, real_field = False)
+        >>> q3.mul(2 + 3*I)
+        (2 + 3*I)*(3 + 4*I) + (2 + 3*I)*(2 + 5*I)*i + 0*j + (2 + 3*I)*(7 + 8*I)*k
+
+        """
+        return self._generic_mul(self, _sympify(other))
+
+    @staticmethod
+    def _generic_mul(q1, q2):
+        """Generic multiplication.
+
+        Parameters
+        ==========
+
+        q1 : Quaternion or symbol
+        q2 : Quaternion or symbol
+
+        It is important to note that if neither q1 nor q2 is a Quaternion,
+        this function simply returns q1 * q2.
+
+        Returns
+        =======
+
+        Quaternion
+            The resultant quaternion after multiplying q1 and q2
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy import Symbol, S
+        >>> q1 = Quaternion(1, 2, 3, 4)
+        >>> q2 = Quaternion(5, 6, 7, 8)
+        >>> Quaternion._generic_mul(q1, q2)
+        (-60) + 12*i + 30*j + 24*k
+        >>> Quaternion._generic_mul(q1, S(2))
+        2 + 4*i + 6*j + 8*k
+        >>> x = Symbol('x', real = True)
+        >>> Quaternion._generic_mul(q1, x)
+        x + 2*x*i + 3*x*j + 4*x*k
+
+        Quaternions over complex fields :
+
+        >>> from sympy import I
+        >>> q3 = Quaternion(3 + 4*I, 2 + 5*I, 0, 7 + 8*I, real_field = False)
+        >>> Quaternion._generic_mul(q3, 2 + 3*I)
+        (2 + 3*I)*(3 + 4*I) + (2 + 3*I)*(2 + 5*I)*i + 0*j + (2 + 3*I)*(7 + 8*I)*k
+
+        """
+        # None is a Quaternion:
+        if not isinstance(q1, Quaternion) and not isinstance(q2, Quaternion):
+            return q1 * q2
+
+        # If q1 is a number or a SymPy expression instead of a quaternion
+        if not isinstance(q1, Quaternion):
+            if q2.real_field and q1.is_complex:
+                return Quaternion(re(q1), im(q1), 0, 0) * q2
+            elif q1.is_commutative:
+                return Quaternion(q1 * q2.a, q1 * q2.b, q1 * q2.c, q1 * q2.d)
+            else:
+                raise ValueError("Only commutative expressions can be multiplied with a Quaternion.")
+
+        # If q2 is a number or a SymPy expression instead of a quaternion
+        if not isinstance(q2, Quaternion):
+            if q1.real_field and q2.is_complex:
+                return q1 * Quaternion(re(q2), im(q2), 0, 0)
+            elif q2.is_commutative:
+                return Quaternion(q2 * q1.a, q2 * q1.b, q2 * q1.c, q2 * q1.d)
+            else:
+                raise ValueError("Only commutative expressions can be multiplied with a Quaternion.")
+
+        # If any of the quaternions has a fixed norm, pre-compute norm
+        if q1._norm is None and q2._norm is None:
+            norm = None
+        else:
+            norm = q1.norm() * q2.norm()
+
+        return Quaternion(-q1.b*q2.b - q1.c*q2.c - q1.d*q2.d + q1.a*q2.a,
+                          q1.b*q2.a + q1.c*q2.d - q1.d*q2.c + q1.a*q2.b,
+                          -q1.b*q2.d + q1.c*q2.a + q1.d*q2.b + q1.a*q2.c,
+                          q1.b*q2.c - q1.c*q2.b + q1.d*q2.a + q1.a * q2.d,
+                          norm=norm)
+
+    def _eval_conjugate(self):
+        """Returns the conjugate of the quaternion."""
+        q = self
+        return Quaternion(q.a, -q.b, -q.c, -q.d, norm=q._norm)
+
+    def norm(self):
+        """Returns the norm of the quaternion."""
+        if self._norm is None:  # check if norm is pre-defined
+            q = self
+            # trigsimp is used to simplify sin(x)^2 + cos(x)^2 (these terms
+            # arise when from_axis_angle is used).
+            return sqrt(trigsimp(q.a**2 + q.b**2 + q.c**2 + q.d**2))
+
+        return self._norm
+
+    def normalize(self):
+        """Returns the normalized form of the quaternion."""
+        q = self
+        return q * (1/q.norm())
+
+    def inverse(self):
+        """Returns the inverse of the quaternion."""
+        q = self
+        if not q.norm():
+            raise ValueError("Cannot compute inverse for a quaternion with zero norm")
+        return conjugate(q) * (1/q.norm()**2)
+
+    def pow(self, p):
+        """Finds the pth power of the quaternion.
+
+        Parameters
+        ==========
+
+        p : int
+            Power to be applied on quaternion.
+
+        Returns
+        =======
+
+        Quaternion
+            Returns the p-th power of the current quaternion.
+            Returns the inverse if p = -1.
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> q = Quaternion(1, 2, 3, 4)
+        >>> q.pow(4)
+        668 + (-224)*i + (-336)*j + (-448)*k
+
+        """
+        try:
+            q, p = self, as_int(p)
+        except ValueError:
+            return NotImplemented
+
+        if p < 0:
+            q, p = q.inverse(), -p
+
+        if p == 1:
+            return q
+
+        res = Quaternion(1, 0, 0, 0)
+        while p > 0:
+            if p & 1:
+                res *= q
+            q *= q
+            p >>= 1
+
+        return res
+
+    def exp(self):
+        """Returns the exponential of $q$, given by $e^q$.
+
+        Returns
+        =======
+
+        Quaternion
+            The exponential of the quaternion.
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> q = Quaternion(1, 2, 3, 4)
+        >>> q.exp()
+        E*cos(sqrt(29))
+        + 2*sqrt(29)*E*sin(sqrt(29))/29*i
+        + 3*sqrt(29)*E*sin(sqrt(29))/29*j
+        + 4*sqrt(29)*E*sin(sqrt(29))/29*k
+
+        """
+        # exp(q) = e^a(cos||v|| + v/||v||*sin||v||)
+        q = self
+        vector_norm = sqrt(q.b**2 + q.c**2 + q.d**2)
+        a = exp(q.a) * cos(vector_norm)
+        b = exp(q.a) * sin(vector_norm) * q.b / vector_norm
+        c = exp(q.a) * sin(vector_norm) * q.c / vector_norm
+        d = exp(q.a) * sin(vector_norm) * q.d / vector_norm
+
+        return Quaternion(a, b, c, d)
+
+    def log(self):
+        r"""Returns the logarithm of the quaternion, given by $\log q$.
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> q = Quaternion(1, 2, 3, 4)
+        >>> q.log()
+        log(sqrt(30))
+        + 2*sqrt(29)*acos(sqrt(30)/30)/29*i
+        + 3*sqrt(29)*acos(sqrt(30)/30)/29*j
+        + 4*sqrt(29)*acos(sqrt(30)/30)/29*k
+
+        """
+        # log(q) = log||q|| + v/||v||*arccos(a/||q||)
+        q = self
+        vector_norm = sqrt(q.b**2 + q.c**2 + q.d**2)
+        q_norm = q.norm()
+        a = ln(q_norm)
+        b = q.b * acos(q.a / q_norm) / vector_norm
+        c = q.c * acos(q.a / q_norm) / vector_norm
+        d = q.d * acos(q.a / q_norm) / vector_norm
+
+        return Quaternion(a, b, c, d)
+
+    def _eval_subs(self, *args):
+        elements = [i.subs(*args) for i in self.args]
+        norm = self._norm
+        if norm is not None:
+            norm = norm.subs(*args)
+        _check_norm(elements, norm)
+        return Quaternion(*elements, norm=norm)
+
+    def _eval_evalf(self, prec):
+        """Returns the floating point approximations (decimal numbers) of the quaternion.
+
+        Returns
+        =======
+
+        Quaternion
+            Floating point approximations of quaternion(self)
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy import sqrt
+        >>> q = Quaternion(1/sqrt(1), 1/sqrt(2), 1/sqrt(3), 1/sqrt(4))
+        >>> q.evalf()
+        1.00000000000000
+        + 0.707106781186547*i
+        + 0.577350269189626*j
+        + 0.500000000000000*k
+
+        """
+        nprec = prec_to_dps(prec)
+        return Quaternion(*[arg.evalf(n=nprec) for arg in self.args])
+
+    def pow_cos_sin(self, p):
+        """Computes the pth power in the cos-sin form.
+
+        Parameters
+        ==========
+
+        p : int
+            Power to be applied on quaternion.
+
+        Returns
+        =======
+
+        Quaternion
+            The p-th power in the cos-sin form.
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> q = Quaternion(1, 2, 3, 4)
+        >>> q.pow_cos_sin(4)
+        900*cos(4*acos(sqrt(30)/30))
+        + 1800*sqrt(29)*sin(4*acos(sqrt(30)/30))/29*i
+        + 2700*sqrt(29)*sin(4*acos(sqrt(30)/30))/29*j
+        + 3600*sqrt(29)*sin(4*acos(sqrt(30)/30))/29*k
+
+        """
+        # q = ||q||*(cos(a) + u*sin(a))
+        # q^p = ||q||^p * (cos(p*a) + u*sin(p*a))
+
+        q = self
+        (v, angle) = q.to_axis_angle()
+        q2 = Quaternion.from_axis_angle(v, p * angle)
+        return q2 * (q.norm()**p)
+
+    def integrate(self, *args):
+        """Computes integration of quaternion.
+
+        Returns
+        =======
+
+        Quaternion
+            Integration of the quaternion(self) with the given variable.
+
+        Examples
+        ========
+
+        Indefinite Integral of quaternion :
+
+        >>> from sympy import Quaternion
+        >>> from sympy.abc import x
+        >>> q = Quaternion(1, 2, 3, 4)
+        >>> q.integrate(x)
+        x + 2*x*i + 3*x*j + 4*x*k
+
+        Definite integral of quaternion :
+
+        >>> from sympy import Quaternion
+        >>> from sympy.abc import x
+        >>> q = Quaternion(1, 2, 3, 4)
+        >>> q.integrate((x, 1, 5))
+        4 + 8*i + 12*j + 16*k
+
+        """
+        # TODO: is this expression correct?
+        return Quaternion(integrate(self.a, *args), integrate(self.b, *args),
+                          integrate(self.c, *args), integrate(self.d, *args))
+
+    @staticmethod
+    def rotate_point(pin, r):
+        """Returns the coordinates of the point pin (a 3 tuple) after rotation.
+
+        Parameters
+        ==========
+
+        pin : tuple
+            A 3-element tuple of coordinates of a point which needs to be
+            rotated.
+        r : Quaternion or tuple
+            Axis and angle of rotation.
+
+            It's important to note that when r is a tuple, it must be of the form
+            (axis, angle)
+
+        Returns
+        =======
+
+        tuple
+            The coordinates of the point after rotation.
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy import symbols, trigsimp, cos, sin
+        >>> x = symbols('x')
+        >>> q = Quaternion(cos(x/2), 0, 0, sin(x/2))
+        >>> trigsimp(Quaternion.rotate_point((1, 1, 1), q))
+        (sqrt(2)*cos(x + pi/4), sqrt(2)*sin(x + pi/4), 1)
+        >>> (axis, angle) = q.to_axis_angle()
+        >>> trigsimp(Quaternion.rotate_point((1, 1, 1), (axis, angle)))
+        (sqrt(2)*cos(x + pi/4), sqrt(2)*sin(x + pi/4), 1)
+
+        """
+        if isinstance(r, tuple):
+            # if r is of the form (vector, angle)
+            q = Quaternion.from_axis_angle(r[0], r[1])
+        else:
+            # if r is a quaternion
+            q = r.normalize()
+        pout = q * Quaternion(0, pin[0], pin[1], pin[2]) * conjugate(q)
+        return (pout.b, pout.c, pout.d)
+
+    def to_axis_angle(self):
+        """Returns the axis and angle of rotation of a quaternion.
+
+        Returns
+        =======
+
+        tuple
+            Tuple of (axis, angle)
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> q = Quaternion(1, 1, 1, 1)
+        >>> (axis, angle) = q.to_axis_angle()
+        >>> axis
+        (sqrt(3)/3, sqrt(3)/3, sqrt(3)/3)
+        >>> angle
+        2*pi/3
+
+        """
+        q = self
+        if q.a.is_negative:
+            q = q * -1
+
+        q = q.normalize()
+        angle = trigsimp(2 * acos(q.a))
+
+        # Since quaternion is normalised, q.a is less than 1.
+        s = sqrt(1 - q.a*q.a)
+
+        x = trigsimp(q.b / s)
+        y = trigsimp(q.c / s)
+        z = trigsimp(q.d / s)
+
+        v = (x, y, z)
+        t = (v, angle)
+
+        return t
+
+    def to_rotation_matrix(self, v=None, homogeneous=True):
+        """Returns the equivalent rotation transformation matrix of the quaternion
+        which represents rotation about the origin if ``v`` is not passed.
+
+        Parameters
+        ==========
+
+        v : tuple or None
+            Default value: None
+        homogeneous : bool
+            When True, gives an expression that may be more efficient for
+            symbolic calculations but less so for direct evaluation. Both
+            formulas are mathematically equivalent.
+            Default value: True
+
+        Returns
+        =======
+
+        tuple
+            Returns the equivalent rotation transformation matrix of the quaternion
+            which represents rotation about the origin if v is not passed.
+
+        Examples
+        ========
+
+        >>> from sympy import Quaternion
+        >>> from sympy import symbols, trigsimp, cos, sin
+        >>> x = symbols('x')
+        >>> q = Quaternion(cos(x/2), 0, 0, sin(x/2))
+        >>> trigsimp(q.to_rotation_matrix())
+        Matrix([
+        [cos(x), -sin(x), 0],
+        [sin(x),  cos(x), 0],
+        [     0,       0, 1]])
+
+        Generates a 4x4 transformation matrix (used for rotation about a point
+        other than the origin) if the point(v) is passed as an argument.
+        """
+
+        q = self
+        s = q.norm()**-2
+
+        # diagonal elements are different according to parameter normal
+        if homogeneous:
+            m00 = s*(q.a**2 + q.b**2 - q.c**2 - q.d**2)
+            m11 = s*(q.a**2 - q.b**2 + q.c**2 - q.d**2)
+            m22 = s*(q.a**2 - q.b**2 - q.c**2 + q.d**2)
+        else:
+            m00 = 1 - 2*s*(q.c**2 + q.d**2)
+            m11 = 1 - 2*s*(q.b**2 + q.d**2)
+            m22 = 1 - 2*s*(q.b**2 + q.c**2)
+
+        m01 = 2*s*(q.b*q.c - q.d*q.a)
+        m02 = 2*s*(q.b*q.d + q.c*q.a)
+
+        m10 = 2*s*(q.b*q.c + q.d*q.a)
+        m12 = 2*s*(q.c*q.d - q.b*q.a)
+
+        m20 = 2*s*(q.b*q.d - q.c*q.a)
+        m21 = 2*s*(q.c*q.d + q.b*q.a)
+
+        if not v:
+            return Matrix([[m00, m01, m02], [m10, m11, m12], [m20, m21, m22]])
+
+        else:
+            (x, y, z) = v
+
+            m03 = x - x*m00 - y*m01 - z*m02
+            m13 = y - x*m10 - y*m11 - z*m12
+            m23 = z - x*m20 - y*m21 - z*m22
+            m30 = m31 = m32 = 0
+            m33 = 1
+
+            return Matrix([[m00, m01, m02, m03], [m10, m11, m12, m13],
+                          [m20, m21, m22, m23], [m30, m31, m32, m33]])
+
+    def scalar_part(self):
+        r"""Returns scalar part($\mathbf{S}(q)$) of the quaternion q.
+
+        Explanation
+        ===========
+
+        Given a quaternion $q = a + bi + cj + dk$, returns $\mathbf{S}(q) = a$.
+
+        Examples
+        ========
+
+        >>> from sympy.algebras.quaternion import Quaternion
+        >>> q = Quaternion(4, 8, 13, 12)
+        >>> q.scalar_part()
+        4
+
+        """
+
+        return self.a
+
+    def vector_part(self):
+        r"""
+        Returns $\mathbf{V}(q)$, the vector part of the quaternion $q$.
+
+        Explanation
+        ===========
+
+        Given a quaternion $q = a + bi + cj + dk$, returns $\mathbf{V}(q) = bi + cj + dk$.
+
+        Examples
+        ========
+
+        >>> from sympy.algebras.quaternion import Quaternion
+        >>> q = Quaternion(1, 1, 1, 1)
+        >>> q.vector_part()
+        0 + 1*i + 1*j + 1*k
+
+        >>> q = Quaternion(4, 8, 13, 12)
+        >>> q.vector_part()
+        0 + 8*i + 13*j + 12*k
+
+        """
+
+        return Quaternion(0, self.b, self.c, self.d)
+
+    def axis(self):
+        r"""
+        Returns $\mathbf{Ax}(q)$, the axis of the quaternion $q$.
+
+        Explanation
+        ===========
+
+        Given a quaternion $q = a + bi + cj + dk$, returns $\mathbf{Ax}(q)$  i.e., the versor of the vector part of that quaternion
+        equal to $\mathbf{U}[\mathbf{V}(q)]$.
+        The axis is always an imaginary unit with square equal to $-1 + 0i + 0j + 0k$.
+
+        Examples
+        ========
+
+        >>> from sympy.algebras.quaternion import Quaternion
+        >>> q = Quaternion(1, 1, 1, 1)
+        >>> q.axis()
+        0 + sqrt(3)/3*i + sqrt(3)/3*j + sqrt(3)/3*k
+
+        See Also
+        ========
+
+        vector_part
+
+        """
+        axis = self.vector_part().normalize()
+
+        return Quaternion(0, axis.b, axis.c, axis.d)
+
+    def is_pure(self):
+        """
+        Returns true if the quaternion is pure, false if the quaternion is not pure
+        or returns none if it is unknown.
+
+        Explanation
+        ===========
+
+        A pure quaternion (also a vector quaternion) is a quaternion with scalar
+        part equal to 0.
+
+        Examples
+        ========
+
+        >>> from sympy.algebras.quaternion import Quaternion
+        >>> q = Quaternion(0, 8, 13, 12)
+        >>> q.is_pure()
+        True
+
+        See Also
+        ========
+        scalar_part
+
+        """
+
+        return self.a.is_zero
+
+    def is_zero_quaternion(self):
+        """
+        Returns true if the quaternion is a zero quaternion or false if it is not a zero quaternion
+        and None if the value is unknown.
+
+        Explanation
+        ===========
+
+        A zero quaternion is a quaternion with both scalar part and
+        vector part equal to 0.
+
+        Examples
+        ========
+
+        >>> from sympy.algebras.quaternion import Quaternion
+        >>> q = Quaternion(1, 0, 0, 0)
+        >>> q.is_zero_quaternion()
+        False
+
+        >>> q = Quaternion(0, 0, 0, 0)
+        >>> q.is_zero_quaternion()
+        True
+
+        See Also
+        ========
+        scalar_part
+        vector_part
+
+        """
+
+        return self.norm().is_zero
+
+    def angle(self):
+        r"""
+        Returns the angle of the quaternion measured in the real-axis plane.
+
+        Explanation
+        ===========
+
+        Given a quaternion $q = a + bi + cj + dk$ where $a$, $b$, $c$ and $d$
+        are real numbers, returns the angle of the quaternion given by
+
+        .. math::
+            \theta := 2 \operatorname{atan_2}\left(\sqrt{b^2 + c^2 + d^2}, {a}\right)
+
+        Examples
+        ========
+
+        >>> from sympy.algebras.quaternion import Quaternion
+        >>> q = Quaternion(1, 4, 4, 4)
+        >>> q.angle()
+        2*atan(4*sqrt(3))
+
+        """
+
+        return 2 * atan2(self.vector_part().norm(), self.scalar_part())
+
+
+    def arc_coplanar(self, other):
+        """
+        Returns True if the transformation arcs represented by the input quaternions happen in the same plane.
+
+        Explanation
+        ===========
+
+        Two quaternions are said to be coplanar (in this arc sense) when their axes are parallel.
+        The plane of a quaternion is the one normal to its axis.
+
+        Parameters
+        ==========
+
+        other : a Quaternion
+
+        Returns
+        =======
+
+        True : if the planes of the two quaternions are the same, apart from its orientation/sign.
+        False : if the planes of the two quaternions are not the same, apart from its orientation/sign.
+        None : if plane of either of the quaternion is unknown.
+
+        Examples
+        ========
+
+        >>> from sympy.algebras.quaternion import Quaternion
+        >>> q1 = Quaternion(1, 4, 4, 4)
+        >>> q2 = Quaternion(3, 8, 8, 8)
+        >>> Quaternion.arc_coplanar(q1, q2)
+        True
+
+        >>> q1 = Quaternion(2, 8, 13, 12)
+        >>> Quaternion.arc_coplanar(q1, q2)
+        False
+
+        See Also
+        ========
+
+        vector_coplanar
+        is_pure
+
+        """
+        if (self.is_zero_quaternion()) or (other.is_zero_quaternion()):
+            raise ValueError('Neither of the given quaternions can be 0')
+
+        return fuzzy_or([(self.axis() - other.axis()).is_zero_quaternion(), (self.axis() + other.axis()).is_zero_quaternion()])
+
+    @classmethod
+    def vector_coplanar(cls, q1, q2, q3):
+        r"""
+        Returns True if the axis of the pure quaternions seen as 3D vectors
+        ``q1``, ``q2``, and ``q3`` are coplanar.
+
+        Explanation
+        ===========
+
+        Three pure quaternions are vector coplanar if the quaternions seen as 3D vectors are coplanar.
+
+        Parameters
+        ==========
+
+        q1
+            A pure Quaternion.
+        q2
+            A pure Quaternion.
+        q3
+            A pure Quaternion.
+
+        Returns
+        =======
+
+        True : if the axis of the pure quaternions seen as 3D vectors
+        q1, q2, and q3 are coplanar.
+        False : if the axis of the pure quaternions seen as 3D vectors
+        q1, q2, and q3 are not coplanar.
+        None : if the axis of the pure quaternions seen as 3D vectors
+        q1, q2, and q3 are coplanar is unknown.
+
+        Examples
+        ========
+
+        >>> from sympy.algebras.quaternion import Quaternion
+        >>> q1 = Quaternion(0, 4, 4, 4)
+        >>> q2 = Quaternion(0, 8, 8, 8)
+        >>> q3 = Quaternion(0, 24, 24, 24)
+        >>> Quaternion.vector_coplanar(q1, q2, q3)
+        True
+
+        >>> q1 = Quaternion(0, 8, 16, 8)
+        >>> q2 = Quaternion(0, 8, 3, 12)
+        >>> Quaternion.vector_coplanar(q1, q2, q3)
+        False
+
+        See Also
+        ========
+
+        axis
+        is_pure
+
+        """
+
+        if fuzzy_not(q1.is_pure()) or fuzzy_not(q2.is_pure()) or fuzzy_not(q3.is_pure()):
+            raise ValueError('The given quaternions must be pure')
+
+        M = Matrix([[q1.b, q1.c, q1.d], [q2.b, q2.c, q2.d], [q3.b, q3.c, q3.d]]).det()
+        return M.is_zero
+
+    def parallel(self, other):
+        """
+        Returns True if the two pure quaternions seen as 3D vectors are parallel.
+
+        Explanation
+        ===========
+
+        Two pure quaternions are called parallel when their vector product is commutative which
+        implies that the quaternions seen as 3D vectors have same direction.
+
+        Parameters
+        ==========
+
+        other : a Quaternion
+
+        Returns
+        =======
+
+        True : if the two pure quaternions seen as 3D vectors are parallel.
+        False : if the two pure quaternions seen as 3D vectors are not parallel.
+        None : if the two pure quaternions seen as 3D vectors are parallel is unknown.
+
+        Examples
+        ========
+
+        >>> from sympy.algebras.quaternion import Quaternion
+        >>> q = Quaternion(0, 4, 4, 4)
+        >>> q1 = Quaternion(0, 8, 8, 8)
+        >>> q.parallel(q1)
+        True
+
+        >>> q1 = Quaternion(0, 8, 13, 12)
+        >>> q.parallel(q1)
+        False
+
+        """
+
+        if fuzzy_not(self.is_pure()) or fuzzy_not(other.is_pure()):
+            raise ValueError('The provided quaternions must be pure')
+
+        return (self*other - other*self).is_zero_quaternion()
+
+    def orthogonal(self, other):
+        """
+        Returns the orthogonality of two quaternions.
+
+        Explanation
+        ===========
+
+        Two pure quaternions are called orthogonal when their product is anti-commutative.
+
+        Parameters
+        ==========
+
+        other : a Quaternion
+
+        Returns
+        =======
+
+        True : if the two pure quaternions seen as 3D vectors are orthogonal.
+        False : if the two pure quaternions seen as 3D vectors are not orthogonal.
+        None : if the two pure quaternions seen as 3D vectors are orthogonal is unknown.
+
+        Examples
+        ========
+
+        >>> from sympy.algebras.quaternion import Quaternion
+        >>> q = Quaternion(0, 4, 4, 4)
+        >>> q1 = Quaternion(0, 8, 8, 8)
+        >>> q.orthogonal(q1)
+        False
+
+        >>> q1 = Quaternion(0, 2, 2, 0)
+        >>> q = Quaternion(0, 2, -2, 0)
+        >>> q.orthogonal(q1)
+        True
+
+        """
+
+        if fuzzy_not(self.is_pure()) or fuzzy_not(other.is_pure()):
+            raise ValueError('The given quaternions must be pure')
+
+        return (self*other + other*self).is_zero_quaternion()
+
+    def index_vector(self):
+        r"""
+        Returns the index vector of the quaternion.
+
+        Explanation
+        ===========
+
+        The index vector is given by $\mathbf{T}(q)$, the norm (or magnitude) of
+        the quaternion $q$, multiplied by $\mathbf{Ax}(q)$, the axis of $q$.
+
+        Returns
+        =======
+
+        Quaternion: representing index vector of the provided quaternion.
+
+        Examples
+        ========
+
+        >>> from sympy.algebras.quaternion import Quaternion
+        >>> q = Quaternion(2, 4, 2, 4)
+        >>> q.index_vector()
+        0 + 4*sqrt(10)/3*i + 2*sqrt(10)/3*j + 4*sqrt(10)/3*k
+
+        See Also
+        ========
+
+        axis
+        norm
+
+        """
+
+        return self.norm() * self.axis()
+
+    def mensor(self):
+        """
+        Returns the natural logarithm of the norm(magnitude) of the quaternion.
+
+        Examples
+        ========
+
+        >>> from sympy.algebras.quaternion import Quaternion
+        >>> q = Quaternion(2, 4, 2, 4)
+        >>> q.mensor()
+        log(2*sqrt(10))
+        >>> q.norm()
+        2*sqrt(10)
+
+        See Also
+        ========
+
+        norm
+
+        """
+
+        return ln(self.norm())

lib/python3.1/site-packages/sympy/algebras/tests/test_quaternion.py

@@ -0,0 +1,428 @@
+from sympy.testing.pytest import slow
+from sympy.core.function import diff
+from sympy.core.function import expand
+from sympy.core.numbers import (E, I, Rational, pi)
+from sympy.core.singleton import S
+from sympy.core.symbol import (Symbol, symbols)
+from sympy.functions.elementary.complexes import (Abs, conjugate, im, re, sign)
+from sympy.functions.elementary.exponential import log
+from sympy.functions.elementary.miscellaneous import sqrt
+from sympy.functions.elementary.trigonometric import (acos, asin, cos, sin, atan2, atan)
+from sympy.integrals.integrals import integrate
+from sympy.matrices.dense import Matrix
+from sympy.simplify import simplify
+from sympy.simplify.trigsimp import trigsimp
+from sympy.algebras.quaternion import Quaternion
+from sympy.testing.pytest import raises
+import math
+from itertools import permutations, product
+
+w, x, y, z = symbols('w:z')
+phi = symbols('phi')
+
+def test_quaternion_construction():
+    q = Quaternion(w, x, y, z)
+    assert q + q == Quaternion(2*w, 2*x, 2*y, 2*z)
+
+    q2 = Quaternion.from_axis_angle((sqrt(3)/3, sqrt(3)/3, sqrt(3)/3),
+                                    pi*Rational(2, 3))
+    assert q2 == Quaternion(S.Half, S.Half,
+                            S.Half, S.Half)
+
+    M = Matrix([[cos(phi), -sin(phi), 0], [sin(phi), cos(phi), 0], [0, 0, 1]])
+    q3 = trigsimp(Quaternion.from_rotation_matrix(M))
+    assert q3 == Quaternion(
+        sqrt(2)*sqrt(cos(phi) + 1)/2, 0, 0, sqrt(2 - 2*cos(phi))*sign(sin(phi))/2)
+
+    nc = Symbol('nc', commutative=False)
+    raises(ValueError, lambda: Quaternion(w, x, nc, z))
+
+
+def test_quaternion_construction_norm():
+    q1 = Quaternion(*symbols('a:d'))
+
+    q2 = Quaternion(w, x, y, z)
+    assert expand((q1*q2).norm()**2 - (q1.norm()**2 * q2.norm()**2)) == 0
+
+    q3 = Quaternion(w, x, y, z, norm=1)
+    assert (q1 * q3).norm() == q1.norm()
+
+
+def test_issue_25254():
+    # calculating the inverse cached the norm which caused problems
+    # when multiplying
+    p = Quaternion(1, 0, 0, 0)
+    q = Quaternion.from_axis_angle((1, 1, 1), 3 * math.pi/4)
+    qi = q.inverse()  # this operation cached the norm
+    test = q * p * qi
+    assert ((test - p).norm() < 1E-10)
+
+
+def test_to_and_from_Matrix():
+    q = Quaternion(w, x, y, z)
+    q_full = Quaternion.from_Matrix(q.to_Matrix())
+    q_vect = Quaternion.from_Matrix(q.to_Matrix(True))
+    assert (q - q_full).is_zero_quaternion()
+    assert (q.vector_part() - q_vect).is_zero_quaternion()
+
+
+def test_product_matrices():
+    q1 = Quaternion(w, x, y, z)
+    q2 = Quaternion(*(symbols("a:d")))
+    assert (q1 * q2).to_Matrix() == q1.product_matrix_left * q2.to_Matrix()
+    assert (q1 * q2).to_Matrix() == q2.product_matrix_right * q1.to_Matrix()
+
+    R1 = (q1.product_matrix_left * q1.product_matrix_right.T)[1:, 1:]
+    R2 = simplify(q1.to_rotation_matrix()*q1.norm()**2)
+    assert R1 == R2
+
+
+def test_quaternion_axis_angle():
+
+    test_data = [ # axis, angle, expected_quaternion
+        ((1, 0, 0), 0, (1, 0, 0, 0)),
+        ((1, 0, 0), pi/2, (sqrt(2)/2, sqrt(2)/2, 0, 0)),
+        ((0, 1, 0), pi/2, (sqrt(2)/2, 0, sqrt(2)/2, 0)),
+        ((0, 0, 1), pi/2, (sqrt(2)/2, 0, 0, sqrt(2)/2)),
+        ((1, 0, 0), pi, (0, 1, 0, 0)),
+        ((0, 1, 0), pi, (0, 0, 1, 0)),
+        ((0, 0, 1), pi, (0, 0, 0, 1)),
+        ((1, 1, 1), pi, (0, 1/sqrt(3),1/sqrt(3),1/sqrt(3))),
+        ((sqrt(3)/3, sqrt(3)/3, sqrt(3)/3), pi*2/3, (S.Half, S.Half, S.Half, S.Half))
+    ]
+
+    for axis, angle, expected in test_data:
+        assert Quaternion.from_axis_angle(axis, angle) == Quaternion(*expected)
+
+
+def test_quaternion_axis_angle_simplification():
+    result = Quaternion.from_axis_angle((1, 2, 3), asin(4))
+    assert result.a == cos(asin(4)/2)
+    assert result.b == sqrt(14)*sin(asin(4)/2)/14
+    assert result.c == sqrt(14)*sin(asin(4)/2)/7
+    assert result.d == 3*sqrt(14)*sin(asin(4)/2)/14
+
+def test_quaternion_complex_real_addition():
+    a = symbols("a", complex=True)
+    b = symbols("b", real=True)
+    # This symbol is not complex:
+    c = symbols("c", commutative=False)
+
+    q = Quaternion(w, x, y, z)
+    assert a + q == Quaternion(w + re(a), x + im(a), y, z)
+    assert 1 + q == Quaternion(1 + w, x, y, z)
+    assert I + q == Quaternion(w, 1 + x, y, z)
+    assert b + q == Quaternion(w + b, x, y, z)
+    raises(ValueError, lambda: c + q)
+    raises(ValueError, lambda: q * c)
+    raises(ValueError, lambda: c * q)
+
+    assert -q == Quaternion(-w, -x, -y, -z)
+
+    q1 = Quaternion(3 + 4*I, 2 + 5*I, 0, 7 + 8*I, real_field = False)
+    q2 = Quaternion(1, 4, 7, 8)
+
+    assert q1 + (2 + 3*I) == Quaternion(5 + 7*I, 2 + 5*I, 0, 7 + 8*I)
+    assert q2 + (2 + 3*I) == Quaternion(3, 7, 7, 8)
+    assert q1 * (2 + 3*I) == \
+    Quaternion((2 + 3*I)*(3 + 4*I), (2 + 3*I)*(2 + 5*I), 0, (2 + 3*I)*(7 + 8*I))
+    assert q2 * (2 + 3*I) == Quaternion(-10, 11, 38, -5)
+
+    q1 = Quaternion(1, 2, 3, 4)
+    q0 = Quaternion(0, 0, 0, 0)
+    assert q1 + q0 == q1
+    assert q1 - q0 == q1
+    assert q1 - q1 == q0
+
+
+def test_quaternion_subs():
+    q = Quaternion.from_axis_angle((0, 0, 1), phi)
+    assert q.subs(phi, 0) == Quaternion(1, 0, 0, 0)
+
+
+def test_quaternion_evalf():
+    assert (Quaternion(sqrt(2), 0, 0, sqrt(3)).evalf() ==
+            Quaternion(sqrt(2).evalf(), 0, 0, sqrt(3).evalf()))
+    assert (Quaternion(1/sqrt(2), 0, 0, 1/sqrt(2)).evalf() ==
+            Quaternion((1/sqrt(2)).evalf(), 0, 0, (1/sqrt(2)).evalf()))
+
+
+def test_quaternion_functions():
+    q = Quaternion(w, x, y, z)
+    q1 = Quaternion(1, 2, 3, 4)
+    q0 = Quaternion(0, 0, 0, 0)
+
+    assert conjugate(q) == Quaternion(w, -x, -y, -z)
+    assert q.norm() == sqrt(w**2 + x**2 + y**2 + z**2)
+    assert q.normalize() == Quaternion(w, x, y, z) / sqrt(w**2 + x**2 + y**2 + z**2)
+    assert q.inverse() == Quaternion(w, -x, -y, -z) / (w**2 + x**2 + y**2 + z**2)
+    assert q.inverse() == q.pow(-1)
+    raises(ValueError, lambda: q0.inverse())
+    assert q.pow(2) == Quaternion(w**2 - x**2 - y**2 - z**2, 2*w*x, 2*w*y, 2*w*z)
+    assert q**(2) == Quaternion(w**2 - x**2 - y**2 - z**2, 2*w*x, 2*w*y, 2*w*z)
+    assert q1.pow(-2) == Quaternion(
+        Rational(-7, 225), Rational(-1, 225), Rational(-1, 150), Rational(-2, 225))
+    assert q1**(-2) == Quaternion(
+        Rational(-7, 225), Rational(-1, 225), Rational(-1, 150), Rational(-2, 225))
+    assert q1.pow(-0.5) == NotImplemented
+    raises(TypeError, lambda: q1**(-0.5))
+
+    assert q1.exp() == \
+    Quaternion(E * cos(sqrt(29)),
+               2 * sqrt(29) * E * sin(sqrt(29)) / 29,
+               3 * sqrt(29) * E * sin(sqrt(29)) / 29,
+               4 * sqrt(29) * E * sin(sqrt(29)) / 29)
+    assert q1.log() == \
+    Quaternion(log(sqrt(30)),
+               2 * sqrt(29) * acos(sqrt(30)/30) / 29,
+               3 * sqrt(29) * acos(sqrt(30)/30) / 29,
+               4 * sqrt(29) * acos(sqrt(30)/30) / 29)
+
+    assert q1.pow_cos_sin(2) == \
+    Quaternion(30 * cos(2 * acos(sqrt(30)/30)),
+               60 * sqrt(29) * sin(2 * acos(sqrt(30)/30)) / 29,
+               90 * sqrt(29) * sin(2 * acos(sqrt(30)/30)) / 29,
+               120 * sqrt(29) * sin(2 * acos(sqrt(30)/30)) / 29)
+
+    assert diff(Quaternion(x, x, x, x), x) == Quaternion(1, 1, 1, 1)
+
+    assert integrate(Quaternion(x, x, x, x), x) == \
+    Quaternion(x**2 / 2, x**2 / 2, x**2 / 2, x**2 / 2)
+
+    assert Quaternion.rotate_point((1, 1, 1), q1) == (S.One / 5, 1, S(7) / 5)
+    n = Symbol('n')
+    raises(TypeError, lambda: q1**n)
+    n = Symbol('n', integer=True)
+    raises(TypeError, lambda: q1**n)
+
+    assert Quaternion(22, 23, 55, 8).scalar_part() == 22
+    assert Quaternion(w, x, y, z).scalar_part() == w
+
+    assert Quaternion(22, 23, 55, 8).vector_part() == Quaternion(0, 23, 55, 8)
+    assert Quaternion(w, x, y, z).vector_part() == Quaternion(0, x, y, z)
+
+    assert q1.axis() == Quaternion(0, 2*sqrt(29)/29, 3*sqrt(29)/29, 4*sqrt(29)/29)
+    assert q1.axis().pow(2) == Quaternion(-1, 0, 0, 0)
+    assert q0.axis().scalar_part() == 0
+    assert (q.axis() == Quaternion(0,
+                                   x/sqrt(x**2 + y**2 + z**2),
+                                   y/sqrt(x**2 + y**2 + z**2),
+                                   z/sqrt(x**2 + y**2 + z**2)))
+
+    assert q0.is_pure() is True
+    assert q1.is_pure() is False
+    assert Quaternion(0, 0, 0, 3).is_pure() is True
+    assert Quaternion(0, 2, 10, 3).is_pure() is True
+    assert Quaternion(w, 2, 10, 3).is_pure() is None
+
+    assert q1.angle() == 2*atan(sqrt(29))
+    assert q.angle() == 2*atan2(sqrt(x**2 + y**2 + z**2), w)
+
+    assert Quaternion.arc_coplanar(q1, Quaternion(2, 4, 6, 8)) is True
+    assert Quaternion.arc_coplanar(q1, Quaternion(1, -2, -3, -4)) is True
+    assert Quaternion.arc_coplanar(q1, Quaternion(1, 8, 12, 16)) is True
+    assert Quaternion.arc_coplanar(q1, Quaternion(1, 2, 3, 4)) is True
+    assert Quaternion.arc_coplanar(q1, Quaternion(w, 4, 6, 8)) is True
+    assert Quaternion.arc_coplanar(q1, Quaternion(2, 7, 4, 1)) is False
+    assert Quaternion.arc_coplanar(q1, Quaternion(w, x, y, z)) is None
+    raises(ValueError, lambda: Quaternion.arc_coplanar(q1, q0))
+
+    assert Quaternion.vector_coplanar(
+        Quaternion(0, 8, 12, 16),
+        Quaternion(0, 4, 6, 8),
+        Quaternion(0, 2, 3, 4)) is True
+    assert Quaternion.vector_coplanar(
+        Quaternion(0, 0, 0, 0), Quaternion(0, 4, 6, 8), Quaternion(0, 2, 3, 4)) is True
+    assert Quaternion.vector_coplanar(
+        Quaternion(0, 8, 2, 6), Quaternion(0, 1, 6, 6), Quaternion(0, 0, 3, 4)) is False
+    assert Quaternion.vector_coplanar(
+        Quaternion(0, 1, 3, 4),
+        Quaternion(0, 4, w, 6),
+        Quaternion(0, 6, 8, 1)) is None
+    raises(ValueError, lambda:
+        Quaternion.vector_coplanar(q0, Quaternion(0, 4, 6, 8), q1))
+
+    assert Quaternion(0, 1, 2, 3).parallel(Quaternion(0, 2, 4, 6)) is True
+    assert Quaternion(0, 1, 2, 3).parallel(Quaternion(0, 2, 2, 6)) is False
+    assert Quaternion(0, 1, 2, 3).parallel(Quaternion(w, x, y, 6)) is None
+    raises(ValueError, lambda: q0.parallel(q1))
+
+    assert Quaternion(0, 1, 2, 3).orthogonal(Quaternion(0, -2, 1, 0)) is True
+    assert Quaternion(0, 2, 4, 7).orthogonal(Quaternion(0, 2, 2, 6)) is False
+    assert Quaternion(0, 2, 4, 7).orthogonal(Quaternion(w, x, y, 6)) is None
+    raises(ValueError, lambda: q0.orthogonal(q1))
+
+    assert q1.index_vector() == Quaternion(
+        0, 2*sqrt(870)/29,
+        3*sqrt(870)/29,
+        4*sqrt(870)/29)
+    assert Quaternion(0, 3, 9, 4).index_vector() == Quaternion(0, 3, 9, 4)
+
+    assert Quaternion(4, 3, 9, 4).mensor() == log(sqrt(122))
+    assert Quaternion(3, 3, 0, 2).mensor() == log(sqrt(22))
+
+    assert q0.is_zero_quaternion() is True
+    assert q1.is_zero_quaternion() is False
+    assert Quaternion(w, 0, 0, 0).is_zero_quaternion() is None
+
+def test_quaternion_conversions():
+    q1 = Quaternion(1, 2, 3, 4)
+
+    assert q1.to_axis_angle() == ((2 * sqrt(29)/29,
+                                   3 * sqrt(29)/29,
+                                   4 * sqrt(29)/29),
+                                   2 * acos(sqrt(30)/30))
+
+    assert (q1.to_rotation_matrix() ==
+            Matrix([[Rational(-2, 3), Rational(2, 15), Rational(11, 15)],
+                    [Rational(2, 3), Rational(-1, 3), Rational(2, 3)],
+                    [Rational(1, 3), Rational(14, 15), Rational(2, 15)]]))
+
+    assert (q1.to_rotation_matrix((1, 1, 1)) ==
+            Matrix([
+                [Rational(-2, 3), Rational(2, 15), Rational(11, 15), Rational(4, 5)],
+                [Rational(2, 3), Rational(-1, 3), Rational(2, 3), S.Zero],
+                [Rational(1, 3), Rational(14, 15), Rational(2, 15), Rational(-2, 5)],
+                [S.Zero, S.Zero, S.Zero, S.One]]))
+
+    theta = symbols("theta", real=True)
+    q2 = Quaternion(cos(theta/2), 0, 0, sin(theta/2))
+
+    assert trigsimp(q2.to_rotation_matrix()) == Matrix([
+                                               [cos(theta), -sin(theta), 0],
+                                               [sin(theta),  cos(theta), 0],
+                                               [0,           0,          1]])
+
+    assert q2.to_axis_angle() == ((0, 0, sin(theta/2)/Abs(sin(theta/2))),
+                                   2*acos(cos(theta/2)))
+
+    assert trigsimp(q2.to_rotation_matrix((1, 1, 1))) == Matrix([
+               [cos(theta), -sin(theta), 0, sin(theta) - cos(theta) + 1],
+               [sin(theta),  cos(theta), 0, -sin(theta) - cos(theta) + 1],
+               [0,           0,          1,  0],
+               [0,           0,          0,  1]])
+
+
+def test_rotation_matrix_homogeneous():
+    q = Quaternion(w, x, y, z)
+    R1 = q.to_rotation_matrix(homogeneous=True) * q.norm()**2
+    R2 = simplify(q.to_rotation_matrix(homogeneous=False) * q.norm()**2)
+    assert R1 == R2
+
+
+def test_quaternion_rotation_iss1593():
+    """
+    There was a sign mistake in the definition,
+    of the rotation matrix. This tests that particular sign mistake.
+    See issue 1593 for reference.
+    See wikipedia
+    https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Quaternion-derived_rotation_matrix
+    for the correct definition
+    """
+    q = Quaternion(cos(phi/2), sin(phi/2), 0, 0)
+    assert(trigsimp(q.to_rotation_matrix()) == Matrix([
+                [1,        0,         0],
+                [0, cos(phi), -sin(phi)],
+                [0, sin(phi),  cos(phi)]]))
+
+
+def test_quaternion_multiplication():
+    q1 = Quaternion(3 + 4*I, 2 + 5*I, 0, 7 + 8*I, real_field = False)
+    q2 = Quaternion(1, 2, 3, 5)
+    q3 = Quaternion(1, 1, 1, y)
+
+    assert Quaternion._generic_mul(S(4), S.One) == 4
+    assert (Quaternion._generic_mul(S(4), q1) ==
+            Quaternion(12 + 16*I, 8 + 20*I, 0, 28 + 32*I))
+    assert q2.mul(2) == Quaternion(2, 4, 6, 10)
+    assert q2.mul(q3) == Quaternion(-5*y - 4, 3*y - 2, 9 - 2*y, y + 4)
+    assert q2.mul(q3) == q2*q3
+
+    z = symbols('z', complex=True)
+    z_quat = Quaternion(re(z), im(z), 0, 0)
+    q = Quaternion(*symbols('q:4', real=True))
+
+    assert z * q == z_quat * q
+    assert q * z == q * z_quat
+
+
+def test_issue_16318():
+    #for rtruediv
+    q0 = Quaternion(0, 0, 0, 0)
+    raises(ValueError, lambda: 1/q0)
+    #for rotate_point
+    q = Quaternion(1, 2, 3, 4)
+    (axis, angle) = q.to_axis_angle()
+    assert Quaternion.rotate_point((1, 1, 1), (axis, angle)) == (S.One / 5, 1, S(7) / 5)
+    #test for to_axis_angle
+    q = Quaternion(-1, 1, 1, 1)
+    axis = (-sqrt(3)/3, -sqrt(3)/3, -sqrt(3)/3)
+    angle = 2*pi/3
+    assert (axis, angle) == q.to_axis_angle()
+
+
+@slow
+def test_to_euler():
+    q = Quaternion(w, x, y, z)
+    q_normalized = q.normalize()
+
+    seqs = ['zxy', 'zyx', 'zyz', 'zxz']
+    seqs += [seq.upper() for seq in seqs]
+
+    for seq in seqs:
+        euler_from_q = q.to_euler(seq)
+        q_back = simplify(Quaternion.from_euler(euler_from_q, seq))
+        assert q_back == q_normalized
+
+
+def test_to_euler_iss24504():
+    """
+    There was a mistake in the degenerate case testing
+    See issue 24504 for reference.
+    """
+    q = Quaternion.from_euler((phi, 0, 0), 'zyz')
+    assert trigsimp(q.to_euler('zyz'), inverse=True) == (phi, 0, 0)
+
+
+def test_to_euler_numerical_singilarities():
+
+    def test_one_case(angles, seq):
+        q = Quaternion.from_euler(angles, seq)
+        assert q.to_euler(seq) == angles
+
+    # symmetric
+    test_one_case((pi/2,  0, 0), 'zyz')
+    test_one_case((pi/2,  0, 0), 'ZYZ')
+    test_one_case((pi/2,  pi, 0), 'zyz')
+    test_one_case((pi/2,  pi, 0), 'ZYZ')
+
+    # asymmetric
+    test_one_case((pi/2,  pi/2, 0), 'zyx')
+    test_one_case((pi/2,  -pi/2, 0), 'zyx')
+    test_one_case((pi/2,  pi/2, 0), 'ZYX')
+    test_one_case((pi/2,  -pi/2, 0), 'ZYX')
+
+
+@slow
+def test_to_euler_options():
+    def test_one_case(q):
+        angles1 = Matrix(q.to_euler(seq, True, True))
+        angles2 = Matrix(q.to_euler(seq, False, False))
+        angle_errors = simplify(angles1-angles2).evalf()
+        for angle_error in angle_errors:
+            # forcing angles to set {-pi, pi}
+            angle_error = (angle_error + pi) % (2 * pi) - pi
+            assert angle_error < 10e-7
+
+    for xyz in ('xyz', 'XYZ'):
+        for seq_tuple in permutations(xyz):
+            for symmetric in (True, False):
+                if symmetric:
+                    seq = ''.join([seq_tuple[0], seq_tuple[1], seq_tuple[0]])
+                else:
+                    seq = ''.join(seq_tuple)
+
+                for elements in product([-1, 0, 1], repeat=4):
+                    q = Quaternion(*elements)
+                    if not q.is_zero_quaternion():
+                        test_one_case(q)

lib/python3.1/site-packages/sympy/codegen/__init__.py

@@ -0,0 +1,24 @@
+""" The ``sympy.codegen`` module contains classes and functions for building
+abstract syntax trees of algorithms. These trees may then be printed by the
+code-printers in ``sympy.printing``.
+
+There are several submodules available:
+- ``sympy.codegen.ast``: AST nodes useful across multiple languages.
+- ``sympy.codegen.cnodes``: AST nodes useful for the C family of languages.
+- ``sympy.codegen.fnodes``: AST nodes useful for Fortran.
+- ``sympy.codegen.cfunctions``: functions specific to C (C99 math functions)
+- ``sympy.codegen.ffunctions``: functions specific to Fortran (e.g. ``kind``).
+
+
+
+"""
+from .ast import (
+    Assignment, aug_assign, CodeBlock, For, Attribute, Variable, Declaration,
+    While, Scope, Print, FunctionPrototype, FunctionDefinition, FunctionCall
+)
+
+__all__ = [
+    'Assignment', 'aug_assign', 'CodeBlock', 'For', 'Attribute', 'Variable',
+    'Declaration', 'While', 'Scope', 'Print', 'FunctionPrototype',
+    'FunctionDefinition', 'FunctionCall',
+]

lib/python3.1/site-packages/sympy/codegen/abstract_nodes.py

@@ -0,0 +1,18 @@
+"""This module provides containers for python objects that are valid
+printing targets but are not a subclass of SymPy's Printable.
+"""
+
+
+from sympy.core.containers import Tuple
+
+
+class List(Tuple):
+    """Represents a (frozen) (Python) list (for code printing purposes)."""
+    def __eq__(self, other):
+        if isinstance(other, list):
+            return self == List(*other)
+        else:
+            return self.args == other
+
+    def __hash__(self):
+        return super().__hash__()

lib/python3.1/site-packages/sympy/codegen/algorithms.py

@@ -0,0 +1,180 @@
+from sympy.core.containers import Tuple
+from sympy.core.numbers import oo
+from sympy.core.relational import (Gt, Lt)
+from sympy.core.symbol import (Dummy, Symbol)
+from sympy.functions.elementary.complexes import Abs
+from sympy.functions.elementary.miscellaneous import Min, Max
+from sympy.logic.boolalg import And
+from sympy.codegen.ast import (
+    Assignment, AddAugmentedAssignment, break_, CodeBlock, Declaration, FunctionDefinition,
+    Print, Return, Scope, While, Variable, Pointer, real
+)
+from sympy.codegen.cfunctions import isnan
+
+""" This module collects functions for constructing ASTs representing algorithms. """
+
+def newtons_method(expr, wrt, atol=1e-12, delta=None, *, rtol=4e-16, debug=False,
+                   itermax=None, counter=None, delta_fn=lambda e, x: -e/e.diff(x),
+                   cse=False, handle_nan=None,
+                   bounds=None):
+    """ Generates an AST for Newton-Raphson method (a root-finding algorithm).
+
+    Explanation
+    ===========
+
+    Returns an abstract syntax tree (AST) based on ``sympy.codegen.ast`` for Netwon's
+    method of root-finding.
+
+    Parameters
+    ==========
+
+    expr : expression
+    wrt : Symbol
+        With respect to, i.e. what is the variable.
+    atol : number or expression
+        Absolute tolerance (stopping criterion)
+    rtol : number or expression
+        Relative tolerance (stopping criterion)
+    delta : Symbol
+        Will be a ``Dummy`` if ``None``.
+    debug : bool
+        Whether to print convergence information during iterations
+    itermax : number or expr
+        Maximum number of iterations.
+    counter : Symbol
+        Will be a ``Dummy`` if ``None``.
+    delta_fn: Callable[[Expr, Symbol], Expr]
+        computes the step, default is newtons method. For e.g. Halley's method
+        use delta_fn=lambda e, x: -2*e*e.diff(x)/(2*e.diff(x)**2 - e*e.diff(x, 2))
+    cse: bool
+        Perform common sub-expression elimination on delta expression
+    handle_nan: Token
+        How to handle occurrence of not-a-number (NaN).
+    bounds: Optional[tuple[Expr, Expr]]
+        Perform optimization within bounds
+
+    Examples
+    ========
+
+    >>> from sympy import symbols, cos
+    >>> from sympy.codegen.ast import Assignment
+    >>> from sympy.codegen.algorithms import newtons_method
+    >>> x, dx, atol = symbols('x dx atol')
+    >>> expr = cos(x) - x**3
+    >>> algo = newtons_method(expr, x, atol=atol, delta=dx)
+    >>> algo.has(Assignment(dx, -expr/expr.diff(x)))
+    True
+
+    References
+    ==========
+
+    .. [1] https://en.wikipedia.org/wiki/Newton%27s_method
+
+    """
+
+    if delta is None:
+        delta = Dummy()
+        Wrapper = Scope
+        name_d = 'delta'
+    else:
+        Wrapper = lambda x: x
+        name_d = delta.name
+
+    delta_expr = delta_fn(expr, wrt)
+    if cse:
+        from sympy.simplify.cse_main import cse
+        cses, (red,) = cse([delta_expr.factor()])
+        whl_bdy = [Assignment(dum, sub_e) for dum, sub_e in cses]
+        whl_bdy += [Assignment(delta, red)]
+    else:
+        whl_bdy = [Assignment(delta, delta_expr)]
+    if handle_nan is not None:
+        whl_bdy += [While(isnan(delta), CodeBlock(handle_nan, break_))]
+    whl_bdy += [AddAugmentedAssignment(wrt, delta)]
+    if bounds is not None:
+        whl_bdy += [Assignment(wrt, Min(Max(wrt, bounds[0]), bounds[1]))]
+    if debug:
+        prnt = Print([wrt, delta], r"{}=%12.5g {}=%12.5g\n".format(wrt.name, name_d))
+        whl_bdy += [prnt]
+    req = Gt(Abs(delta), atol + rtol*Abs(wrt))
+    declars = [Declaration(Variable(delta, type=real, value=oo))]
+    if itermax is not None:
+        counter = counter or Dummy(integer=True)
+        v_counter = Variable.deduced(counter, 0)
+        declars.append(Declaration(v_counter))
+        whl_bdy.append(AddAugmentedAssignment(counter, 1))
+        req = And(req, Lt(counter, itermax))
+    whl = While(req, CodeBlock(*whl_bdy))
+    blck = declars
+    if debug:
+        blck.append(Print([wrt], r"{}=%12.5g\n".format(wrt.name)))
+    blck += [whl]
+    return Wrapper(CodeBlock(*blck))
+
+
+def _symbol_of(arg):
+    if isinstance(arg, Declaration):
+        arg = arg.variable.symbol
+    elif isinstance(arg, Variable):
+        arg = arg.symbol
+    return arg
+
+
+def newtons_method_function(expr, wrt, params=None, func_name="newton", attrs=Tuple(), *, delta=None, **kwargs):
+    """ Generates an AST for a function implementing the Newton-Raphson method.
+
+    Parameters
+    ==========
+
+    expr : expression
+    wrt : Symbol
+        With respect to, i.e. what is the variable
+    params : iterable of symbols
+        Symbols appearing in expr that are taken as constants during the iterations
+        (these will be accepted as parameters to the generated function).
+    func_name : str
+        Name of the generated function.
+    attrs : Tuple
+        Attribute instances passed as ``attrs`` to ``FunctionDefinition``.
+    \\*\\*kwargs :
+        Keyword arguments passed to :func:`sympy.codegen.algorithms.newtons_method`.
+
+    Examples
+    ========
+
+    >>> from sympy import symbols, cos
+    >>> from sympy.codegen.algorithms import newtons_method_function
+    >>> from sympy.codegen.pyutils import render_as_module
+    >>> x = symbols('x')
+    >>> expr = cos(x) - x**3
+    >>> func = newtons_method_function(expr, x)
+    >>> py_mod = render_as_module(func)  # source code as string
+    >>> namespace = {}
+    >>> exec(py_mod, namespace, namespace)
+    >>> res = eval('newton(0.5)', namespace)
+    >>> abs(res - 0.865474033102) < 1e-12
+    True
+
+    See Also
+    ========
+
+    sympy.codegen.algorithms.newtons_method
+
+    """
+    if params is None:
+        params = (wrt,)
+    pointer_subs = {p.symbol: Symbol('(*%s)' % p.symbol.name)
+                    for p in params if isinstance(p, Pointer)}
+    if delta is None:
+        delta = Symbol('d_' + wrt.name)
+        if expr.has(delta):
+            delta = None  # will use Dummy
+    algo = newtons_method(expr, wrt, delta=delta, **kwargs).xreplace(pointer_subs)
+    if isinstance(algo, Scope):
+        algo = algo.body
+    not_in_params = expr.free_symbols.difference({_symbol_of(p) for p in params})
+    if not_in_params:
+        raise ValueError("Missing symbols in params: %s" % ', '.join(map(str, not_in_params)))
+    declars = tuple(Variable(p, real) for p in params)
+    body = CodeBlock(algo, Return(wrt))
+    return FunctionDefinition(real, func_name, declars, body, attrs=attrs)

lib/python3.1/site-packages/sympy/codegen/approximations.py

@@ -0,0 +1,187 @@
+import math
+from sympy.sets.sets import Interval
+from sympy.calculus.singularities import is_increasing, is_decreasing
+from sympy.codegen.rewriting import Optimization
+from sympy.core.function import UndefinedFunction
+
+"""
+This module collects classes useful for approimate rewriting of expressions.
+This can be beneficial when generating numeric code for which performance is
+of greater importance than precision (e.g. for preconditioners used in iterative
+methods).
+"""
+
+class SumApprox(Optimization):
+    """
+    Approximates sum by neglecting small terms.
+
+    Explanation
+    ===========
+
+    If terms are expressions which can be determined to be monotonic, then
+    bounds for those expressions are added.
+
+    Parameters
+    ==========
+
+    bounds : dict
+        Mapping expressions to length 2 tuple of bounds (low, high).
+    reltol : number
+        Threshold for when to ignore a term. Taken relative to the largest
+        lower bound among bounds.
+
+    Examples
+    ========
+
+    >>> from sympy import exp
+    >>> from sympy.abc import x, y, z
+    >>> from sympy.codegen.rewriting import optimize
+    >>> from sympy.codegen.approximations import SumApprox
+    >>> bounds = {x: (-1, 1), y: (1000, 2000), z: (-10, 3)}
+    >>> sum_approx3 = SumApprox(bounds, reltol=1e-3)
+    >>> sum_approx2 = SumApprox(bounds, reltol=1e-2)
+    >>> sum_approx1 = SumApprox(bounds, reltol=1e-1)
+    >>> expr = 3*(x + y + exp(z))
+    >>> optimize(expr, [sum_approx3])
+    3*(x + y + exp(z))
+    >>> optimize(expr, [sum_approx2])
+    3*y + 3*exp(z)
+    >>> optimize(expr, [sum_approx1])
+    3*y
+
+    """
+
+    def __init__(self, bounds, reltol, **kwargs):
+        super().__init__(**kwargs)
+        self.bounds = bounds
+        self.reltol = reltol
+
+    def __call__(self, expr):
+        return expr.factor().replace(self.query, lambda arg: self.value(arg))
+
+    def query(self, expr):
+        return expr.is_Add
+
+    def value(self, add):
+        for term in add.args:
+            if term.is_number or term in self.bounds or len(term.free_symbols) != 1:
+                continue
+            fs, = term.free_symbols
+            if fs not in self.bounds:
+                continue
+            intrvl = Interval(*self.bounds[fs])
+            if is_increasing(term, intrvl, fs):
+                self.bounds[term] = (
+                    term.subs({fs: self.bounds[fs][0]}),
+                    term.subs({fs: self.bounds[fs][1]})
+                )
+            elif is_decreasing(term, intrvl, fs):
+                self.bounds[term] = (
+                    term.subs({fs: self.bounds[fs][1]}),
+                    term.subs({fs: self.bounds[fs][0]})
+                )
+            else:
+                return add
+
+        if all(term.is_number or term in self.bounds for term in add.args):
+            bounds = [(term, term) if term.is_number else self.bounds[term] for term in add.args]
+            largest_abs_guarantee = 0
+            for lo, hi in bounds:
+                if lo <= 0 <= hi:
+                    continue
+                largest_abs_guarantee = max(largest_abs_guarantee,
+                                            min(abs(lo), abs(hi)))
+            new_terms = []
+            for term, (lo, hi) in zip(add.args, bounds):
+                if max(abs(lo), abs(hi)) >= largest_abs_guarantee*self.reltol:
+                    new_terms.append(term)
+            return add.func(*new_terms)
+        else:
+            return add
+
+
+class SeriesApprox(Optimization):
+    """ Approximates functions by expanding them as a series.
+
+    Parameters
+    ==========
+
+    bounds : dict
+        Mapping expressions to length 2 tuple of bounds (low, high).
+    reltol : number
+        Threshold for when to ignore a term. Taken relative to the largest
+        lower bound among bounds.
+    max_order : int
+        Largest order to include in series expansion
+    n_point_checks : int (even)
+        The validity of an expansion (with respect to reltol) is checked at
+        discrete points (linearly spaced over the bounds of the variable). The
+        number of points used in this numerical check is given by this number.
+
+    Examples
+    ========
+
+    >>> from sympy import sin, pi
+    >>> from sympy.abc import x, y
+    >>> from sympy.codegen.rewriting import optimize
+    >>> from sympy.codegen.approximations import SeriesApprox
+    >>> bounds = {x: (-.1, .1), y: (pi-1, pi+1)}
+    >>> series_approx2 = SeriesApprox(bounds, reltol=1e-2)
+    >>> series_approx3 = SeriesApprox(bounds, reltol=1e-3)
+    >>> series_approx8 = SeriesApprox(bounds, reltol=1e-8)
+    >>> expr = sin(x)*sin(y)
+    >>> optimize(expr, [series_approx2])
+    x*(-y + (y - pi)**3/6 + pi)
+    >>> optimize(expr, [series_approx3])
+    (-x**3/6 + x)*sin(y)
+    >>> optimize(expr, [series_approx8])
+    sin(x)*sin(y)
+
+    """
+    def __init__(self, bounds, reltol, max_order=4, n_point_checks=4, **kwargs):
+        super().__init__(**kwargs)
+        self.bounds = bounds
+        self.reltol = reltol
+        self.max_order = max_order
+        if n_point_checks % 2 == 1:
+            raise ValueError("Checking the solution at expansion point is not helpful")
+        self.n_point_checks = n_point_checks
+        self._prec = math.ceil(-math.log10(self.reltol))
+
+    def __call__(self, expr):
+        return expr.factor().replace(self.query, lambda arg: self.value(arg))
+
+    def query(self, expr):
+        return (expr.is_Function and not isinstance(expr, UndefinedFunction)
+                and len(expr.args) == 1)
+
+    def value(self, fexpr):
+        free_symbols = fexpr.free_symbols
+        if len(free_symbols) != 1:
+            return fexpr
+        symb, = free_symbols
+        if symb not in self.bounds:
+            return fexpr
+        lo, hi = self.bounds[symb]
+        x0 = (lo + hi)/2
+        cheapest = None
+        for n in range(self.max_order+1, 0, -1):
+            fseri = fexpr.series(symb, x0=x0, n=n).removeO()
+            n_ok = True
+            for idx in range(self.n_point_checks):
+                x = lo + idx*(hi - lo)/(self.n_point_checks - 1)
+                val = fseri.xreplace({symb: x})
+                ref = fexpr.xreplace({symb: x})
+                if abs((1 - val/ref).evalf(self._prec)) > self.reltol:
+                    n_ok = False
+                    break
+
+            if n_ok:
+                cheapest = fseri
+            else:
+                break
+
+        if cheapest is None:
+            return fexpr
+        else:
+            return cheapest

lib/python3.1/site-packages/sympy/codegen/ast.py

@@ -0,0 +1,1906 @@
+"""
+Types used to represent a full function/module as an Abstract Syntax Tree.
+
+Most types are small, and are merely used as tokens in the AST. A tree diagram
+has been included below to illustrate the relationships between the AST types.
+
+
+AST Type Tree
+-------------
+::
+
+  *Basic*
+       |
+       |
+   CodegenAST
+       |
+       |--->AssignmentBase
+       |             |--->Assignment
+       |             |--->AugmentedAssignment
+       |                                    |--->AddAugmentedAssignment
+       |                                    |--->SubAugmentedAssignment
+       |                                    |--->MulAugmentedAssignment
+       |                                    |--->DivAugmentedAssignment
+       |                                    |--->ModAugmentedAssignment
+       |
+       |--->CodeBlock
+       |
+       |
+       |--->Token
+                |--->Attribute
+                |--->For
+                |--->String
+                |       |--->QuotedString
+                |       |--->Comment
+                |--->Type
+                |       |--->IntBaseType
+                |       |              |--->_SizedIntType
+                |       |                               |--->SignedIntType
+                |       |                               |--->UnsignedIntType
+                |       |--->FloatBaseType
+                |                        |--->FloatType
+                |                        |--->ComplexBaseType
+                |                                           |--->ComplexType
+                |--->Node
+                |       |--->Variable
+                |       |           |---> Pointer
+                |       |--->FunctionPrototype
+                |                            |--->FunctionDefinition
+                |--->Element
+                |--->Declaration
+                |--->While
+                |--->Scope
+                |--->Stream
+                |--->Print
+                |--->FunctionCall
+                |--->BreakToken
+                |--->ContinueToken
+                |--->NoneToken
+                |--->Return
+
+
+Predefined types
+----------------
+
+A number of ``Type`` instances are provided in the ``sympy.codegen.ast`` module
+for convenience. Perhaps the two most common ones for code-generation (of numeric
+codes) are ``float32`` and ``float64`` (known as single and double precision respectively).
+There are also precision generic versions of Types (for which the codeprinters selects the
+underlying data type at time of printing): ``real``, ``integer``, ``complex_``, ``bool_``.
+
+The other ``Type`` instances defined are:
+
+- ``intc``: Integer type used by C's "int".
+- ``intp``: Integer type used by C's "unsigned".
+- ``int8``, ``int16``, ``int32``, ``int64``: n-bit integers.
+- ``uint8``, ``uint16``, ``uint32``, ``uint64``: n-bit unsigned integers.
+- ``float80``: known as "extended precision" on modern x86/amd64 hardware.
+- ``complex64``: Complex number represented by two ``float32`` numbers
+- ``complex128``: Complex number represented by two ``float64`` numbers
+
+Using the nodes
+---------------
+
+It is possible to construct simple algorithms using the AST nodes. Let's construct a loop applying
+Newton's method::
+
+    >>> from sympy import symbols, cos
+    >>> from sympy.codegen.ast import While, Assignment, aug_assign, Print, QuotedString
+    >>> t, dx, x = symbols('tol delta val')
+    >>> expr = cos(x) - x**3
+    >>> whl = While(abs(dx) > t, [
+    ...     Assignment(dx, -expr/expr.diff(x)),
+    ...     aug_assign(x, '+', dx),
+    ...     Print([x])
+    ... ])
+    >>> from sympy import pycode
+    >>> py_str = pycode(whl)
+    >>> print(py_str)
+    while (abs(delta) > tol):
+        delta = (val**3 - math.cos(val))/(-3*val**2 - math.sin(val))
+        val += delta
+        print(val)
+    >>> import math
+    >>> tol, val, delta = 1e-5, 0.5, float('inf')
+    >>> exec(py_str)
+    1.1121416371
+    0.909672693737
+    0.867263818209
+    0.865477135298
+    0.865474033111
+    >>> print('%3.1g' % (math.cos(val) - val**3))
+    -3e-11
+
+If we want to generate Fortran code for the same while loop we simple call ``fcode``::
+
+    >>> from sympy import fcode
+    >>> print(fcode(whl, standard=2003, source_format='free'))
+    do while (abs(delta) > tol)
+       delta = (val**3 - cos(val))/(-3*val**2 - sin(val))
+       val = val + delta
+       print *, val
+    end do
+
+There is a function constructing a loop (or a complete function) like this in
+:mod:`sympy.codegen.algorithms`.
+
+"""
+
+from __future__ import annotations
+from typing import Any
+
+from collections import defaultdict
+
+from sympy.core.relational import (Ge, Gt, Le, Lt)
+from sympy.core import Symbol, Tuple, Dummy
+from sympy.core.basic import Basic
+from sympy.core.expr import Expr, Atom
+from sympy.core.numbers import Float, Integer, oo
+from sympy.core.sympify import _sympify, sympify, SympifyError
+from sympy.utilities.iterables import (iterable, topological_sort,
+                                       numbered_symbols, filter_symbols)
+
+
+def _mk_Tuple(args):
+    """
+    Create a SymPy Tuple object from an iterable, converting Python strings to
+    AST strings.
+
+    Parameters
+    ==========
+
+    args: iterable
+        Arguments to :class:`sympy.Tuple`.
+
+    Returns
+    =======
+
+    sympy.Tuple
+    """
+    args = [String(arg) if isinstance(arg, str) else arg for arg in args]
+    return Tuple(*args)
+
+
+class CodegenAST(Basic):
+    __slots__ = ()
+
+
+class Token(CodegenAST):
+    """ Base class for the AST types.
+
+    Explanation
+    ===========
+
+    Defining fields are set in ``_fields``. Attributes (defined in _fields)
+    are only allowed to contain instances of Basic (unless atomic, see
+    ``String``). The arguments to ``__new__()`` correspond to the attributes in
+    the order defined in ``_fields`. The ``defaults`` class attribute is a
+    dictionary mapping attribute names to their default values.
+
+    Subclasses should not need to override the ``__new__()`` method. They may
+    define a class or static method named ``_construct_<attr>`` for each
+    attribute to process the value passed to ``__new__()``. Attributes listed
+    in the class attribute ``not_in_args`` are not passed to :class:`~.Basic`.
+    """
+
+    __slots__: tuple[str, ...] = ()
+    _fields = __slots__
+    defaults: dict[str, Any] = {}
+    not_in_args: list[str] = []
+    indented_args = ['body']
+
+    @property
+    def is_Atom(self):
+        return len(self._fields) == 0
+
+    @classmethod
+    def _get_constructor(cls, attr):
+        """ Get the constructor function for an attribute by name. """
+        return getattr(cls, '_construct_%s' % attr, lambda x: x)
+
+    @classmethod
+    def _construct(cls, attr, arg):
+        """ Construct an attribute value from argument passed to ``__new__()``. """
+        # arg may be ``NoneToken()``, so comparison is done using == instead of ``is`` operator
+        if arg == None:
+            return cls.defaults.get(attr, none)
+        else:
+            if isinstance(arg, Dummy):  # SymPy's replace uses Dummy instances
+                return arg
+            else:
+                return cls._get_constructor(attr)(arg)
+
+    def __new__(cls, *args, **kwargs):
+        # Pass through existing instances when given as sole argument
+        if len(args) == 1 and not kwargs and isinstance(args[0], cls):
+            return args[0]
+
+        if len(args) > len(cls._fields):
+            raise ValueError("Too many arguments (%d), expected at most %d" % (len(args), len(cls._fields)))
+
+        attrvals = []
+
+        # Process positional arguments
+        for attrname, argval in zip(cls._fields, args):
+            if attrname in kwargs:
+                raise TypeError('Got multiple values for attribute %r' % attrname)
+
+            attrvals.append(cls._construct(attrname, argval))
+
+        # Process keyword arguments
+        for attrname in cls._fields[len(args):]:
+            if attrname in kwargs:
+                argval = kwargs.pop(attrname)
+
+            elif attrname in cls.defaults:
+                argval = cls.defaults[attrname]
+
+            else:
+                raise TypeError('No value for %r given and attribute has no default' % attrname)
+
+            attrvals.append(cls._construct(attrname, argval))
+
+        if kwargs:
+            raise ValueError("Unknown keyword arguments: %s" % ' '.join(kwargs))
+
+        # Parent constructor
+        basic_args = [
+            val for attr, val in zip(cls._fields, attrvals)
+            if attr not in cls.not_in_args
+        ]
+        obj = CodegenAST.__new__(cls, *basic_args)
+
+        # Set attributes
+        for attr, arg in zip(cls._fields, attrvals):
+            setattr(obj, attr, arg)
+
+        return obj
+
+    def __eq__(self, other):
+        if not isinstance(other, self.__class__):
+            return False
+        for attr in self._fields:
+            if getattr(self, attr) != getattr(other, attr):
+                return False
+        return True
+
+    def _hashable_content(self):
+        return tuple([getattr(self, attr) for attr in self._fields])
+
+    def __hash__(self):
+        return super().__hash__()
+
+    def _joiner(self, k, indent_level):
+        return (',\n' + ' '*indent_level) if k in self.indented_args else ', '
+
+    def _indented(self, printer, k, v, *args, **kwargs):
+        il = printer._context['indent_level']
+        def _print(arg):
+            if isinstance(arg, Token):
+                return printer._print(arg, *args, joiner=self._joiner(k, il), **kwargs)
+            else:
+                return printer._print(arg, *args, **kwargs)
+
+        if isinstance(v, Tuple):
+            joined = self._joiner(k, il).join([_print(arg) for arg in v.args])
+            if k in self.indented_args:
+                return '(\n' + ' '*il + joined + ',\n' + ' '*(il - 4) + ')'
+            else:
+                return ('({0},)' if len(v.args) == 1 else '({0})').format(joined)
+        else:
+            return _print(v)
+
+    def _sympyrepr(self, printer, *args, joiner=', ', **kwargs):
+        from sympy.printing.printer import printer_context
+        exclude = kwargs.get('exclude', ())
+        values = [getattr(self, k) for k in self._fields]
+        indent_level = printer._context.get('indent_level', 0)
+
+        arg_reprs = []
+
+        for i, (attr, value) in enumerate(zip(self._fields, values)):
+            if attr in exclude:
+                continue
+
+            # Skip attributes which have the default value
+            if attr in self.defaults and value == self.defaults[attr]:
+                continue
+
+            ilvl = indent_level + 4 if attr in self.indented_args else 0
+            with printer_context(printer, indent_level=ilvl):
+                indented = self._indented(printer, attr, value, *args, **kwargs)
+            arg_reprs.append(('{1}' if i == 0 else '{0}={1}').format(attr, indented.lstrip()))
+
+        return "{}({})".format(self.__class__.__name__, joiner.join(arg_reprs))
+
+    _sympystr = _sympyrepr
+
+    def __repr__(self):  # sympy.core.Basic.__repr__ uses sstr
+        from sympy.printing import srepr
+        return srepr(self)
+
+    def kwargs(self, exclude=(), apply=None):
+        """ Get instance's attributes as dict of keyword arguments.
+
+        Parameters
+        ==========
+
+        exclude : collection of str
+            Collection of keywords to exclude.
+
+        apply : callable, optional
+            Function to apply to all values.
+        """
+        kwargs = {k: getattr(self, k) for k in self._fields if k not in exclude}
+        if apply is not None:
+            return {k: apply(v) for k, v in kwargs.items()}
+        else:
+            return kwargs
+
+class BreakToken(Token):
+    """ Represents 'break' in C/Python ('exit' in Fortran).
+
+    Use the premade instance ``break_`` or instantiate manually.
+
+    Examples
+    ========
+
+    >>> from sympy import ccode, fcode
+    >>> from sympy.codegen.ast import break_
+    >>> ccode(break_)
+    'break'
+    >>> fcode(break_, source_format='free')
+    'exit'
+    """
+
+break_ = BreakToken()
+
+
+class ContinueToken(Token):
+    """ Represents 'continue' in C/Python ('cycle' in Fortran)
+
+    Use the premade instance ``continue_`` or instantiate manually.
+
+    Examples
+    ========
+
+    >>> from sympy import ccode, fcode
+    >>> from sympy.codegen.ast import continue_
+    >>> ccode(continue_)
+    'continue'
+    >>> fcode(continue_, source_format='free')
+    'cycle'
+    """
+
+continue_ = ContinueToken()
+
+class NoneToken(Token):
+    """ The AST equivalence of Python's NoneType
+
+    The corresponding instance of Python's ``None`` is ``none``.
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.ast import none, Variable
+    >>> from sympy import pycode
+    >>> print(pycode(Variable('x').as_Declaration(value=none)))
+    x = None
+
+    """
+    def __eq__(self, other):
+        return other is None or isinstance(other, NoneToken)
+
+    def _hashable_content(self):
+        return ()
+
+    def __hash__(self):
+        return super().__hash__()
+
+
+none = NoneToken()
+
+
+class AssignmentBase(CodegenAST):
+    """ Abstract base class for Assignment and AugmentedAssignment.
+
+    Attributes:
+    ===========
+
+    op : str
+        Symbol for assignment operator, e.g. "=", "+=", etc.
+    """
+
+    def __new__(cls, lhs, rhs):
+        lhs = _sympify(lhs)
+        rhs = _sympify(rhs)
+
+        cls._check_args(lhs, rhs)
+
+        return super().__new__(cls, lhs, rhs)
+
+    @property
+    def lhs(self):
+        return self.args[0]
+
+    @property
+    def rhs(self):
+        return self.args[1]
+
+    @classmethod
+    def _check_args(cls, lhs, rhs):
+        """ Check arguments to __new__ and raise exception if any problems found.
+
+        Derived classes may wish to override this.
+        """
+        from sympy.matrices.expressions.matexpr import (
+            MatrixElement, MatrixSymbol)
+        from sympy.tensor.indexed import Indexed
+        from sympy.tensor.array.expressions import ArrayElement
+
+        # Tuple of things that can be on the lhs of an assignment
+        assignable = (Symbol, MatrixSymbol, MatrixElement, Indexed, Element, Variable,
+                ArrayElement)
+        if not isinstance(lhs, assignable):
+            raise TypeError("Cannot assign to lhs of type %s." % type(lhs))
+
+        # Indexed types implement shape, but don't define it until later. This
+        # causes issues in assignment validation. For now, matrices are defined
+        # as anything with a shape that is not an Indexed
+        lhs_is_mat = hasattr(lhs, 'shape') and not isinstance(lhs, Indexed)
+        rhs_is_mat = hasattr(rhs, 'shape') and not isinstance(rhs, Indexed)
+
+        # If lhs and rhs have same structure, then this assignment is ok
+        if lhs_is_mat:
+            if not rhs_is_mat:
+                raise ValueError("Cannot assign a scalar to a matrix.")
+            elif lhs.shape != rhs.shape:
+                raise ValueError("Dimensions of lhs and rhs do not align.")
+        elif rhs_is_mat and not lhs_is_mat:
+            raise ValueError("Cannot assign a matrix to a scalar.")
+
+
+class Assignment(AssignmentBase):
+    """
+    Represents variable assignment for code generation.
+
+    Parameters
+    ==========
+
+    lhs : Expr
+        SymPy object representing the lhs of the expression. These should be
+        singular objects, such as one would use in writing code. Notable types
+        include Symbol, MatrixSymbol, MatrixElement, and Indexed. Types that
+        subclass these types are also supported.
+
+    rhs : Expr
+        SymPy object representing the rhs of the expression. This can be any
+        type, provided its shape corresponds to that of the lhs. For example,
+        a Matrix type can be assigned to MatrixSymbol, but not to Symbol, as
+        the dimensions will not align.
+
+    Examples
+    ========
+
+    >>> from sympy import symbols, MatrixSymbol, Matrix
+    >>> from sympy.codegen.ast import Assignment
+    >>> x, y, z = symbols('x, y, z')
+    >>> Assignment(x, y)
+    Assignment(x, y)
+    >>> Assignment(x, 0)
+    Assignment(x, 0)
+    >>> A = MatrixSymbol('A', 1, 3)
+    >>> mat = Matrix([x, y, z]).T
+    >>> Assignment(A, mat)
+    Assignment(A, Matrix([[x, y, z]]))
+    >>> Assignment(A[0, 1], x)
+    Assignment(A[0, 1], x)
+    """
+
+    op = ':='
+
+
+class AugmentedAssignment(AssignmentBase):
+    """
+    Base class for augmented assignments.
+
+    Attributes:
+    ===========
+
+    binop : str
+       Symbol for binary operation being applied in the assignment, such as "+",
+       "*", etc.
+    """
+    binop = None  # type: str
+
+    @property
+    def op(self):
+        return self.binop + '='
+
+
+class AddAugmentedAssignment(AugmentedAssignment):
+    binop = '+'
+
+
+class SubAugmentedAssignment(AugmentedAssignment):
+    binop = '-'
+
+
+class MulAugmentedAssignment(AugmentedAssignment):
+    binop = '*'
+
+
+class DivAugmentedAssignment(AugmentedAssignment):
+    binop = '/'
+
+
+class ModAugmentedAssignment(AugmentedAssignment):
+    binop = '%'
+
+
+# Mapping from binary op strings to AugmentedAssignment subclasses
+augassign_classes = {
+    cls.binop: cls for cls in [
+        AddAugmentedAssignment, SubAugmentedAssignment, MulAugmentedAssignment,
+        DivAugmentedAssignment, ModAugmentedAssignment
+    ]
+}
+
+
+def aug_assign(lhs, op, rhs):
+    """
+    Create 'lhs op= rhs'.
+
+    Explanation
+    ===========
+
+    Represents augmented variable assignment for code generation. This is a
+    convenience function. You can also use the AugmentedAssignment classes
+    directly, like AddAugmentedAssignment(x, y).
+
+    Parameters
+    ==========
+
+    lhs : Expr
+        SymPy object representing the lhs of the expression. These should be
+        singular objects, such as one would use in writing code. Notable types
+        include Symbol, MatrixSymbol, MatrixElement, and Indexed. Types that
+        subclass these types are also supported.
+
+    op : str
+        Operator (+, -, /, \\*, %).
+
+    rhs : Expr
+        SymPy object representing the rhs of the expression. This can be any
+        type, provided its shape corresponds to that of the lhs. For example,
+        a Matrix type can be assigned to MatrixSymbol, but not to Symbol, as
+        the dimensions will not align.
+
+    Examples
+    ========
+
+    >>> from sympy import symbols
+    >>> from sympy.codegen.ast import aug_assign
+    >>> x, y = symbols('x, y')
+    >>> aug_assign(x, '+', y)
+    AddAugmentedAssignment(x, y)
+    """
+    if op not in augassign_classes:
+        raise ValueError("Unrecognized operator %s" % op)
+    return augassign_classes[op](lhs, rhs)
+
+
+class CodeBlock(CodegenAST):
+    """
+    Represents a block of code.
+
+    Explanation
+    ===========
+
+    For now only assignments are supported. This restriction will be lifted in
+    the future.
+
+    Useful attributes on this object are:
+
+    ``left_hand_sides``:
+        Tuple of left-hand sides of assignments, in order.
+    ``left_hand_sides``:
+        Tuple of right-hand sides of assignments, in order.
+    ``free_symbols``: Free symbols of the expressions in the right-hand sides
+        which do not appear in the left-hand side of an assignment.
+
+    Useful methods on this object are:
+
+    ``topological_sort``:
+        Class method. Return a CodeBlock with assignments
+        sorted so that variables are assigned before they
+        are used.
+    ``cse``:
+        Return a new CodeBlock with common subexpressions eliminated and
+        pulled out as assignments.
+
+    Examples
+    ========
+
+    >>> from sympy import symbols, ccode
+    >>> from sympy.codegen.ast import CodeBlock, Assignment
+    >>> x, y = symbols('x y')
+    >>> c = CodeBlock(Assignment(x, 1), Assignment(y, x + 1))
+    >>> print(ccode(c))
+    x = 1;
+    y = x + 1;
+
+    """
+    def __new__(cls, *args):
+        left_hand_sides = []
+        right_hand_sides = []
+        for i in args:
+            if isinstance(i, Assignment):
+                lhs, rhs = i.args
+                left_hand_sides.append(lhs)
+                right_hand_sides.append(rhs)
+
+        obj = CodegenAST.__new__(cls, *args)
+
+        obj.left_hand_sides = Tuple(*left_hand_sides)
+        obj.right_hand_sides = Tuple(*right_hand_sides)
+        return obj
+
+    def __iter__(self):
+        return iter(self.args)
+
+    def _sympyrepr(self, printer, *args, **kwargs):
+        il = printer._context.get('indent_level', 0)
+        joiner = ',\n' + ' '*il
+        joined = joiner.join(map(printer._print, self.args))
+        return ('{}(\n'.format(' '*(il-4) + self.__class__.__name__,) +
+                ' '*il + joined + '\n' + ' '*(il - 4) + ')')
+
+    _sympystr = _sympyrepr
+
+    @property
+    def free_symbols(self):
+        return super().free_symbols - set(self.left_hand_sides)
+
+    @classmethod
+    def topological_sort(cls, assignments):
+        """
+        Return a CodeBlock with topologically sorted assignments so that
+        variables are assigned before they are used.
+
+        Examples
+        ========
+
+        The existing order of assignments is preserved as much as possible.
+
+        This function assumes that variables are assigned to only once.
+
+        This is a class constructor so that the default constructor for
+        CodeBlock can error when variables are used before they are assigned.
+
+        >>> from sympy import symbols
+        >>> from sympy.codegen.ast import CodeBlock, Assignment
+        >>> x, y, z = symbols('x y z')
+
+        >>> assignments = [
+        ...     Assignment(x, y + z),
+        ...     Assignment(y, z + 1),
+        ...     Assignment(z, 2),
+        ... ]
+        >>> CodeBlock.topological_sort(assignments)
+        CodeBlock(
+            Assignment(z, 2),
+            Assignment(y, z + 1),
+            Assignment(x, y + z)
+        )
+
+        """
+
+        if not all(isinstance(i, Assignment) for i in assignments):
+            # Will support more things later
+            raise NotImplementedError("CodeBlock.topological_sort only supports Assignments")
+
+        if any(isinstance(i, AugmentedAssignment) for i in assignments):
+            raise NotImplementedError("CodeBlock.topological_sort does not yet work with AugmentedAssignments")
+
+        # Create a graph where the nodes are assignments and there is a directed edge
+        # between nodes that use a variable and nodes that assign that
+        # variable, like
+
+        # [(x := 1, y := x + 1), (x := 1, z := y + z), (y := x + 1, z := y + z)]
+
+        # If we then topologically sort these nodes, they will be in
+        # assignment order, like
+
+        # x := 1
+        # y := x + 1
+        # z := y + z
+
+        # A = The nodes
+        #
+        # enumerate keeps nodes in the same order they are already in if
+        # possible. It will also allow us to handle duplicate assignments to
+        # the same variable when those are implemented.
+        A = list(enumerate(assignments))
+
+        # var_map = {variable: [nodes for which this variable is assigned to]}
+        # like {x: [(1, x := y + z), (4, x := 2 * w)], ...}
+        var_map = defaultdict(list)
+        for node in A:
+            i, a = node
+            var_map[a.lhs].append(node)
+
+        # E = Edges in the graph
+        E = []
+        for dst_node in A:
+            i, a = dst_node
+            for s in a.rhs.free_symbols:
+                for src_node in var_map[s]:
+                    E.append((src_node, dst_node))
+
+        ordered_assignments = topological_sort([A, E])
+
+        # De-enumerate the result
+        return cls(*[a for i, a in ordered_assignments])
+
+    def cse(self, symbols=None, optimizations=None, postprocess=None,
+        order='canonical'):
+        """
+        Return a new code block with common subexpressions eliminated.
+
+        Explanation
+        ===========
+
+        See the docstring of :func:`sympy.simplify.cse_main.cse` for more
+        information.
+
+        Examples
+        ========
+
+        >>> from sympy import symbols, sin
+        >>> from sympy.codegen.ast import CodeBlock, Assignment
+        >>> x, y, z = symbols('x y z')
+
+        >>> c = CodeBlock(
+        ...     Assignment(x, 1),
+        ...     Assignment(y, sin(x) + 1),
+        ...     Assignment(z, sin(x) - 1),
+        ... )
+        ...
+        >>> c.cse()
+        CodeBlock(
+            Assignment(x, 1),
+            Assignment(x0, sin(x)),
+            Assignment(y, x0 + 1),
+            Assignment(z, x0 - 1)
+        )
+
+        """
+        from sympy.simplify.cse_main import cse
+
+        # Check that the CodeBlock only contains assignments to unique variables
+        if not all(isinstance(i, Assignment) for i in self.args):
+            # Will support more things later
+            raise NotImplementedError("CodeBlock.cse only supports Assignments")
+
+        if any(isinstance(i, AugmentedAssignment) for i in self.args):
+            raise NotImplementedError("CodeBlock.cse does not yet work with AugmentedAssignments")
+
+        for i, lhs in enumerate(self.left_hand_sides):
+            if lhs in self.left_hand_sides[:i]:
+                raise NotImplementedError("Duplicate assignments to the same "
+                    "variable are not yet supported (%s)" % lhs)
+
+        # Ensure new symbols for subexpressions do not conflict with existing
+        existing_symbols = self.atoms(Symbol)
+        if symbols is None:
+            symbols = numbered_symbols()
+        symbols = filter_symbols(symbols, existing_symbols)
+
+        replacements, reduced_exprs = cse(list(self.right_hand_sides),
+            symbols=symbols, optimizations=optimizations, postprocess=postprocess,
+            order=order)
+
+        new_block = [Assignment(var, expr) for var, expr in
+            zip(self.left_hand_sides, reduced_exprs)]
+        new_assignments = [Assignment(var, expr) for var, expr in replacements]
+        return self.topological_sort(new_assignments + new_block)
+
+
+class For(Token):
+    """Represents a 'for-loop' in the code.
+
+    Expressions are of the form:
+        "for target in iter:
+            body..."
+
+    Parameters
+    ==========
+
+    target : symbol
+    iter : iterable
+    body : CodeBlock or iterable
+!        When passed an iterable it is used to instantiate a CodeBlock.
+
+    Examples
+    ========
+
+    >>> from sympy import symbols, Range
+    >>> from sympy.codegen.ast import aug_assign, For
+    >>> x, i, j, k = symbols('x i j k')
+    >>> for_i = For(i, Range(10), [aug_assign(x, '+', i*j*k)])
+    >>> for_i  # doctest: -NORMALIZE_WHITESPACE
+    For(i, iterable=Range(0, 10, 1), body=CodeBlock(
+        AddAugmentedAssignment(x, i*j*k)
+    ))
+    >>> for_ji = For(j, Range(7), [for_i])
+    >>> for_ji  # doctest: -NORMALIZE_WHITESPACE
+    For(j, iterable=Range(0, 7, 1), body=CodeBlock(
+        For(i, iterable=Range(0, 10, 1), body=CodeBlock(
+            AddAugmentedAssignment(x, i*j*k)
+        ))
+    ))
+    >>> for_kji =For(k, Range(5), [for_ji])
+    >>> for_kji  # doctest: -NORMALIZE_WHITESPACE
+    For(k, iterable=Range(0, 5, 1), body=CodeBlock(
+        For(j, iterable=Range(0, 7, 1), body=CodeBlock(
+            For(i, iterable=Range(0, 10, 1), body=CodeBlock(
+                AddAugmentedAssignment(x, i*j*k)
+            ))
+        ))
+    ))
+    """
+    __slots__ = _fields = ('target', 'iterable', 'body')
+    _construct_target = staticmethod(_sympify)
+
+    @classmethod
+    def _construct_body(cls, itr):
+        if isinstance(itr, CodeBlock):
+            return itr
+        else:
+            return CodeBlock(*itr)
+
+    @classmethod
+    def _construct_iterable(cls, itr):
+        if not iterable(itr):
+            raise TypeError("iterable must be an iterable")
+        if isinstance(itr, list):  # _sympify errors on lists because they are mutable
+            itr = tuple(itr)
+        return _sympify(itr)
+
+
+class String(Atom, Token):
+    """ SymPy object representing a string.
+
+    Atomic object which is not an expression (as opposed to Symbol).
+
+    Parameters
+    ==========
+
+    text : str
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.ast import String
+    >>> f = String('foo')
+    >>> f
+    foo
+    >>> str(f)
+    'foo'
+    >>> f.text
+    'foo'
+    >>> print(repr(f))
+    String('foo')
+
+    """
+    __slots__ = _fields = ('text',)
+    not_in_args = ['text']
+    is_Atom = True
+
+    @classmethod
+    def _construct_text(cls, text):
+        if not isinstance(text, str):
+            raise TypeError("Argument text is not a string type.")
+        return text
+
+    def _sympystr(self, printer, *args, **kwargs):
+        return self.text
+
+    def kwargs(self, exclude = (), apply = None):
+        return {}
+
+    #to be removed when Atom is given a suitable func
+    @property
+    def func(self):
+        return lambda: self
+
+    def _latex(self, printer):
+        from sympy.printing.latex import latex_escape
+        return r'\texttt{{"{}"}}'.format(latex_escape(self.text))
+
+class QuotedString(String):
+    """ Represents a string which should be printed with quotes. """
+
+class Comment(String):
+    """ Represents a comment. """
+
+class Node(Token):
+    """ Subclass of Token, carrying the attribute 'attrs' (Tuple)
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.ast import Node, value_const, pointer_const
+    >>> n1 = Node([value_const])
+    >>> n1.attr_params('value_const')  # get the parameters of attribute (by name)
+    ()
+    >>> from sympy.codegen.fnodes import dimension
+    >>> n2 = Node([value_const, dimension(5, 3)])
+    >>> n2.attr_params(value_const)  # get the parameters of attribute (by Attribute instance)
+    ()
+    >>> n2.attr_params('dimension')  # get the parameters of attribute (by name)
+    (5, 3)
+    >>> n2.attr_params(pointer_const) is None
+    True
+
+    """
+
+    __slots__: tuple[str, ...] = ('attrs',)
+    _fields = __slots__
+
+    defaults: dict[str, Any] = {'attrs': Tuple()}
+
+    _construct_attrs = staticmethod(_mk_Tuple)
+
+    def attr_params(self, looking_for):
+        """ Returns the parameters of the Attribute with name ``looking_for`` in self.attrs """
+        for attr in self.attrs:
+            if str(attr.name) == str(looking_for):
+                return attr.parameters
+
+
+class Type(Token):
+    """ Represents a type.
+
+    Explanation
+    ===========
+
+    The naming is a super-set of NumPy naming. Type has a classmethod
+    ``from_expr`` which offer type deduction. It also has a method
+    ``cast_check`` which casts the argument to its type, possibly raising an
+    exception if rounding error is not within tolerances, or if the value is not
+    representable by the underlying data type (e.g. unsigned integers).
+
+    Parameters
+    ==========
+
+    name : str
+        Name of the type, e.g. ``object``, ``int16``, ``float16`` (where the latter two
+        would use the ``Type`` sub-classes ``IntType`` and ``FloatType`` respectively).
+        If a ``Type`` instance is given, the said instance is returned.
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.ast import Type
+    >>> t = Type.from_expr(42)
+    >>> t
+    integer
+    >>> print(repr(t))
+    IntBaseType(String('integer'))
+    >>> from sympy.codegen.ast import uint8
+    >>> uint8.cast_check(-1)   # doctest: +ELLIPSIS
+    Traceback (most recent call last):
+      ...
+    ValueError: Minimum value for data type bigger than new value.
+    >>> from sympy.codegen.ast import float32
+    >>> v6 = 0.123456
+    >>> float32.cast_check(v6)
+    0.123456
+    >>> v10 = 12345.67894
+    >>> float32.cast_check(v10)  # doctest: +ELLIPSIS
+    Traceback (most recent call last):
+      ...
+    ValueError: Casting gives a significantly different value.
+    >>> boost_mp50 = Type('boost::multiprecision::cpp_dec_float_50')
+    >>> from sympy import cxxcode
+    >>> from sympy.codegen.ast import Declaration, Variable
+    >>> cxxcode(Declaration(Variable('x', type=boost_mp50)))
+    'boost::multiprecision::cpp_dec_float_50 x'
+
+    References
+    ==========
+
+    .. [1] https://numpy.org/doc/stable/user/basics.types.html
+
+    """
+    __slots__: tuple[str, ...] = ('name',)
+    _fields = __slots__
+
+    _construct_name = String
+
+    def _sympystr(self, printer, *args, **kwargs):
+        return str(self.name)
+
+    @classmethod
+    def from_expr(cls, expr):
+        """ Deduces type from an expression or a ``Symbol``.
+
+        Parameters
+        ==========
+
+        expr : number or SymPy object
+            The type will be deduced from type or properties.
+
+        Examples
+        ========
+
+        >>> from sympy.codegen.ast import Type, integer, complex_
+        >>> Type.from_expr(2) == integer
+        True
+        >>> from sympy import Symbol
+        >>> Type.from_expr(Symbol('z', complex=True)) == complex_
+        True
+        >>> Type.from_expr(sum)  # doctest: +ELLIPSIS
+        Traceback (most recent call last):
+          ...
+        ValueError: Could not deduce type from expr.
+
+        Raises
+        ======
+
+        ValueError when type deduction fails.
+
+        """
+        if isinstance(expr, (float, Float)):
+            return real
+        if isinstance(expr, (int, Integer)) or getattr(expr, 'is_integer', False):
+            return integer
+        if getattr(expr, 'is_real', False):
+            return real
+        if isinstance(expr, complex) or getattr(expr, 'is_complex', False):
+            return complex_
+        if isinstance(expr, bool) or getattr(expr, 'is_Relational', False):
+            return bool_
+        else:
+            raise ValueError("Could not deduce type from expr.")
+
+    def _check(self, value):
+        pass
+
+    def cast_check(self, value, rtol=None, atol=0, precision_targets=None):
+        """ Casts a value to the data type of the instance.
+
+        Parameters
+        ==========
+
+        value : number
+        rtol : floating point number
+            Relative tolerance. (will be deduced if not given).
+        atol : floating point number
+            Absolute tolerance (in addition to ``rtol``).
+        type_aliases : dict
+            Maps substitutions for Type, e.g. {integer: int64, real: float32}
+
+        Examples
+        ========
+
+        >>> from sympy.codegen.ast import integer, float32, int8
+        >>> integer.cast_check(3.0) == 3
+        True
+        >>> float32.cast_check(1e-40)  # doctest: +ELLIPSIS
+        Traceback (most recent call last):
+          ...
+        ValueError: Minimum value for data type bigger than new value.
+        >>> int8.cast_check(256)  # doctest: +ELLIPSIS
+        Traceback (most recent call last):
+          ...
+        ValueError: Maximum value for data type smaller than new value.
+        >>> v10 = 12345.67894
+        >>> float32.cast_check(v10)  # doctest: +ELLIPSIS
+        Traceback (most recent call last):
+          ...
+        ValueError: Casting gives a significantly different value.
+        >>> from sympy.codegen.ast import float64
+        >>> float64.cast_check(v10)
+        12345.67894
+        >>> from sympy import Float
+        >>> v18 = Float('0.123456789012345646')
+        >>> float64.cast_check(v18)
+        Traceback (most recent call last):
+          ...
+        ValueError: Casting gives a significantly different value.
+        >>> from sympy.codegen.ast import float80
+        >>> float80.cast_check(v18)
+        0.123456789012345649
+
+        """
+        val = sympify(value)
+
+        ten = Integer(10)
+        exp10 = getattr(self, 'decimal_dig', None)
+
+        if rtol is None:
+            rtol = 1e-15 if exp10 is None else 2.0*ten**(-exp10)
+
+        def tol(num):
+            return atol + rtol*abs(num)
+
+        new_val = self.cast_nocheck(value)
+        self._check(new_val)
+
+        delta = new_val - val
+        if abs(delta) > tol(val):  # rounding, e.g. int(3.5) != 3.5
+            raise ValueError("Casting gives a significantly different value.")
+
+        return new_val
+
+    def _latex(self, printer):
+        from sympy.printing.latex import latex_escape
+        type_name = latex_escape(self.__class__.__name__)
+        name = latex_escape(self.name.text)
+        return r"\text{{{}}}\left(\texttt{{{}}}\right)".format(type_name, name)
+
+
+class IntBaseType(Type):
+    """ Integer base type, contains no size information. """
+    __slots__ = ()
+    cast_nocheck = lambda self, i: Integer(int(i))
+
+
+class _SizedIntType(IntBaseType):
+    __slots__ = ('nbits',)
+    _fields = Type._fields + __slots__
+
+    _construct_nbits = Integer
+
+    def _check(self, value):
+        if value < self.min:
+            raise ValueError("Value is too small: %d < %d" % (value, self.min))
+        if value > self.max:
+            raise ValueError("Value is too big: %d > %d" % (value, self.max))
+
+
+class SignedIntType(_SizedIntType):
+    """ Represents a signed integer type. """
+    __slots__ = ()
+    @property
+    def min(self):
+        return -2**(self.nbits-1)
+
+    @property
+    def max(self):
+        return 2**(self.nbits-1) - 1
+
+
+class UnsignedIntType(_SizedIntType):
+    """ Represents an unsigned integer type. """
+    __slots__ = ()
+    @property
+    def min(self):
+        return 0
+
+    @property
+    def max(self):
+        return 2**self.nbits - 1
+
+two = Integer(2)
+
+class FloatBaseType(Type):
+    """ Represents a floating point number type. """
+    __slots__ = ()
+    cast_nocheck = Float
+
+class FloatType(FloatBaseType):
+    """ Represents a floating point type with fixed bit width.
+
+    Base 2 & one sign bit is assumed.
+
+    Parameters
+    ==========
+
+    name : str
+        Name of the type.
+    nbits : integer
+        Number of bits used (storage).
+    nmant : integer
+        Number of bits used to represent the mantissa.
+    nexp : integer
+        Number of bits used to represent the mantissa.
+
+    Examples
+    ========
+
+    >>> from sympy import S
+    >>> from sympy.codegen.ast import FloatType
+    >>> half_precision = FloatType('f16', nbits=16, nmant=10, nexp=5)
+    >>> half_precision.max
+    65504
+    >>> half_precision.tiny == S(2)**-14
+    True
+    >>> half_precision.eps == S(2)**-10
+    True
+    >>> half_precision.dig == 3
+    True
+    >>> half_precision.decimal_dig == 5
+    True
+    >>> half_precision.cast_check(1.0)
+    1.0
+    >>> half_precision.cast_check(1e5)  # doctest: +ELLIPSIS
+    Traceback (most recent call last):
+      ...
+    ValueError: Maximum value for data type smaller than new value.
+    """
+
+    __slots__ = ('nbits', 'nmant', 'nexp',)
+    _fields = Type._fields + __slots__
+
+    _construct_nbits = _construct_nmant = _construct_nexp = Integer
+
+
+    @property
+    def max_exponent(self):
+        """ The largest positive number n, such that 2**(n - 1) is a representable finite value. """
+        # cf. C++'s ``std::numeric_limits::max_exponent``
+        return two**(self.nexp - 1)
+
+    @property
+    def min_exponent(self):
+        """ The lowest negative number n, such that 2**(n - 1) is a valid normalized number. """
+        # cf. C++'s ``std::numeric_limits::min_exponent``
+        return 3 - self.max_exponent
+
+    @property
+    def max(self):
+        """ Maximum value representable. """
+        return (1 - two**-(self.nmant+1))*two**self.max_exponent
+
+    @property
+    def tiny(self):
+        """ The minimum positive normalized value. """
+        # See C macros: FLT_MIN, DBL_MIN, LDBL_MIN
+        # or C++'s ``std::numeric_limits::min``
+        # or numpy.finfo(dtype).tiny
+        return two**(self.min_exponent - 1)
+
+
+    @property
+    def eps(self):
+        """ Difference between 1.0 and the next representable value. """
+        return two**(-self.nmant)
+
+    @property
+    def dig(self):
+        """ Number of decimal digits that are guaranteed to be preserved in text.
+
+        When converting text -> float -> text, you are guaranteed that at least ``dig``
+        number of digits are preserved with respect to rounding or overflow.
+        """
+        from sympy.functions import floor, log
+        return floor(self.nmant * log(2)/log(10))
+
+    @property
+    def decimal_dig(self):
+        """ Number of digits needed to store & load without loss.
+
+        Explanation
+        ===========
+
+        Number of decimal digits needed to guarantee that two consecutive conversions
+        (float -> text -> float) to be idempotent. This is useful when one do not want
+        to loose precision due to rounding errors when storing a floating point value
+        as text.
+        """
+        from sympy.functions import ceiling, log
+        return ceiling((self.nmant + 1) * log(2)/log(10) + 1)
+
+    def cast_nocheck(self, value):
+        """ Casts without checking if out of bounds or subnormal. """
+        if value == oo:  # float(oo) or oo
+            return float(oo)
+        elif value == -oo:  # float(-oo) or -oo
+            return float(-oo)
+        return Float(str(sympify(value).evalf(self.decimal_dig)), self.decimal_dig)
+
+    def _check(self, value):
+        if value < -self.max:
+            raise ValueError("Value is too small: %d < %d" % (value, -self.max))
+        if value > self.max:
+            raise ValueError("Value is too big: %d > %d" % (value, self.max))
+        if abs(value) < self.tiny:
+            raise ValueError("Smallest (absolute) value for data type bigger than new value.")
+
+class ComplexBaseType(FloatBaseType):
+
+    __slots__ = ()
+
+    def cast_nocheck(self, value):
+        """ Casts without checking if out of bounds or subnormal. """
+        from sympy.functions import re, im
+        return (
+            super().cast_nocheck(re(value)) +
+            super().cast_nocheck(im(value))*1j
+        )
+
+    def _check(self, value):
+        from sympy.functions import re, im
+        super()._check(re(value))
+        super()._check(im(value))
+
+
+class ComplexType(ComplexBaseType, FloatType):
+    """ Represents a complex floating point number. """
+    __slots__ = ()
+
+
+# NumPy types:
+intc = IntBaseType('intc')
+intp = IntBaseType('intp')
+int8 = SignedIntType('int8', 8)
+int16 = SignedIntType('int16', 16)
+int32 = SignedIntType('int32', 32)
+int64 = SignedIntType('int64', 64)
+uint8 = UnsignedIntType('uint8', 8)
+uint16 = UnsignedIntType('uint16', 16)
+uint32 = UnsignedIntType('uint32', 32)
+uint64 = UnsignedIntType('uint64', 64)
+float16 = FloatType('float16', 16, nexp=5, nmant=10)  # IEEE 754 binary16, Half precision
+float32 = FloatType('float32', 32, nexp=8, nmant=23)  # IEEE 754 binary32, Single precision
+float64 = FloatType('float64', 64, nexp=11, nmant=52)  # IEEE 754 binary64, Double precision
+float80 = FloatType('float80', 80, nexp=15, nmant=63)  # x86 extended precision (1 integer part bit), "long double"
+float128 = FloatType('float128', 128, nexp=15, nmant=112)  # IEEE 754 binary128, Quadruple precision
+float256 = FloatType('float256', 256, nexp=19, nmant=236)  # IEEE 754 binary256, Octuple precision
+
+complex64 = ComplexType('complex64', nbits=64, **float32.kwargs(exclude=('name', 'nbits')))
+complex128 = ComplexType('complex128', nbits=128, **float64.kwargs(exclude=('name', 'nbits')))
+
+# Generic types (precision may be chosen by code printers):
+untyped = Type('untyped')
+real = FloatBaseType('real')
+integer = IntBaseType('integer')
+complex_ = ComplexBaseType('complex')
+bool_ = Type('bool')
+
+
+class Attribute(Token):
+    """ Attribute (possibly parametrized)
+
+    For use with :class:`sympy.codegen.ast.Node` (which takes instances of
+    ``Attribute`` as ``attrs``).
+
+    Parameters
+    ==========
+
+    name : str
+    parameters : Tuple
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.ast import Attribute
+    >>> volatile = Attribute('volatile')
+    >>> volatile
+    volatile
+    >>> print(repr(volatile))
+    Attribute(String('volatile'))
+    >>> a = Attribute('foo', [1, 2, 3])
+    >>> a
+    foo(1, 2, 3)
+    >>> a.parameters == (1, 2, 3)
+    True
+    """
+    __slots__ = _fields = ('name', 'parameters')
+    defaults = {'parameters': Tuple()}
+
+    _construct_name = String
+    _construct_parameters = staticmethod(_mk_Tuple)
+
+    def _sympystr(self, printer, *args, **kwargs):
+        result = str(self.name)
+        if self.parameters:
+            result += '(%s)' % ', '.join((printer._print(
+                arg, *args, **kwargs) for arg in self.parameters))
+        return result
+
+value_const = Attribute('value_const')
+pointer_const = Attribute('pointer_const')
+
+
+class Variable(Node):
+    """ Represents a variable.
+
+    Parameters
+    ==========
+
+    symbol : Symbol
+    type : Type (optional)
+        Type of the variable.
+    attrs : iterable of Attribute instances
+        Will be stored as a Tuple.
+
+    Examples
+    ========
+
+    >>> from sympy import Symbol
+    >>> from sympy.codegen.ast import Variable, float32, integer
+    >>> x = Symbol('x')
+    >>> v = Variable(x, type=float32)
+    >>> v.attrs
+    ()
+    >>> v == Variable('x')
+    False
+    >>> v == Variable('x', type=float32)
+    True
+    >>> v
+    Variable(x, type=float32)
+
+    One may also construct a ``Variable`` instance with the type deduced from
+    assumptions about the symbol using the ``deduced`` classmethod:
+
+    >>> i = Symbol('i', integer=True)
+    >>> v = Variable.deduced(i)
+    >>> v.type == integer
+    True
+    >>> v == Variable('i')
+    False
+    >>> from sympy.codegen.ast import value_const
+    >>> value_const in v.attrs
+    False
+    >>> w = Variable('w', attrs=[value_const])
+    >>> w
+    Variable(w, attrs=(value_const,))
+    >>> value_const in w.attrs
+    True
+    >>> w.as_Declaration(value=42)
+    Declaration(Variable(w, value=42, attrs=(value_const,)))
+
+    """
+
+    __slots__ = ('symbol', 'type', 'value')
+    _fields = __slots__ + Node._fields
+
+    defaults = Node.defaults.copy()
+    defaults.update({'type': untyped, 'value': none})
+
+    _construct_symbol = staticmethod(sympify)
+    _construct_value = staticmethod(sympify)
+
+    @classmethod
+    def deduced(cls, symbol, value=None, attrs=Tuple(), cast_check=True):
+        """ Alt. constructor with type deduction from ``Type.from_expr``.
+
+        Deduces type primarily from ``symbol``, secondarily from ``value``.
+
+        Parameters
+        ==========
+
+        symbol : Symbol
+        value : expr
+            (optional) value of the variable.
+        attrs : iterable of Attribute instances
+        cast_check : bool
+            Whether to apply ``Type.cast_check`` on ``value``.
+
+        Examples
+        ========
+
+        >>> from sympy import Symbol
+        >>> from sympy.codegen.ast import Variable, complex_
+        >>> n = Symbol('n', integer=True)
+        >>> str(Variable.deduced(n).type)
+        'integer'
+        >>> x = Symbol('x', real=True)
+        >>> v = Variable.deduced(x)
+        >>> v.type
+        real
+        >>> z = Symbol('z', complex=True)
+        >>> Variable.deduced(z).type == complex_
+        True
+
+        """
+        if isinstance(symbol, Variable):
+            return symbol
+
+        try:
+            type_ = Type.from_expr(symbol)
+        except ValueError:
+            type_ = Type.from_expr(value)
+
+        if value is not None and cast_check:
+            value = type_.cast_check(value)
+        return cls(symbol, type=type_, value=value, attrs=attrs)
+
+    def as_Declaration(self, **kwargs):
+        """ Convenience method for creating a Declaration instance.
+
+        Explanation
+        ===========
+
+        If the variable of the Declaration need to wrap a modified
+        variable keyword arguments may be passed (overriding e.g.
+        the ``value`` of the Variable instance).
+
+        Examples
+        ========
+
+        >>> from sympy.codegen.ast import Variable, NoneToken
+        >>> x = Variable('x')
+        >>> decl1 = x.as_Declaration()
+        >>> # value is special NoneToken() which must be tested with == operator
+        >>> decl1.variable.value is None  # won't work
+        False
+        >>> decl1.variable.value == None  # not PEP-8 compliant
+        True
+        >>> decl1.variable.value == NoneToken()  # OK
+        True
+        >>> decl2 = x.as_Declaration(value=42.0)
+        >>> decl2.variable.value == 42.0
+        True
+
+        """
+        kw = self.kwargs()
+        kw.update(kwargs)
+        return Declaration(self.func(**kw))
+
+    def _relation(self, rhs, op):
+        try:
+            rhs = _sympify(rhs)
+        except SympifyError:
+            raise TypeError("Invalid comparison %s < %s" % (self, rhs))
+        return op(self, rhs, evaluate=False)
+
+    __lt__ = lambda self, other: self._relation(other, Lt)
+    __le__ = lambda self, other: self._relation(other, Le)
+    __ge__ = lambda self, other: self._relation(other, Ge)
+    __gt__ = lambda self, other: self._relation(other, Gt)
+
+class Pointer(Variable):
+    """ Represents a pointer. See ``Variable``.
+
+    Examples
+    ========
+
+    Can create instances of ``Element``:
+
+    >>> from sympy import Symbol
+    >>> from sympy.codegen.ast import Pointer
+    >>> i = Symbol('i', integer=True)
+    >>> p = Pointer('x')
+    >>> p[i+1]
+    Element(x, indices=(i + 1,))
+
+    """
+    __slots__ = ()
+
+    def __getitem__(self, key):
+        try:
+            return Element(self.symbol, key)
+        except TypeError:
+            return Element(self.symbol, (key,))
+
+
+class Element(Token):
+    """ Element in (a possibly N-dimensional) array.
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.ast import Element
+    >>> elem = Element('x', 'ijk')
+    >>> elem.symbol.name == 'x'
+    True
+    >>> elem.indices
+    (i, j, k)
+    >>> from sympy import ccode
+    >>> ccode(elem)
+    'x[i][j][k]'
+    >>> ccode(Element('x', 'ijk', strides='lmn', offset='o'))
+    'x[i*l + j*m + k*n + o]'
+
+    """
+    __slots__ = _fields = ('symbol', 'indices', 'strides', 'offset')
+    defaults = {'strides': none, 'offset': none}
+    _construct_symbol = staticmethod(sympify)
+    _construct_indices = staticmethod(lambda arg: Tuple(*arg))
+    _construct_strides = staticmethod(lambda arg: Tuple(*arg))
+    _construct_offset = staticmethod(sympify)
+
+
+class Declaration(Token):
+    """ Represents a variable declaration
+
+    Parameters
+    ==========
+
+    variable : Variable
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.ast import Declaration, NoneToken, untyped
+    >>> z = Declaration('z')
+    >>> z.variable.type == untyped
+    True
+    >>> # value is special NoneToken() which must be tested with == operator
+    >>> z.variable.value is None  # won't work
+    False
+    >>> z.variable.value == None  # not PEP-8 compliant
+    True
+    >>> z.variable.value == NoneToken()  # OK
+    True
+    """
+    __slots__ = _fields = ('variable',)
+    _construct_variable = Variable
+
+
+class While(Token):
+    """ Represents a 'for-loop' in the code.
+
+    Expressions are of the form:
+        "while condition:
+             body..."
+
+    Parameters
+    ==========
+
+    condition : expression convertible to Boolean
+    body : CodeBlock or iterable
+        When passed an iterable it is used to instantiate a CodeBlock.
+
+    Examples
+    ========
+
+    >>> from sympy import symbols, Gt, Abs
+    >>> from sympy.codegen import aug_assign, Assignment, While
+    >>> x, dx = symbols('x dx')
+    >>> expr = 1 - x**2
+    >>> whl = While(Gt(Abs(dx), 1e-9), [
+    ...     Assignment(dx, -expr/expr.diff(x)),
+    ...     aug_assign(x, '+', dx)
+    ... ])
+
+    """
+    __slots__ = _fields = ('condition', 'body')
+    _construct_condition = staticmethod(lambda cond: _sympify(cond))
+
+    @classmethod
+    def _construct_body(cls, itr):
+        if isinstance(itr, CodeBlock):
+            return itr
+        else:
+            return CodeBlock(*itr)
+
+
+class Scope(Token):
+    """ Represents a scope in the code.
+
+    Parameters
+    ==========
+
+    body : CodeBlock or iterable
+        When passed an iterable it is used to instantiate a CodeBlock.
+
+    """
+    __slots__ = _fields = ('body',)
+
+    @classmethod
+    def _construct_body(cls, itr):
+        if isinstance(itr, CodeBlock):
+            return itr
+        else:
+            return CodeBlock(*itr)
+
+
+class Stream(Token):
+    """ Represents a stream.
+
+    There are two predefined Stream instances ``stdout`` & ``stderr``.
+
+    Parameters
+    ==========
+
+    name : str
+
+    Examples
+    ========
+
+    >>> from sympy import pycode, Symbol
+    >>> from sympy.codegen.ast import Print, stderr, QuotedString
+    >>> print(pycode(Print(['x'], file=stderr)))
+    print(x, file=sys.stderr)
+    >>> x = Symbol('x')
+    >>> print(pycode(Print([QuotedString('x')], file=stderr)))  # print literally "x"
+    print("x", file=sys.stderr)
+
+    """
+    __slots__ = _fields = ('name',)
+    _construct_name = String
+
+stdout = Stream('stdout')
+stderr = Stream('stderr')
+
+
+class Print(Token):
+    r""" Represents print command in the code.
+
+    Parameters
+    ==========
+
+    formatstring : str
+    *args : Basic instances (or convertible to such through sympify)
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.ast import Print
+    >>> from sympy import pycode
+    >>> print(pycode(Print('x y'.split(), "coordinate: %12.5g %12.5g\\n")))
+    print("coordinate: %12.5g %12.5g\n" % (x, y), end="")
+
+    """
+
+    __slots__ = _fields = ('print_args', 'format_string', 'file')
+    defaults = {'format_string': none, 'file': none}
+
+    _construct_print_args = staticmethod(_mk_Tuple)
+    _construct_format_string = QuotedString
+    _construct_file = Stream
+
+
+class FunctionPrototype(Node):
+    """ Represents a function prototype
+
+    Allows the user to generate forward declaration in e.g. C/C++.
+
+    Parameters
+    ==========
+
+    return_type : Type
+    name : str
+    parameters: iterable of Variable instances
+    attrs : iterable of Attribute instances
+
+    Examples
+    ========
+
+    >>> from sympy import ccode, symbols
+    >>> from sympy.codegen.ast import real, FunctionPrototype
+    >>> x, y = symbols('x y', real=True)
+    >>> fp = FunctionPrototype(real, 'foo', [x, y])
+    >>> ccode(fp)
+    'double foo(double x, double y)'
+
+    """
+
+    __slots__ = ('return_type', 'name', 'parameters')
+    _fields: tuple[str, ...] = __slots__ + Node._fields
+
+    _construct_return_type = Type
+    _construct_name = String
+
+    @staticmethod
+    def _construct_parameters(args):
+        def _var(arg):
+            if isinstance(arg, Declaration):
+                return arg.variable
+            elif isinstance(arg, Variable):
+                return arg
+            else:
+                return Variable.deduced(arg)
+        return Tuple(*map(_var, args))
+
+    @classmethod
+    def from_FunctionDefinition(cls, func_def):
+        if not isinstance(func_def, FunctionDefinition):
+            raise TypeError("func_def is not an instance of FunctionDefinition")
+        return cls(**func_def.kwargs(exclude=('body',)))
+
+
+class FunctionDefinition(FunctionPrototype):
+    """ Represents a function definition in the code.
+
+    Parameters
+    ==========
+
+    return_type : Type
+    name : str
+    parameters: iterable of Variable instances
+    body : CodeBlock or iterable
+    attrs : iterable of Attribute instances
+
+    Examples
+    ========
+
+    >>> from sympy import ccode, symbols
+    >>> from sympy.codegen.ast import real, FunctionPrototype
+    >>> x, y = symbols('x y', real=True)
+    >>> fp = FunctionPrototype(real, 'foo', [x, y])
+    >>> ccode(fp)
+    'double foo(double x, double y)'
+    >>> from sympy.codegen.ast import FunctionDefinition, Return
+    >>> body = [Return(x*y)]
+    >>> fd = FunctionDefinition.from_FunctionPrototype(fp, body)
+    >>> print(ccode(fd))
+    double foo(double x, double y){
+        return x*y;
+    }
+    """
+
+    __slots__ = ('body', )
+    _fields = FunctionPrototype._fields[:-1] + __slots__ + Node._fields
+
+    @classmethod
+    def _construct_body(cls, itr):
+        if isinstance(itr, CodeBlock):
+            return itr
+        else:
+            return CodeBlock(*itr)
+
+    @classmethod
+    def from_FunctionPrototype(cls, func_proto, body):
+        if not isinstance(func_proto, FunctionPrototype):
+            raise TypeError("func_proto is not an instance of FunctionPrototype")
+        return cls(body=body, **func_proto.kwargs())
+
+
+class Return(Token):
+    """ Represents a return command in the code.
+
+    Parameters
+    ==========
+
+    return : Basic
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.ast import Return
+    >>> from sympy.printing.pycode import pycode
+    >>> from sympy import Symbol
+    >>> x = Symbol('x')
+    >>> print(pycode(Return(x)))
+    return x
+
+    """
+    __slots__ = _fields = ('return',)
+    _construct_return=staticmethod(_sympify)
+
+
+class FunctionCall(Token, Expr):
+    """ Represents a call to a function in the code.
+
+    Parameters
+    ==========
+
+    name : str
+    function_args : Tuple
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.ast import FunctionCall
+    >>> from sympy import pycode
+    >>> fcall = FunctionCall('foo', 'bar baz'.split())
+    >>> print(pycode(fcall))
+    foo(bar, baz)
+
+    """
+    __slots__ = _fields = ('name', 'function_args')
+
+    _construct_name = String
+    _construct_function_args = staticmethod(lambda args: Tuple(*args))
+
+
+class Raise(Token):
+    """ Prints as 'raise ...' in Python, 'throw ...' in C++"""
+    __slots__ = _fields = ('exception',)
+
+
+class RuntimeError_(Token):
+    """ Represents 'std::runtime_error' in C++ and 'RuntimeError' in Python.
+
+    Note that the latter is uncommon, and you might want to use e.g. ValueError.
+    """
+    __slots__ = _fields = ('message',)
+    _construct_message = String

lib/python3.1/site-packages/sympy/codegen/cfunctions.py

@@ -0,0 +1,536 @@
+"""
+This module contains SymPy functions mathcin corresponding to special math functions in the
+C standard library (since C99, also available in C++11).
+
+The functions defined in this module allows the user to express functions such as ``expm1``
+as a SymPy function for symbolic manipulation.
+
+"""
+from sympy.core.function import ArgumentIndexError, Function
+from sympy.core.numbers import Rational
+from sympy.core.power import Pow
+from sympy.core.singleton import S
+from sympy.functions.elementary.exponential import exp, log
+from sympy.functions.elementary.miscellaneous import sqrt
+
+
+def _expm1(x):
+    return exp(x) - S.One
+
+
+class expm1(Function):
+    """
+    Represents the exponential function minus one.
+
+    Explanation
+    ===========
+
+    The benefit of using ``expm1(x)`` over ``exp(x) - 1``
+    is that the latter is prone to cancellation under finite precision
+    arithmetic when x is close to zero.
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x
+    >>> from sympy.codegen.cfunctions import expm1
+    >>> '%.0e' % expm1(1e-99).evalf()
+    '1e-99'
+    >>> from math import exp
+    >>> exp(1e-99) - 1
+    0.0
+    >>> expm1(x).diff(x)
+    exp(x)
+
+    See Also
+    ========
+
+    log1p
+    """
+    nargs = 1
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex == 1:
+            return exp(*self.args)
+        else:
+            raise ArgumentIndexError(self, argindex)
+
+    def _eval_expand_func(self, **hints):
+        return _expm1(*self.args)
+
+    def _eval_rewrite_as_exp(self, arg, **kwargs):
+        return exp(arg) - S.One
+
+    _eval_rewrite_as_tractable = _eval_rewrite_as_exp
+
+    @classmethod
+    def eval(cls, arg):
+        exp_arg = exp.eval(arg)
+        if exp_arg is not None:
+            return exp_arg - S.One
+
+    def _eval_is_real(self):
+        return self.args[0].is_real
+
+    def _eval_is_finite(self):
+        return self.args[0].is_finite
+
+
+def _log1p(x):
+    return log(x + S.One)
+
+
+class log1p(Function):
+    """
+    Represents the natural logarithm of a number plus one.
+
+    Explanation
+    ===========
+
+    The benefit of using ``log1p(x)`` over ``log(x + 1)``
+    is that the latter is prone to cancellation under finite precision
+    arithmetic when x is close to zero.
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x
+    >>> from sympy.codegen.cfunctions import log1p
+    >>> from sympy import expand_log
+    >>> '%.0e' % expand_log(log1p(1e-99)).evalf()
+    '1e-99'
+    >>> from math import log
+    >>> log(1 + 1e-99)
+    0.0
+    >>> log1p(x).diff(x)
+    1/(x + 1)
+
+    See Also
+    ========
+
+    expm1
+    """
+    nargs = 1
+
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex == 1:
+            return S.One/(self.args[0] + S.One)
+        else:
+            raise ArgumentIndexError(self, argindex)
+
+
+    def _eval_expand_func(self, **hints):
+        return _log1p(*self.args)
+
+    def _eval_rewrite_as_log(self, arg, **kwargs):
+        return _log1p(arg)
+
+    _eval_rewrite_as_tractable = _eval_rewrite_as_log
+
+    @classmethod
+    def eval(cls, arg):
+        if arg.is_Rational:
+            return log(arg + S.One)
+        elif not arg.is_Float:  # not safe to add 1 to Float
+            return log.eval(arg + S.One)
+        elif arg.is_number:
+            return log(Rational(arg) + S.One)
+
+    def _eval_is_real(self):
+        return (self.args[0] + S.One).is_nonnegative
+
+    def _eval_is_finite(self):
+        if (self.args[0] + S.One).is_zero:
+            return False
+        return self.args[0].is_finite
+
+    def _eval_is_positive(self):
+        return self.args[0].is_positive
+
+    def _eval_is_zero(self):
+        return self.args[0].is_zero
+
+    def _eval_is_nonnegative(self):
+        return self.args[0].is_nonnegative
+
+_Two = S(2)
+
+def _exp2(x):
+    return Pow(_Two, x)
+
+class exp2(Function):
+    """
+    Represents the exponential function with base two.
+
+    Explanation
+    ===========
+
+    The benefit of using ``exp2(x)`` over ``2**x``
+    is that the latter is not as efficient under finite precision
+    arithmetic.
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x
+    >>> from sympy.codegen.cfunctions import exp2
+    >>> exp2(2).evalf() == 4.0
+    True
+    >>> exp2(x).diff(x)
+    log(2)*exp2(x)
+
+    See Also
+    ========
+
+    log2
+    """
+    nargs = 1
+
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex == 1:
+            return self*log(_Two)
+        else:
+            raise ArgumentIndexError(self, argindex)
+
+    def _eval_rewrite_as_Pow(self, arg, **kwargs):
+        return _exp2(arg)
+
+    _eval_rewrite_as_tractable = _eval_rewrite_as_Pow
+
+    def _eval_expand_func(self, **hints):
+        return _exp2(*self.args)
+
+    @classmethod
+    def eval(cls, arg):
+        if arg.is_number:
+            return _exp2(arg)
+
+
+def _log2(x):
+    return log(x)/log(_Two)
+
+
+class log2(Function):
+    """
+    Represents the logarithm function with base two.
+
+    Explanation
+    ===========
+
+    The benefit of using ``log2(x)`` over ``log(x)/log(2)``
+    is that the latter is not as efficient under finite precision
+    arithmetic.
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x
+    >>> from sympy.codegen.cfunctions import log2
+    >>> log2(4).evalf() == 2.0
+    True
+    >>> log2(x).diff(x)
+    1/(x*log(2))
+
+    See Also
+    ========
+
+    exp2
+    log10
+    """
+    nargs = 1
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex == 1:
+            return S.One/(log(_Two)*self.args[0])
+        else:
+            raise ArgumentIndexError(self, argindex)
+
+
+    @classmethod
+    def eval(cls, arg):
+        if arg.is_number:
+            result = log.eval(arg, base=_Two)
+            if result.is_Atom:
+                return result
+        elif arg.is_Pow and arg.base == _Two:
+            return arg.exp
+
+    def _eval_evalf(self, *args, **kwargs):
+        return self.rewrite(log).evalf(*args, **kwargs)
+
+    def _eval_expand_func(self, **hints):
+        return _log2(*self.args)
+
+    def _eval_rewrite_as_log(self, arg, **kwargs):
+        return _log2(arg)
+
+    _eval_rewrite_as_tractable = _eval_rewrite_as_log
+
+
+def _fma(x, y, z):
+    return x*y + z
+
+
+class fma(Function):
+    """
+    Represents "fused multiply add".
+
+    Explanation
+    ===========
+
+    The benefit of using ``fma(x, y, z)`` over ``x*y + z``
+    is that, under finite precision arithmetic, the former is
+    supported by special instructions on some CPUs.
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x, y, z
+    >>> from sympy.codegen.cfunctions import fma
+    >>> fma(x, y, z).diff(x)
+    y
+
+    """
+    nargs = 3
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex in (1, 2):
+            return self.args[2 - argindex]
+        elif argindex == 3:
+            return S.One
+        else:
+            raise ArgumentIndexError(self, argindex)
+
+
+    def _eval_expand_func(self, **hints):
+        return _fma(*self.args)
+
+    def _eval_rewrite_as_tractable(self, arg, limitvar=None, **kwargs):
+        return _fma(arg)
+
+
+_Ten = S(10)
+
+
+def _log10(x):
+    return log(x)/log(_Ten)
+
+
+class log10(Function):
+    """
+    Represents the logarithm function with base ten.
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x
+    >>> from sympy.codegen.cfunctions import log10
+    >>> log10(100).evalf() == 2.0
+    True
+    >>> log10(x).diff(x)
+    1/(x*log(10))
+
+    See Also
+    ========
+
+    log2
+    """
+    nargs = 1
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex == 1:
+            return S.One/(log(_Ten)*self.args[0])
+        else:
+            raise ArgumentIndexError(self, argindex)
+
+
+    @classmethod
+    def eval(cls, arg):
+        if arg.is_number:
+            result = log.eval(arg, base=_Ten)
+            if result.is_Atom:
+                return result
+        elif arg.is_Pow and arg.base == _Ten:
+            return arg.exp
+
+    def _eval_expand_func(self, **hints):
+        return _log10(*self.args)
+
+    def _eval_rewrite_as_log(self, arg, **kwargs):
+        return _log10(arg)
+
+    _eval_rewrite_as_tractable = _eval_rewrite_as_log
+
+
+def _Sqrt(x):
+    return Pow(x, S.Half)
+
+
+class Sqrt(Function):  # 'sqrt' already defined in sympy.functions.elementary.miscellaneous
+    """
+    Represents the square root function.
+
+    Explanation
+    ===========
+
+    The reason why one would use ``Sqrt(x)`` over ``sqrt(x)``
+    is that the latter is internally represented as ``Pow(x, S.Half)`` which
+    may not be what one wants when doing code-generation.
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x
+    >>> from sympy.codegen.cfunctions import Sqrt
+    >>> Sqrt(x)
+    Sqrt(x)
+    >>> Sqrt(x).diff(x)
+    1/(2*sqrt(x))
+
+    See Also
+    ========
+
+    Cbrt
+    """
+    nargs = 1
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex == 1:
+            return Pow(self.args[0], Rational(-1, 2))/_Two
+        else:
+            raise ArgumentIndexError(self, argindex)
+
+    def _eval_expand_func(self, **hints):
+        return _Sqrt(*self.args)
+
+    def _eval_rewrite_as_Pow(self, arg, **kwargs):
+        return _Sqrt(arg)
+
+    _eval_rewrite_as_tractable = _eval_rewrite_as_Pow
+
+
+def _Cbrt(x):
+    return Pow(x, Rational(1, 3))
+
+
+class Cbrt(Function):  # 'cbrt' already defined in sympy.functions.elementary.miscellaneous
+    """
+    Represents the cube root function.
+
+    Explanation
+    ===========
+
+    The reason why one would use ``Cbrt(x)`` over ``cbrt(x)``
+    is that the latter is internally represented as ``Pow(x, Rational(1, 3))`` which
+    may not be what one wants when doing code-generation.
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x
+    >>> from sympy.codegen.cfunctions import Cbrt
+    >>> Cbrt(x)
+    Cbrt(x)
+    >>> Cbrt(x).diff(x)
+    1/(3*x**(2/3))
+
+    See Also
+    ========
+
+    Sqrt
+    """
+    nargs = 1
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex == 1:
+            return Pow(self.args[0], Rational(-_Two/3))/3
+        else:
+            raise ArgumentIndexError(self, argindex)
+
+
+    def _eval_expand_func(self, **hints):
+        return _Cbrt(*self.args)
+
+    def _eval_rewrite_as_Pow(self, arg, **kwargs):
+        return _Cbrt(arg)
+
+    _eval_rewrite_as_tractable = _eval_rewrite_as_Pow
+
+
+def _hypot(x, y):
+    return sqrt(Pow(x, 2) + Pow(y, 2))
+
+
+class hypot(Function):
+    """
+    Represents the hypotenuse function.
+
+    Explanation
+    ===========
+
+    The hypotenuse function is provided by e.g. the math library
+    in the C99 standard, hence one may want to represent the function
+    symbolically when doing code-generation.
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x, y
+    >>> from sympy.codegen.cfunctions import hypot
+    >>> hypot(3, 4).evalf() == 5.0
+    True
+    >>> hypot(x, y)
+    hypot(x, y)
+    >>> hypot(x, y).diff(x)
+    x/hypot(x, y)
+
+    """
+    nargs = 2
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex in (1, 2):
+            return 2*self.args[argindex-1]/(_Two*self.func(*self.args))
+        else:
+            raise ArgumentIndexError(self, argindex)
+
+
+    def _eval_expand_func(self, **hints):
+        return _hypot(*self.args)
+
+    def _eval_rewrite_as_Pow(self, arg, **kwargs):
+        return _hypot(arg)
+
+    _eval_rewrite_as_tractable = _eval_rewrite_as_Pow
+
+
+class isnan(Function):
+    nargs = 1

lib/python3.1/site-packages/sympy/codegen/cnodes.py

@@ -0,0 +1,156 @@
+"""
+AST nodes specific to the C family of languages
+"""
+
+from sympy.codegen.ast import (
+    Attribute, Declaration, Node, String, Token, Type, none,
+    FunctionCall, CodeBlock
+    )
+from sympy.core.basic import Basic
+from sympy.core.containers import Tuple
+from sympy.core.sympify import sympify
+
+void = Type('void')
+
+restrict = Attribute('restrict')  # guarantees no pointer aliasing
+volatile = Attribute('volatile')
+static = Attribute('static')
+
+
+def alignof(arg):
+    """ Generate of FunctionCall instance for calling 'alignof' """
+    return FunctionCall('alignof', [String(arg) if isinstance(arg, str) else arg])
+
+
+def sizeof(arg):
+    """ Generate of FunctionCall instance for calling 'sizeof'
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.ast import real
+    >>> from sympy.codegen.cnodes import sizeof
+    >>> from sympy import ccode
+    >>> ccode(sizeof(real))
+    'sizeof(double)'
+    """
+    return FunctionCall('sizeof', [String(arg) if isinstance(arg, str) else arg])
+
+
+class CommaOperator(Basic):
+    """ Represents the comma operator in C """
+    def __new__(cls, *args):
+        return Basic.__new__(cls, *[sympify(arg) for arg in args])
+
+
+class Label(Node):
+    """ Label for use with e.g. goto statement.
+
+    Examples
+    ========
+
+    >>> from sympy import ccode, Symbol
+    >>> from sympy.codegen.cnodes import Label, PreIncrement
+    >>> print(ccode(Label('foo')))
+    foo:
+    >>> print(ccode(Label('bar', [PreIncrement(Symbol('a'))])))
+    bar:
+    ++(a);
+
+    """
+    __slots__ = _fields = ('name', 'body')
+    defaults = {'body': none}
+    _construct_name = String
+
+    @classmethod
+    def _construct_body(cls, itr):
+        if isinstance(itr, CodeBlock):
+            return itr
+        else:
+            return CodeBlock(*itr)
+
+
+class goto(Token):
+    """ Represents goto in C """
+    __slots__ = _fields = ('label',)
+    _construct_label = Label
+
+
+class PreDecrement(Basic):
+    """ Represents the pre-decrement operator
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x
+    >>> from sympy.codegen.cnodes import PreDecrement
+    >>> from sympy import ccode
+    >>> ccode(PreDecrement(x))
+    '--(x)'
+
+    """
+    nargs = 1
+
+
+class PostDecrement(Basic):
+    """ Represents the post-decrement operator
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x
+    >>> from sympy.codegen.cnodes import PostDecrement
+    >>> from sympy import ccode
+    >>> ccode(PostDecrement(x))
+    '(x)--'
+
+    """
+    nargs = 1
+
+
+class PreIncrement(Basic):
+    """ Represents the pre-increment operator
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x
+    >>> from sympy.codegen.cnodes import PreIncrement
+    >>> from sympy import ccode
+    >>> ccode(PreIncrement(x))
+    '++(x)'
+
+    """
+    nargs = 1
+
+
+class PostIncrement(Basic):
+    """ Represents the post-increment operator
+
+    Examples
+    ========
+
+    >>> from sympy.abc import x
+    >>> from sympy.codegen.cnodes import PostIncrement
+    >>> from sympy import ccode
+    >>> ccode(PostIncrement(x))
+    '(x)++'
+
+    """
+    nargs = 1
+
+
+class struct(Node):
+    """ Represents a struct in C """
+    __slots__ = _fields = ('name', 'declarations')
+    defaults = {'name': none}
+    _construct_name = String
+
+    @classmethod
+    def _construct_declarations(cls, args):
+        return Tuple(*[Declaration(arg) for arg in args])
+
+
+class union(struct):
+    """ Represents a union in C """
+    __slots__ = ()

lib/python3.1/site-packages/sympy/codegen/cutils.py

@@ -0,0 +1,8 @@
+from sympy.printing.c import C99CodePrinter
+
+def render_as_source_file(content, Printer=C99CodePrinter, settings=None):
+    """ Renders a C source file (with required #include statements) """
+    printer = Printer(settings or {})
+    code_str = printer.doprint(content)
+    includes = '\n'.join(['#include <%s>' % h for h in printer.headers])
+    return includes + '\n\n' + code_str

lib/python3.1/site-packages/sympy/codegen/cxxnodes.py

@@ -0,0 +1,14 @@
+"""
+AST nodes specific to C++.
+"""
+
+from sympy.codegen.ast import Attribute, String, Token, Type, none
+
+class using(Token):
+    """ Represents a 'using' statement in C++ """
+    __slots__ = _fields = ('type', 'alias')
+    defaults = {'alias': none}
+    _construct_type = Type
+    _construct_alias = String
+
+constexpr = Attribute('constexpr')

lib/python3.1/site-packages/sympy/codegen/fnodes.py

@@ -0,0 +1,657 @@
+"""
+AST nodes specific to Fortran.
+
+The functions defined in this module allows the user to express functions such as ``dsign``
+as a SymPy function for symbolic manipulation.
+"""
+
+from sympy.codegen.ast import (
+    Attribute, CodeBlock, FunctionCall, Node, none, String,
+    Token, _mk_Tuple, Variable
+)
+from sympy.core.basic import Basic
+from sympy.core.containers import Tuple
+from sympy.core.expr import Expr
+from sympy.core.function import Function
+from sympy.core.numbers import Float, Integer
+from sympy.core.symbol import Str
+from sympy.core.sympify import sympify
+from sympy.logic import true, false
+from sympy.utilities.iterables import iterable
+
+
+
+pure = Attribute('pure')
+elemental = Attribute('elemental')  # (all elemental procedures are also pure)
+
+intent_in = Attribute('intent_in')
+intent_out = Attribute('intent_out')
+intent_inout = Attribute('intent_inout')
+
+allocatable = Attribute('allocatable')
+
+class Program(Token):
+    """ Represents a 'program' block in Fortran.
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.ast import Print
+    >>> from sympy.codegen.fnodes import Program
+    >>> prog = Program('myprogram', [Print([42])])
+    >>> from sympy import fcode
+    >>> print(fcode(prog, source_format='free'))
+    program myprogram
+        print *, 42
+    end program
+
+    """
+    __slots__ = _fields = ('name', 'body')
+    _construct_name = String
+    _construct_body = staticmethod(lambda body: CodeBlock(*body))
+
+
+class use_rename(Token):
+    """ Represents a renaming in a use statement in Fortran.
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.fnodes import use_rename, use
+    >>> from sympy import fcode
+    >>> ren = use_rename("thingy", "convolution2d")
+    >>> print(fcode(ren, source_format='free'))
+    thingy => convolution2d
+    >>> full = use('signallib', only=['snr', ren])
+    >>> print(fcode(full, source_format='free'))
+    use signallib, only: snr, thingy => convolution2d
+
+    """
+    __slots__ = _fields = ('local', 'original')
+    _construct_local = String
+    _construct_original = String
+
+def _name(arg):
+    if hasattr(arg, 'name'):
+        return arg.name
+    else:
+        return String(arg)
+
+class use(Token):
+    """ Represents a use statement in Fortran.
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.fnodes import use
+    >>> from sympy import fcode
+    >>> fcode(use('signallib'), source_format='free')
+    'use signallib'
+    >>> fcode(use('signallib', [('metric', 'snr')]), source_format='free')
+    'use signallib, metric => snr'
+    >>> fcode(use('signallib', only=['snr', 'convolution2d']), source_format='free')
+    'use signallib, only: snr, convolution2d'
+
+    """
+    __slots__ = _fields = ('namespace', 'rename', 'only')
+    defaults = {'rename': none, 'only': none}
+    _construct_namespace = staticmethod(_name)
+    _construct_rename = staticmethod(lambda args: Tuple(*[arg if isinstance(arg, use_rename) else use_rename(*arg) for arg in args]))
+    _construct_only = staticmethod(lambda args: Tuple(*[arg if isinstance(arg, use_rename) else _name(arg) for arg in args]))
+
+
+class Module(Token):
+    """ Represents a module in Fortran.
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.fnodes import Module
+    >>> from sympy import fcode
+    >>> print(fcode(Module('signallib', ['implicit none'], []), source_format='free'))
+    module signallib
+    implicit none
+    <BLANKLINE>
+    contains
+    <BLANKLINE>
+    <BLANKLINE>
+    end module
+
+    """
+    __slots__ = _fields = ('name', 'declarations', 'definitions')
+    defaults = {'declarations': Tuple()}
+    _construct_name = String
+
+    @classmethod
+    def _construct_declarations(cls, args):
+        args = [Str(arg) if isinstance(arg, str) else arg for arg in args]
+        return CodeBlock(*args)
+
+    _construct_definitions = staticmethod(lambda arg: CodeBlock(*arg))
+
+
+class Subroutine(Node):
+    """ Represents a subroutine in Fortran.
+
+    Examples
+    ========
+
+    >>> from sympy import fcode, symbols
+    >>> from sympy.codegen.ast import Print
+    >>> from sympy.codegen.fnodes import Subroutine
+    >>> x, y = symbols('x y', real=True)
+    >>> sub = Subroutine('mysub', [x, y], [Print([x**2 + y**2, x*y])])
+    >>> print(fcode(sub, source_format='free', standard=2003))
+    subroutine mysub(x, y)
+    real*8 :: x
+    real*8 :: y
+    print *, x**2 + y**2, x*y
+    end subroutine
+
+    """
+    __slots__ = ('name', 'parameters', 'body')
+    _fields = __slots__ + Node._fields
+    _construct_name = String
+    _construct_parameters = staticmethod(lambda params: Tuple(*map(Variable.deduced, params)))
+
+    @classmethod
+    def _construct_body(cls, itr):
+        if isinstance(itr, CodeBlock):
+            return itr
+        else:
+            return CodeBlock(*itr)
+
+class SubroutineCall(Token):
+    """ Represents a call to a subroutine in Fortran.
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.fnodes import SubroutineCall
+    >>> from sympy import fcode
+    >>> fcode(SubroutineCall('mysub', 'x y'.split()))
+    '       call mysub(x, y)'
+
+    """
+    __slots__ = _fields = ('name', 'subroutine_args')
+    _construct_name = staticmethod(_name)
+    _construct_subroutine_args = staticmethod(_mk_Tuple)
+
+
+class Do(Token):
+    """ Represents a Do loop in in Fortran.
+
+    Examples
+    ========
+
+    >>> from sympy import fcode, symbols
+    >>> from sympy.codegen.ast import aug_assign, Print
+    >>> from sympy.codegen.fnodes import Do
+    >>> i, n = symbols('i n', integer=True)
+    >>> r = symbols('r', real=True)
+    >>> body = [aug_assign(r, '+', 1/i), Print([i, r])]
+    >>> do1 = Do(body, i, 1, n)
+    >>> print(fcode(do1, source_format='free'))
+    do i = 1, n
+        r = r + 1d0/i
+        print *, i, r
+    end do
+    >>> do2 = Do(body, i, 1, n, 2)
+    >>> print(fcode(do2, source_format='free'))
+    do i = 1, n, 2
+        r = r + 1d0/i
+        print *, i, r
+    end do
+
+    """
+
+    __slots__ = _fields = ('body', 'counter', 'first', 'last', 'step', 'concurrent')
+    defaults = {'step': Integer(1), 'concurrent': false}
+    _construct_body = staticmethod(lambda body: CodeBlock(*body))
+    _construct_counter = staticmethod(sympify)
+    _construct_first = staticmethod(sympify)
+    _construct_last = staticmethod(sympify)
+    _construct_step = staticmethod(sympify)
+    _construct_concurrent = staticmethod(lambda arg: true if arg else false)
+
+
+class ArrayConstructor(Token):
+    """ Represents an array constructor.
+
+    Examples
+    ========
+
+    >>> from sympy import fcode
+    >>> from sympy.codegen.fnodes import ArrayConstructor
+    >>> ac = ArrayConstructor([1, 2, 3])
+    >>> fcode(ac, standard=95, source_format='free')
+    '(/1, 2, 3/)'
+    >>> fcode(ac, standard=2003, source_format='free')
+    '[1, 2, 3]'
+
+    """
+    __slots__ = _fields = ('elements',)
+    _construct_elements = staticmethod(_mk_Tuple)
+
+
+class ImpliedDoLoop(Token):
+    """ Represents an implied do loop in Fortran.
+
+    Examples
+    ========
+
+    >>> from sympy import Symbol, fcode
+    >>> from sympy.codegen.fnodes import ImpliedDoLoop, ArrayConstructor
+    >>> i = Symbol('i', integer=True)
+    >>> idl = ImpliedDoLoop(i**3, i, -3, 3, 2)  # -27, -1, 1, 27
+    >>> ac = ArrayConstructor([-28, idl, 28]) # -28, -27, -1, 1, 27, 28
+    >>> fcode(ac, standard=2003, source_format='free')
+    '[-28, (i**3, i = -3, 3, 2), 28]'
+
+    """
+    __slots__ = _fields = ('expr', 'counter', 'first', 'last', 'step')
+    defaults = {'step': Integer(1)}
+    _construct_expr = staticmethod(sympify)
+    _construct_counter = staticmethod(sympify)
+    _construct_first = staticmethod(sympify)
+    _construct_last = staticmethod(sympify)
+    _construct_step = staticmethod(sympify)
+
+
+class Extent(Basic):
+    """ Represents a dimension extent.
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.fnodes import Extent
+    >>> e = Extent(-3, 3)  # -3, -2, -1, 0, 1, 2, 3
+    >>> from sympy import fcode
+    >>> fcode(e, source_format='free')
+    '-3:3'
+    >>> from sympy.codegen.ast import Variable, real
+    >>> from sympy.codegen.fnodes import dimension, intent_out
+    >>> dim = dimension(e, e)
+    >>> arr = Variable('x', real, attrs=[dim, intent_out])
+    >>> fcode(arr.as_Declaration(), source_format='free', standard=2003)
+    'real*8, dimension(-3:3, -3:3), intent(out) :: x'
+
+    """
+    def __new__(cls, *args):
+        if len(args) == 2:
+            low, high = args
+            return Basic.__new__(cls, sympify(low), sympify(high))
+        elif len(args) == 0 or (len(args) == 1 and args[0] in (':', None)):
+            return Basic.__new__(cls)  # assumed shape
+        else:
+            raise ValueError("Expected 0 or 2 args (or one argument == None or ':')")
+
+    def _sympystr(self, printer):
+        if len(self.args) == 0:
+            return ':'
+        return ":".join(str(arg) for arg in self.args)
+
+assumed_extent = Extent() # or Extent(':'), Extent(None)
+
+
+def dimension(*args):
+    """ Creates a 'dimension' Attribute with (up to 7) extents.
+
+    Examples
+    ========
+
+    >>> from sympy import fcode
+    >>> from sympy.codegen.fnodes import dimension, intent_in
+    >>> dim = dimension('2', ':')  # 2 rows, runtime determined number of columns
+    >>> from sympy.codegen.ast import Variable, integer
+    >>> arr = Variable('a', integer, attrs=[dim, intent_in])
+    >>> fcode(arr.as_Declaration(), source_format='free', standard=2003)
+    'integer*4, dimension(2, :), intent(in) :: a'
+
+    """
+    if len(args) > 7:
+        raise ValueError("Fortran only supports up to 7 dimensional arrays")
+    parameters = []
+    for arg in args:
+        if isinstance(arg, Extent):
+            parameters.append(arg)
+        elif isinstance(arg, str):
+            if arg == ':':
+                parameters.append(Extent())
+            else:
+                parameters.append(String(arg))
+        elif iterable(arg):
+            parameters.append(Extent(*arg))
+        else:
+            parameters.append(sympify(arg))
+    if len(args) == 0:
+        raise ValueError("Need at least one dimension")
+    return Attribute('dimension', parameters)
+
+
+assumed_size = dimension('*')
+
+def array(symbol, dim, intent=None, *, attrs=(), value=None, type=None):
+    """ Convenience function for creating a Variable instance for a Fortran array.
+
+    Parameters
+    ==========
+
+    symbol : symbol
+    dim : Attribute or iterable
+        If dim is an ``Attribute`` it need to have the name 'dimension'. If it is
+        not an ``Attribute``, then it is passed to :func:`dimension` as ``*dim``
+    intent : str
+        One of: 'in', 'out', 'inout' or None
+    \\*\\*kwargs:
+        Keyword arguments for ``Variable`` ('type' & 'value')
+
+    Examples
+    ========
+
+    >>> from sympy import fcode
+    >>> from sympy.codegen.ast import integer, real
+    >>> from sympy.codegen.fnodes import array
+    >>> arr = array('a', '*', 'in', type=integer)
+    >>> print(fcode(arr.as_Declaration(), source_format='free', standard=2003))
+    integer*4, dimension(*), intent(in) :: a
+    >>> x = array('x', [3, ':', ':'], intent='out', type=real)
+    >>> print(fcode(x.as_Declaration(value=1), source_format='free', standard=2003))
+    real*8, dimension(3, :, :), intent(out) :: x = 1
+
+    """
+    if isinstance(dim, Attribute):
+        if str(dim.name) != 'dimension':
+            raise ValueError("Got an unexpected Attribute argument as dim: %s" % str(dim))
+    else:
+        dim = dimension(*dim)
+
+    attrs = list(attrs) + [dim]
+    if intent is not None:
+        if intent not in (intent_in, intent_out, intent_inout):
+            intent = {'in': intent_in, 'out': intent_out, 'inout': intent_inout}[intent]
+        attrs.append(intent)
+    if type is None:
+        return Variable.deduced(symbol, value=value, attrs=attrs)
+    else:
+        return Variable(symbol, type, value=value, attrs=attrs)
+
+def _printable(arg):
+    return String(arg) if isinstance(arg, str) else sympify(arg)
+
+
+def allocated(array):
+    """ Creates an AST node for a function call to Fortran's "allocated(...)"
+
+    Examples
+    ========
+
+    >>> from sympy import fcode
+    >>> from sympy.codegen.fnodes import allocated
+    >>> alloc = allocated('x')
+    >>> fcode(alloc, source_format='free')
+    'allocated(x)'
+
+    """
+    return FunctionCall('allocated', [_printable(array)])
+
+
+def lbound(array, dim=None, kind=None):
+    """ Creates an AST node for a function call to Fortran's "lbound(...)"
+
+    Parameters
+    ==========
+
+    array : Symbol or String
+    dim : expr
+    kind : expr
+
+    Examples
+    ========
+
+    >>> from sympy import fcode
+    >>> from sympy.codegen.fnodes import lbound
+    >>> lb = lbound('arr', dim=2)
+    >>> fcode(lb, source_format='free')
+    'lbound(arr, 2)'
+
+    """
+    return FunctionCall(
+        'lbound',
+        [_printable(array)] +
+        ([_printable(dim)] if dim else []) +
+        ([_printable(kind)] if kind else [])
+    )
+
+
+def ubound(array, dim=None, kind=None):
+    return FunctionCall(
+        'ubound',
+        [_printable(array)] +
+        ([_printable(dim)] if dim else []) +
+        ([_printable(kind)] if kind else [])
+    )
+
+
+def shape(source, kind=None):
+    """ Creates an AST node for a function call to Fortran's "shape(...)"
+
+    Parameters
+    ==========
+
+    source : Symbol or String
+    kind : expr
+
+    Examples
+    ========
+
+    >>> from sympy import fcode
+    >>> from sympy.codegen.fnodes import shape
+    >>> shp = shape('x')
+    >>> fcode(shp, source_format='free')
+    'shape(x)'
+
+    """
+    return FunctionCall(
+        'shape',
+        [_printable(source)] +
+        ([_printable(kind)] if kind else [])
+    )
+
+
+def size(array, dim=None, kind=None):
+    """ Creates an AST node for a function call to Fortran's "size(...)"
+
+    Examples
+    ========
+
+    >>> from sympy import fcode, Symbol
+    >>> from sympy.codegen.ast import FunctionDefinition, real, Return
+    >>> from sympy.codegen.fnodes import array, sum_, size
+    >>> a = Symbol('a', real=True)
+    >>> body = [Return((sum_(a**2)/size(a))**.5)]
+    >>> arr = array(a, dim=[':'], intent='in')
+    >>> fd = FunctionDefinition(real, 'rms', [arr], body)
+    >>> print(fcode(fd, source_format='free', standard=2003))
+    real*8 function rms(a)
+    real*8, dimension(:), intent(in) :: a
+    rms = sqrt(sum(a**2)*1d0/size(a))
+    end function
+
+    """
+    return FunctionCall(
+        'size',
+        [_printable(array)] +
+        ([_printable(dim)] if dim else []) +
+        ([_printable(kind)] if kind else [])
+    )
+
+
+def reshape(source, shape, pad=None, order=None):
+    """ Creates an AST node for a function call to Fortran's "reshape(...)"
+
+    Parameters
+    ==========
+
+    source : Symbol or String
+    shape : ArrayExpr
+
+    """
+    return FunctionCall(
+        'reshape',
+        [_printable(source), _printable(shape)] +
+        ([_printable(pad)] if pad else []) +
+        ([_printable(order)] if pad else [])
+    )
+
+
+def bind_C(name=None):
+    """ Creates an Attribute ``bind_C`` with a name.
+
+    Parameters
+    ==========
+
+    name : str
+
+    Examples
+    ========
+
+    >>> from sympy import fcode, Symbol
+    >>> from sympy.codegen.ast import FunctionDefinition, real, Return
+    >>> from sympy.codegen.fnodes import array, sum_, bind_C
+    >>> a = Symbol('a', real=True)
+    >>> s = Symbol('s', integer=True)
+    >>> arr = array(a, dim=[s], intent='in')
+    >>> body = [Return((sum_(a**2)/s)**.5)]
+    >>> fd = FunctionDefinition(real, 'rms', [arr, s], body, attrs=[bind_C('rms')])
+    >>> print(fcode(fd, source_format='free', standard=2003))
+    real*8 function rms(a, s) bind(C, name="rms")
+    real*8, dimension(s), intent(in) :: a
+    integer*4 :: s
+    rms = sqrt(sum(a**2)/s)
+    end function
+
+    """
+    return Attribute('bind_C', [String(name)] if name else [])
+
+class GoTo(Token):
+    """ Represents a goto statement in Fortran
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.fnodes import GoTo
+    >>> go = GoTo([10, 20, 30], 'i')
+    >>> from sympy import fcode
+    >>> fcode(go, source_format='free')
+    'go to (10, 20, 30), i'
+
+    """
+    __slots__ = _fields = ('labels', 'expr')
+    defaults = {'expr': none}
+    _construct_labels = staticmethod(_mk_Tuple)
+    _construct_expr = staticmethod(sympify)
+
+
+class FortranReturn(Token):
+    """ AST node explicitly mapped to a fortran "return".
+
+    Explanation
+    ===========
+
+    Because a return statement in fortran is different from C, and
+    in order to aid reuse of our codegen ASTs the ordinary
+    ``.codegen.ast.Return`` is interpreted as assignment to
+    the result variable of the function. If one for some reason needs
+    to generate a fortran RETURN statement, this node should be used.
+
+    Examples
+    ========
+
+    >>> from sympy.codegen.fnodes import FortranReturn
+    >>> from sympy import fcode
+    >>> fcode(FortranReturn('x'))
+    '       return x'
+
+    """
+    __slots__ = _fields = ('return_value',)
+    defaults = {'return_value': none}
+    _construct_return_value = staticmethod(sympify)
+
+
+class FFunction(Function):
+    _required_standard = 77
+
+    def _fcode(self, printer):
+        name = self.__class__.__name__
+        if printer._settings['standard'] < self._required_standard:
+            raise NotImplementedError("%s requires Fortran %d or newer" %
+                                      (name, self._required_standard))
+        return '{}({})'.format(name, ', '.join(map(printer._print, self.args)))
+
+
+class F95Function(FFunction):
+    _required_standard = 95
+
+
+class isign(FFunction):
+    """ Fortran sign intrinsic for integer arguments. """
+    nargs = 2
+
+
+class dsign(FFunction):
+    """ Fortran sign intrinsic for double precision arguments. """
+    nargs = 2
+
+
+class cmplx(FFunction):
+    """ Fortran complex conversion function. """
+    nargs = 2  # may be extended to (2, 3) at a later point
+
+
+class kind(FFunction):
+    """ Fortran kind function. """
+    nargs = 1
+
+
+class merge(F95Function):
+    """ Fortran merge function """
+    nargs = 3
+
+
+class _literal(Float):
+    _token = None  # type: str
+    _decimals = None  # type: int
+
+    def _fcode(self, printer, *args, **kwargs):
+        mantissa, sgnd_ex = ('%.{}e'.format(self._decimals) % self).split('e')
+        mantissa = mantissa.strip('0').rstrip('.')
+        ex_sgn, ex_num = sgnd_ex[0], sgnd_ex[1:].lstrip('0')
+        ex_sgn = '' if ex_sgn == '+' else ex_sgn
+        return (mantissa or '0') + self._token + ex_sgn + (ex_num or '0')
+
+
+class literal_sp(_literal):
+    """ Fortran single precision real literal """
+    _token = 'e'
+    _decimals = 9
+
+
+class literal_dp(_literal):
+    """ Fortran double precision real literal """
+    _token = 'd'
+    _decimals = 17
+
+
+class sum_(Token, Expr):
+    __slots__ = _fields = ('array', 'dim', 'mask')
+    defaults = {'dim': none, 'mask': none}
+    _construct_array = staticmethod(sympify)
+    _construct_dim = staticmethod(sympify)
+
+
+class product_(Token, Expr):
+    __slots__ = _fields = ('array', 'dim', 'mask')
+    defaults = {'dim': none, 'mask': none}
+    _construct_array = staticmethod(sympify)
+    _construct_dim = staticmethod(sympify)

lib/python3.1/site-packages/sympy/codegen/futils.py

@@ -0,0 +1,40 @@
+from itertools import chain
+from sympy.codegen.fnodes import Module
+from sympy.core.symbol import Dummy
+from sympy.printing.fortran import FCodePrinter
+
+""" This module collects utilities for rendering Fortran code. """
+
+
+def render_as_module(definitions, name, declarations=(), printer_settings=None):
+    """ Creates a ``Module`` instance and renders it as a string.
+
+    This generates Fortran source code for a module with the correct ``use`` statements.
+
+    Parameters
+    ==========
+
+    definitions : iterable
+        Passed to :class:`sympy.codegen.fnodes.Module`.
+    name : str
+        Passed to :class:`sympy.codegen.fnodes.Module`.
+    declarations : iterable
+        Passed to :class:`sympy.codegen.fnodes.Module`. It will be extended with
+        use statements, 'implicit none' and public list generated from ``definitions``.
+    printer_settings : dict
+        Passed to ``FCodePrinter`` (default: ``{'standard': 2003, 'source_format': 'free'}``).
+
+    """
+    printer_settings = printer_settings or {'standard': 2003, 'source_format': 'free'}
+    printer = FCodePrinter(printer_settings)
+    dummy = Dummy()
+    if isinstance(definitions, Module):
+        raise ValueError("This function expects to construct a module on its own.")
+    mod = Module(name, chain(declarations, [dummy]), definitions)
+    fstr = printer.doprint(mod)
+    module_use_str = '   %s\n' % '   \n'.join(['use %s, only: %s' % (k, ', '.join(v)) for
+                                                k, v in printer.module_uses.items()])
+    module_use_str += '   implicit none\n'
+    module_use_str += '   private\n'
+    module_use_str += '   public %s\n' % ', '.join([str(node.name) for node in definitions if getattr(node, 'name', None)])
+    return fstr.replace(printer.doprint(dummy), module_use_str)

lib/python3.1/site-packages/sympy/codegen/matrix_nodes.py

@@ -0,0 +1,71 @@
+"""
+Additional AST nodes for operations on matrices. The nodes in this module
+are meant to represent optimization of matrix expressions within codegen's
+target languages that cannot be represented by SymPy expressions.
+
+As an example, we can use :meth:`sympy.codegen.rewriting.optimize` and the
+``matin_opt`` optimization provided in :mod:`sympy.codegen.rewriting` to
+transform matrix multiplication under certain assumptions:
+
+    >>> from sympy import symbols, MatrixSymbol
+    >>> n = symbols('n', integer=True)
+    >>> A = MatrixSymbol('A', n, n)
+    >>> x = MatrixSymbol('x', n, 1)
+    >>> expr = A**(-1) * x
+    >>> from sympy import assuming, Q
+    >>> from sympy.codegen.rewriting import matinv_opt, optimize
+    >>> with assuming(Q.fullrank(A)):
+    ...     optimize(expr, [matinv_opt])
+    MatrixSolve(A, vector=x)
+"""
+
+from .ast import Token
+from sympy.matrices import MatrixExpr
+from sympy.core.sympify import sympify
+
+
+class MatrixSolve(Token, MatrixExpr):
+    """Represents an operation to solve a linear matrix equation.
+
+    Parameters
+    ==========
+
+    matrix : MatrixSymbol
+
+      Matrix representing the coefficients of variables in the linear
+      equation. This matrix must be square and full-rank (i.e. all columns must
+      be linearly independent) for the solving operation to be valid.
+
+    vector : MatrixSymbol
+
+      One-column matrix representing the solutions to the equations
+      represented in ``matrix``.
+
+    Examples
+    ========
+
+    >>> from sympy import symbols, MatrixSymbol
+    >>> from sympy.codegen.matrix_nodes import MatrixSolve
+    >>> n = symbols('n', integer=True)
+    >>> A = MatrixSymbol('A', n, n)
+    >>> x = MatrixSymbol('x', n, 1)
+    >>> from sympy.printing.numpy import NumPyPrinter
+    >>> NumPyPrinter().doprint(MatrixSolve(A, x))
+    'numpy.linalg.solve(A, x)'
+    >>> from sympy import octave_code
+    >>> octave_code(MatrixSolve(A, x))
+    'A \\\\ x'
+
+    """
+    __slots__ = _fields = ('matrix', 'vector')
+
+    _construct_matrix = staticmethod(sympify)
+    _construct_vector = staticmethod(sympify)
+
+    @property
+    def shape(self):
+        return self.vector.shape
+
+    def _eval_derivative(self, x):
+        A, b = self.matrix, self.vector
+        return MatrixSolve(A, b.diff(x) - A.diff(x) * MatrixSolve(A, b))

lib/python3.1/site-packages/sympy/codegen/numpy_nodes.py

@@ -0,0 +1,110 @@
+from sympy.core.function import Add, ArgumentIndexError, Function
+from sympy.core.power import Pow
+from sympy.core.singleton import S
+from sympy.core.sorting import default_sort_key
+from sympy.functions.elementary.exponential import exp, log
+
+
+def _logaddexp(x1, x2, *, evaluate=True):
+    return log(Add(exp(x1, evaluate=evaluate), exp(x2, evaluate=evaluate), evaluate=evaluate))
+
+
+_two = S.One*2
+_ln2 = log(_two)
+
+
+def _lb(x, *, evaluate=True):
+    return log(x, evaluate=evaluate)/_ln2
+
+
+def _exp2(x, *, evaluate=True):
+    return Pow(_two, x, evaluate=evaluate)
+
+
+def _logaddexp2(x1, x2, *, evaluate=True):
+    return _lb(Add(_exp2(x1, evaluate=evaluate),
+                   _exp2(x2, evaluate=evaluate), evaluate=evaluate))
+
+
+class logaddexp(Function):
+    """ Logarithm of the sum of exponentiations of the inputs.
+
+    Helper class for use with e.g. numpy.logaddexp
+
+    See Also
+    ========
+
+    https://numpy.org/doc/stable/reference/generated/numpy.logaddexp.html
+    """
+    nargs = 2
+
+    def __new__(cls, *args):
+        return Function.__new__(cls, *sorted(args, key=default_sort_key))
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex == 1:
+            wrt, other = self.args
+        elif argindex == 2:
+            other, wrt = self.args
+        else:
+            raise ArgumentIndexError(self, argindex)
+        return S.One/(S.One + exp(other-wrt))
+
+    def _eval_rewrite_as_log(self, x1, x2, **kwargs):
+        return _logaddexp(x1, x2)
+
+    def _eval_evalf(self, *args, **kwargs):
+        return self.rewrite(log).evalf(*args, **kwargs)
+
+    def _eval_simplify(self, *args, **kwargs):
+        a, b = (x.simplify(**kwargs) for x in self.args)
+        candidate = _logaddexp(a, b)
+        if candidate != _logaddexp(a, b, evaluate=False):
+            return candidate
+        else:
+            return logaddexp(a, b)
+
+
+class logaddexp2(Function):
+    """ Logarithm of the sum of exponentiations of the inputs in base-2.
+
+    Helper class for use with e.g. numpy.logaddexp2
+
+    See Also
+    ========
+
+    https://numpy.org/doc/stable/reference/generated/numpy.logaddexp2.html
+    """
+    nargs = 2
+
+    def __new__(cls, *args):
+        return Function.__new__(cls, *sorted(args, key=default_sort_key))
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex == 1:
+            wrt, other = self.args
+        elif argindex == 2:
+            other, wrt = self.args
+        else:
+            raise ArgumentIndexError(self, argindex)
+        return S.One/(S.One + _exp2(other-wrt))
+
+    def _eval_rewrite_as_log(self, x1, x2, **kwargs):
+        return _logaddexp2(x1, x2)
+
+    def _eval_evalf(self, *args, **kwargs):
+        return self.rewrite(log).evalf(*args, **kwargs)
+
+    def _eval_simplify(self, *args, **kwargs):
+        a, b = (x.simplify(**kwargs).factor() for x in self.args)
+        candidate = _logaddexp2(a, b)
+        if candidate != _logaddexp2(a, b, evaluate=False):
+            return candidate
+        else:
+            return logaddexp2(a, b)

lib/python3.1/site-packages/sympy/codegen/pynodes.py

@@ -0,0 +1,11 @@
+from .abstract_nodes import List as AbstractList
+from .ast import Token
+
+
+class List(AbstractList):
+    pass
+
+
+class NumExprEvaluate(Token):
+    """represents a call to :class:`numexpr`s :func:`evaluate`"""
+    __slots__ = _fields = ('expr',)

lib/python3.1/site-packages/sympy/codegen/pyutils.py

@@ -0,0 +1,24 @@
+from sympy.printing.pycode import PythonCodePrinter
+
+""" This module collects utilities for rendering Python code. """
+
+
+def render_as_module(content, standard='python3'):
+    """Renders Python code as a module (with the required imports).
+
+    Parameters
+    ==========
+
+    standard :
+        See the parameter ``standard`` in
+        :meth:`sympy.printing.pycode.pycode`
+    """
+
+    printer = PythonCodePrinter({'standard':standard})
+    pystr = printer.doprint(content)
+    if printer._settings['fully_qualified_modules']:
+        module_imports_str = '\n'.join('import %s' % k for k in printer.module_imports)
+    else:
+        module_imports_str = '\n'.join(['from %s import %s' % (k, ', '.join(v)) for
+                                        k, v in printer.module_imports.items()])
+    return module_imports_str + '\n\n' + pystr

lib/python3.1/site-packages/sympy/codegen/rewriting.py

@@ -0,0 +1,357 @@
+"""
+Classes and functions useful for rewriting expressions for optimized code
+generation. Some languages (or standards thereof), e.g. C99, offer specialized
+math functions for better performance and/or precision.
+
+Using the ``optimize`` function in this module, together with a collection of
+rules (represented as instances of ``Optimization``), one can rewrite the
+expressions for this purpose::
+
+    >>> from sympy import Symbol, exp, log
+    >>> from sympy.codegen.rewriting import optimize, optims_c99
+    >>> x = Symbol('x')
+    >>> optimize(3*exp(2*x) - 3, optims_c99)
+    3*expm1(2*x)
+    >>> optimize(exp(2*x) - 1 - exp(-33), optims_c99)
+    expm1(2*x) - exp(-33)
+    >>> optimize(log(3*x + 3), optims_c99)
+    log1p(x) + log(3)
+    >>> optimize(log(2*x + 3), optims_c99)
+    log(2*x + 3)
+
+The ``optims_c99`` imported above is tuple containing the following instances
+(which may be imported from ``sympy.codegen.rewriting``):
+
+- ``expm1_opt``
+- ``log1p_opt``
+- ``exp2_opt``
+- ``log2_opt``
+- ``log2const_opt``
+
+
+"""
+from sympy.core.function import expand_log
+from sympy.core.singleton import S
+from sympy.core.symbol import Wild
+from sympy.functions.elementary.complexes import sign
+from sympy.functions.elementary.exponential import (exp, log)
+from sympy.functions.elementary.miscellaneous import (Max, Min)
+from sympy.functions.elementary.trigonometric import (cos, sin, sinc)
+from sympy.assumptions import Q, ask
+from sympy.codegen.cfunctions import log1p, log2, exp2, expm1
+from sympy.codegen.matrix_nodes import MatrixSolve
+from sympy.core.expr import UnevaluatedExpr
+from sympy.core.power import Pow
+from sympy.codegen.numpy_nodes import logaddexp, logaddexp2
+from sympy.codegen.scipy_nodes import cosm1, powm1
+from sympy.core.mul import Mul
+from sympy.matrices.expressions.matexpr import MatrixSymbol
+from sympy.utilities.iterables import sift
+
+
+class Optimization:
+    """ Abstract base class for rewriting optimization.
+
+    Subclasses should implement ``__call__`` taking an expression
+    as argument.
+
+    Parameters
+    ==========
+    cost_function : callable returning number
+    priority : number
+
+    """
+    def __init__(self, cost_function=None, priority=1):
+        self.cost_function = cost_function
+        self.priority=priority
+
+    def cheapest(self, *args):
+        return min(args, key=self.cost_function)
+
+
+class ReplaceOptim(Optimization):
+    """ Rewriting optimization calling replace on expressions.
+
+    Explanation
+    ===========
+
+    The instance can be used as a function on expressions for which
+    it will apply the ``replace`` method (see
+    :meth:`sympy.core.basic.Basic.replace`).
+
+    Parameters
+    ==========
+
+    query :
+        First argument passed to replace.
+    value :
+        Second argument passed to replace.
+
+    Examples
+    ========
+
+    >>> from sympy import Symbol
+    >>> from sympy.codegen.rewriting import ReplaceOptim
+    >>> from sympy.codegen.cfunctions import exp2
+    >>> x = Symbol('x')
+    >>> exp2_opt = ReplaceOptim(lambda p: p.is_Pow and p.base == 2,
+    ...     lambda p: exp2(p.exp))
+    >>> exp2_opt(2**x)
+    exp2(x)
+
+    """
+
+    def __init__(self, query, value, **kwargs):
+        super().__init__(**kwargs)
+        self.query = query
+        self.value = value
+
+    def __call__(self, expr):
+        return expr.replace(self.query, self.value)
+
+
+def optimize(expr, optimizations):
+    """ Apply optimizations to an expression.
+
+    Parameters
+    ==========
+
+    expr : expression
+    optimizations : iterable of ``Optimization`` instances
+        The optimizations will be sorted with respect to ``priority`` (highest first).
+
+    Examples
+    ========
+
+    >>> from sympy import log, Symbol
+    >>> from sympy.codegen.rewriting import optims_c99, optimize
+    >>> x = Symbol('x')
+    >>> optimize(log(x+3)/log(2) + log(x**2 + 1), optims_c99)
+    log1p(x**2) + log2(x + 3)
+
+    """
+
+    for optim in sorted(optimizations, key=lambda opt: opt.priority, reverse=True):
+        new_expr = optim(expr)
+        if optim.cost_function is None:
+            expr = new_expr
+        else:
+            expr = optim.cheapest(expr, new_expr)
+    return expr
+
+
+exp2_opt = ReplaceOptim(
+    lambda p: p.is_Pow and p.base == 2,
+    lambda p: exp2(p.exp)
+)
+
+
+_d = Wild('d', properties=[lambda x: x.is_Dummy])
+_u = Wild('u', properties=[lambda x: not x.is_number and not x.is_Add])
+_v = Wild('v')
+_w = Wild('w')
+_n = Wild('n', properties=[lambda x: x.is_number])
+
+sinc_opt1 = ReplaceOptim(
+    sin(_w)/_w, sinc(_w)
+)
+sinc_opt2 = ReplaceOptim(
+    sin(_n*_w)/_w, _n*sinc(_n*_w)
+)
+sinc_opts = (sinc_opt1, sinc_opt2)
+
+log2_opt = ReplaceOptim(_v*log(_w)/log(2), _v*log2(_w), cost_function=lambda expr: expr.count(
+    lambda e: (  # division & eval of transcendentals are expensive floating point operations...
+        e.is_Pow and e.exp.is_negative  # division
+        or (isinstance(e, (log, log2)) and not e.args[0].is_number))  # transcendental
+    )
+)
+
+log2const_opt = ReplaceOptim(log(2)*log2(_w), log(_w))
+
+logsumexp_2terms_opt = ReplaceOptim(
+    lambda l: (isinstance(l, log)
+               and l.args[0].is_Add
+               and len(l.args[0].args) == 2
+               and all(isinstance(t, exp) for t in l.args[0].args)),
+    lambda l: (
+        Max(*[e.args[0] for e in l.args[0].args]) +
+        log1p(exp(Min(*[e.args[0] for e in l.args[0].args])))
+    )
+)
+
+
+class FuncMinusOneOptim(ReplaceOptim):
+    """Specialization of ReplaceOptim for functions evaluating "f(x) - 1".
+
+    Explanation
+    ===========
+
+    Numerical functions which go toward one as x go toward zero is often best
+    implemented by a dedicated function in order to avoid catastrophic
+    cancellation. One such example is ``expm1(x)`` in the C standard library
+    which evaluates ``exp(x) - 1``. Such functions preserves many more
+    significant digits when its argument is much smaller than one, compared
+    to subtracting one afterwards.
+
+    Parameters
+    ==========
+
+    func :
+        The function which is subtracted by one.
+    func_m_1 :
+        The specialized function evaluating ``func(x) - 1``.
+    opportunistic : bool
+        When ``True``, apply the transformation as long as the magnitude of the
+        remaining number terms decreases. When ``False``, only apply the
+        transformation if it completely eliminates the number term.
+
+    Examples
+    ========
+
+    >>> from sympy import symbols, exp
+    >>> from sympy.codegen.rewriting import FuncMinusOneOptim
+    >>> from sympy.codegen.cfunctions import expm1
+    >>> x, y = symbols('x y')
+    >>> expm1_opt = FuncMinusOneOptim(exp, expm1)
+    >>> expm1_opt(exp(x) + 2*exp(5*y) - 3)
+    expm1(x) + 2*expm1(5*y)
+
+
+    """
+
+    def __init__(self, func, func_m_1, opportunistic=True):
+        weight = 10  # <-- this is an arbitrary number (heuristic)
+        super().__init__(lambda e: e.is_Add, self.replace_in_Add,
+                         cost_function=lambda expr: expr.count_ops() - weight*expr.count(func_m_1))
+        self.func = func
+        self.func_m_1 = func_m_1
+        self.opportunistic = opportunistic
+
+    def _group_Add_terms(self, add):
+        numbers, non_num = sift(add.args, lambda arg: arg.is_number, binary=True)
+        numsum = sum(numbers)
+        terms_with_func, other = sift(non_num, lambda arg: arg.has(self.func), binary=True)
+        return numsum, terms_with_func, other
+
+    def replace_in_Add(self, e):
+        """ passed as second argument to Basic.replace(...) """
+        numsum, terms_with_func, other_non_num_terms = self._group_Add_terms(e)
+        if numsum == 0:
+            return e
+        substituted, untouched = [], []
+        for with_func in terms_with_func:
+            if with_func.is_Mul:
+                func, coeff = sift(with_func.args, lambda arg: arg.func == self.func, binary=True)
+                if len(func) == 1 and len(coeff) == 1:
+                    func, coeff = func[0], coeff[0]
+                else:
+                    coeff = None
+            elif with_func.func == self.func:
+                func, coeff = with_func, S.One
+            else:
+                coeff = None
+
+            if coeff is not None and coeff.is_number and sign(coeff) == -sign(numsum):
+                if self.opportunistic:
+                    do_substitute = abs(coeff+numsum) < abs(numsum)
+                else:
+                    do_substitute = coeff+numsum == 0
+
+                if do_substitute:  # advantageous substitution
+                    numsum += coeff
+                    substituted.append(coeff*self.func_m_1(*func.args))
+                    continue
+            untouched.append(with_func)
+
+        return e.func(numsum, *substituted, *untouched, *other_non_num_terms)
+
+    def __call__(self, expr):
+        alt1 = super().__call__(expr)
+        alt2 = super().__call__(expr.factor())
+        return self.cheapest(alt1, alt2)
+
+
+expm1_opt = FuncMinusOneOptim(exp, expm1)
+cosm1_opt = FuncMinusOneOptim(cos, cosm1)
+powm1_opt = FuncMinusOneOptim(Pow, powm1)
+
+log1p_opt = ReplaceOptim(
+    lambda e: isinstance(e, log),
+    lambda l: expand_log(l.replace(
+        log, lambda arg: log(arg.factor())
+    )).replace(log(_u+1), log1p(_u))
+)
+
+def create_expand_pow_optimization(limit, *, base_req=lambda b: b.is_symbol):
+    """ Creates an instance of :class:`ReplaceOptim` for expanding ``Pow``.
+
+    Explanation
+    ===========
+
+    The requirements for expansions are that the base needs to be a symbol
+    and the exponent needs to be an Integer (and be less than or equal to
+    ``limit``).
+
+    Parameters
+    ==========
+
+    limit : int
+         The highest power which is expanded into multiplication.
+    base_req : function returning bool
+         Requirement on base for expansion to happen, default is to return
+         the ``is_symbol`` attribute of the base.
+
+    Examples
+    ========
+
+    >>> from sympy import Symbol, sin
+    >>> from sympy.codegen.rewriting import create_expand_pow_optimization
+    >>> x = Symbol('x')
+    >>> expand_opt = create_expand_pow_optimization(3)
+    >>> expand_opt(x**5 + x**3)
+    x**5 + x*x*x
+    >>> expand_opt(x**5 + x**3 + sin(x)**3)
+    x**5 + sin(x)**3 + x*x*x
+    >>> opt2 = create_expand_pow_optimization(3, base_req=lambda b: not b.is_Function)
+    >>> opt2((x+1)**2 + sin(x)**2)
+    sin(x)**2 + (x + 1)*(x + 1)
+
+    """
+    return ReplaceOptim(
+        lambda e: e.is_Pow and base_req(e.base) and e.exp.is_Integer and abs(e.exp) <= limit,
+        lambda p: (
+            UnevaluatedExpr(Mul(*([p.base]*+p.exp), evaluate=False)) if p.exp > 0 else
+            1/UnevaluatedExpr(Mul(*([p.base]*-p.exp), evaluate=False))
+        ))
+
+# Optimization procedures for turning A**(-1) * x into MatrixSolve(A, x)
+def _matinv_predicate(expr):
+    # TODO: We should be able to support more than 2 elements
+    if expr.is_MatMul and len(expr.args) == 2:
+        left, right = expr.args
+        if left.is_Inverse and right.shape[1] == 1:
+            inv_arg = left.arg
+            if isinstance(inv_arg, MatrixSymbol):
+                return bool(ask(Q.fullrank(left.arg)))
+
+    return False
+
+def _matinv_transform(expr):
+    left, right = expr.args
+    inv_arg = left.arg
+    return MatrixSolve(inv_arg, right)
+
+
+matinv_opt = ReplaceOptim(_matinv_predicate, _matinv_transform)
+
+
+logaddexp_opt = ReplaceOptim(log(exp(_v)+exp(_w)), logaddexp(_v, _w))
+logaddexp2_opt = ReplaceOptim(log(Pow(2, _v)+Pow(2, _w)), logaddexp2(_v, _w)*log(2))
+
+# Collections of optimizations:
+optims_c99 = (expm1_opt, log1p_opt, exp2_opt, log2_opt, log2const_opt)
+
+optims_numpy = optims_c99 + (logaddexp_opt, logaddexp2_opt,) + sinc_opts
+
+optims_scipy = (cosm1_opt, powm1_opt)

lib/python3.1/site-packages/sympy/codegen/scipy_nodes.py

@@ -0,0 +1,79 @@
+from sympy.core.function import Add, ArgumentIndexError, Function
+from sympy.core.power import Pow
+from sympy.core.singleton import S
+from sympy.functions.elementary.exponential import log
+from sympy.functions.elementary.trigonometric import cos, sin
+
+
+def _cosm1(x, *, evaluate=True):
+    return Add(cos(x, evaluate=evaluate), -S.One, evaluate=evaluate)
+
+
+class cosm1(Function):
+    """ Minus one plus cosine of x, i.e. cos(x) - 1. For use when x is close to zero.
+
+    Helper class for use with e.g. scipy.special.cosm1
+    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.cosm1.html
+    """
+    nargs = 1
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex == 1:
+            return -sin(*self.args)
+        else:
+            raise ArgumentIndexError(self, argindex)
+
+    def _eval_rewrite_as_cos(self, x, **kwargs):
+        return _cosm1(x)
+
+    def _eval_evalf(self, *args, **kwargs):
+        return self.rewrite(cos).evalf(*args, **kwargs)
+
+    def _eval_simplify(self, **kwargs):
+        x, = self.args
+        candidate = _cosm1(x.simplify(**kwargs))
+        if candidate != _cosm1(x, evaluate=False):
+            return candidate
+        else:
+            return cosm1(x)
+
+
+def _powm1(x, y, *, evaluate=True):
+    return Add(Pow(x, y, evaluate=evaluate), -S.One, evaluate=evaluate)
+
+
+class powm1(Function):
+    """ Minus one plus x to the power of y, i.e. x**y - 1. For use when x is close to one or y is close to zero.
+
+    Helper class for use with e.g. scipy.special.powm1
+    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.powm1.html
+    """
+    nargs = 2
+
+    def fdiff(self, argindex=1):
+        """
+        Returns the first derivative of this function.
+        """
+        if argindex == 1:
+            return Pow(self.args[0], self.args[1])*self.args[1]/self.args[0]
+        elif argindex == 2:
+            return log(self.args[0])*Pow(*self.args)
+        else:
+            raise ArgumentIndexError(self, argindex)
+
+    def _eval_rewrite_as_Pow(self, x, y, **kwargs):
+        return _powm1(x, y)
+
+    def _eval_evalf(self, *args, **kwargs):
+        return self.rewrite(Pow).evalf(*args, **kwargs)
+
+    def _eval_simplify(self, **kwargs):
+        x, y = self.args
+        candidate = _powm1(x.simplify(**kwargs), y.simplify(**kwargs))
+        if candidate != _powm1(x, y, evaluate=False):
+            return candidate
+        else:
+            return powm1(x, y)