Split mdpoly.algebra into multiple files

17 files changed, 854 insertions, 731 deletions

diff --git a/mdpoly/__init__.py b/mdpoly/__init__.py
index 15ab8c9..1ac17a5 100644
--- a/mdpoly/__init__.py
+++ b/mdpoly/__init__.py
@@ -61,9 +61,9 @@ Expression
 
 Within the code, we will often call something like :math:`x^2 + 1` and
 *expression* (for example :py:class:`mdpoly.abc.Expr` or
-:py:class:`mdpoly.algebra.PolyRingExpr`) instead of polynomial. Expressions are
-more general and can include any mathematical operation. Of course in the end
-and expression should result in a polynomial.
+:py:class:`mdpoly.expressions.poly.PolyRingExpr`) instead of polynomial.
+Expressions are more general and can include any mathematical operation. Of
+course in the end and expression should result in a polynomial.
 
 
 Decision Variables, Program
@@ -107,12 +107,8 @@ TODO: data structure(s) that represent the polynomial
 
 from .abc import (Shape as _Shape)
 
-from .algebra import (PolyConst as _PolyConst,
-                      PolyVar as _PolyVar,
-                      PolyParam as _PolyParam,
-                      MatConst as _MatConst,
-                      MatVar as _MatVar,
-                      MatParam as _MatParam)
+from .expressions.poly import (PolyConst as _PolyConst, PolyVar as _PolyVar, PolyParam as _PolyParam)
+from .expressions.matrix import (MatConst as _MatConst, MatVar as _MatVar, MatParam as _MatParam)
 
 from .state import State as _State
 
diff --git a/mdpoly/abc.py b/mdpoly/abc.py
index 16b155b..cd90f85 100644
--- a/mdpoly/abc.py
+++ b/mdpoly/abc.py
@@ -2,23 +2,21 @@
 from __future__ import annotations
 from typing import TYPE_CHECKING
 
-from .index import Number, Shape, MatrixIndex, PolyIndex
-from .constants import NUMERICS_EPS
-from .errors import AlgebraicError
-from .util import iszero
-
-if TYPE_CHECKING:
-    from .state import State
-
 from typing import Self, Type, TypeVar, Generic, Any, Iterable, Sequence
 from enum import Enum, auto
 from copy import copy
 from abc import ABC, abstractmethod
-from dataclassabc import dataclassabc
 from hashlib import sha256
+from dataclassabc import dataclassabc
 
+from .constants import NUMERICS_EPS
+from .errors import AlgebraicError
+from .util import iszero
+
+if TYPE_CHECKING:
+    from .index import Number, Shape, MatrixIndex, PolyIndex
+    from .state import State
 
-ReprT = TypeVar("ReprT", bound="Repr")
 
 # ┏━╸╻ ╻┏━┓┏━┓┏━╸┏━┓┏━┓╻┏━┓┏┓╻┏━┓
 # ┣╸ ┏╋┛┣━┛┣┳┛┣╸ ┗━┓┗━┓┃┃ ┃┃┗┫┗━┓
@@ -380,6 +378,9 @@ class Rel(ABC):
 # ╹┗╸┗━╸╹  ╹┗╸┗━╸┗━┛┗━╸╹ ╹ ╹ ╹ ╹ ╹ ╹┗━┛╹ ╹┗━┛
 
 
+ReprT = TypeVar("ReprT", bound="Repr")
+
+
 class Repr(ABC):
     r""" Representation of a multivariate matrix polynomial expression.
 
diff --git a/mdpoly/algebra.py b/mdpoly/algebra.py
deleted file mode 100644
index 03a6fd5..0000000
--- a/mdpoly/algebra.py
+++ /dev/null
@@ -1,701 +0,0 @@
-""" Algebraic Structures for Expressions """
-from __future__ import annotations
-
-from .abc import Algebra, Expr, ReprT, Nothing, Const, Var, Param
-from .errors import AlgebraicError, InvalidShape, MissingParameters
-from .state import State
-from .index import Shape, MatrixIndex, PolyIndex, PolyVarIndex, Number
-
-from typing import cast, Sequence, Iterable, Type, TypeVar, Self
-from functools import reduce
-from itertools import product, chain, combinations_with_replacement
-from abc import abstractmethod
-from dataclasses import field, dataclass
-import operator
-
-from dataclassabc import dataclassabc
-
-
-@dataclassabc(eq=False)
-class BinaryOp(Expr):
-    left: Expr = field()
-    right: Expr = field()
-
-
-@dataclassabc(eq=False)
-class UnaryOp(Expr):
-    right: Expr = field()
-
-    @property
-    def inner(self) -> Expr:
-        """ Inner expression on which the operator is acting, alias for right. """
-        return self.right
-
-    @property
-    def left(self) -> Expr:
-        return Nothing()
-
-    @left.setter
-    def left(self, left) -> None:
-        if not isinstance(left, Nothing):
-            raise ValueError("Cannot set left of left-acting unary operator "
-                             "to something that is not of type Nothing.")
-
-
-class Reducible(Expr):
-    """ Reducible Expression.
-
-    Algebraic expression that can be written in terms of other (existing)
-    expression objects, i.e. can be reduced to another expression made of
-    simpler operations. For example subtraction can be written in term of
-    addition and multiplication.
-    """
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
-        return self.reduce().to_repr(repr_type, state)
-
-    @abstractmethod
-    def reduce(self) -> Expr:
-        """ Reduce the expression to its basic elements """
-
-
-# TODO: review this idea before implementing
-# class Simplifiable(Expr):
-#     """ Simplifiable Expression. """
-# 
-#     @abstractmethod
-#     def simplify(self) -> Expr:
-#         """ Simplify the expression """
-
-
-# ┏━┓┏━┓╻  ╻ ╻┏┓╻┏━┓┏┳┓╻┏━┓╻     ┏━┓╻  ┏━╸┏━╸┏┓ ┏━┓┏━┓
-# ┣━┛┃ ┃┃  ┗┳┛┃┗┫┃ ┃┃┃┃┃┣━┫┃     ┣━┫┃  ┃╺┓┣╸ ┣┻┓┣┳┛┣━┫
-# ╹  ┗━┛┗━╸ ╹ ╹ ╹┗━┛╹ ╹╹╹ ╹┗━╸   ╹ ╹┗━╸┗━┛┗━╸┗━┛╹┗╸╹ ╹
-
-
-class PolyRingExpr(Expr):
-    r""" Expression with the algebraic behaviour of a polynomial ring.
-
-    This is the algebra of :math:`\mathbb{R}[x_1, \ldots, x_n]`. Note that the
-    polynomials are scalars.
-    """
-
-    # -- Properties ---
-
-    @property
-    def algebra(self) -> Algebra:
-        return Algebra.poly_ring
-
-    @property
-    def shape(self) -> Shape:
-        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
-        if self.left.shape != self.right.shape:
-            raise InvalidShape(f"Cannot perform operation {repr(self)} with "
-                                f"shapes {self.left.shape} and {self.right.shape}.")
-        return self.left.shape
-
-    # --- Helper methods for construction ---
-
-    @classmethod
-    def from_powers(cls, variable: PolyRingExpr, exponents: Iterable[int]) -> PolyRingExpr:
-        """ Given a variable, say :math:`x`, and a list of exponents ``[0, 3, 5]``
-        returns a polynomial :math:`x^5 + x^3 + 1`.
-        """
-        nonzero_exponents = filter(lambda e: e != 0, exponents)
-        # Ignore typecheck because dataclassabc has no typing stub
-        constant_term = PolyConst(1 if 0 in exponents else 0) # type: ignore[call-arg]
-        # FIXME: remove annoying 0
-        return sum((variable ** e for e in nonzero_exponents), constant_term)
-
-
-    @classmethod
-    def make_combinations(cls, variables: Iterable[PolyVar], max_degree: int) -> Iterable[PolyRingExpr]:
-        """ Make combinations of terms.
-
-        For example given :math:`x, y` and *max_degree* of 2 generates
-        :math:`x^2, xy, x, y^2, y, 1`.
-        """
-        # Ignore typecheck because dataclassabc has no typing stub
-        vars_and_const = chain(variables, (PolyConst(1),)) # type: ignore[call-arg]
-        for comb in combinations_with_replacement(vars_and_const, max_degree):
-            yield cast(PolyRingExpr, reduce(operator.mul, comb))
-
-    # --- Operator Overloading ---
-
-    def __add__(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, PolyConst, other)
-        return PolyAdd(self, other)
-
-    def __radd__(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, PolyConst, other)
-        return PolyAdd(other, self)
-
-    def __sub__(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, PolyConst, other)
-        return PolySub(self, other)
-
-    def __rsub__(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, PolyConst, other)
-        return PolyAdd(other, self)
-
-    def __neg__(self):
-        # FIXME: Create PolyNeg?
-        return PolyMul(PolyConst(-1), self)
-
-    def __mul__(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, PolyConst, other)
-        return PolyMul(self, other)
-
-    def __rmul__(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, PolyConst, other)
-        return PolyMul(other, self)
-
-    def __matmul__(self, other):
-        raise AlgebraicError("Cannot perform matrix multiplication in polynomial ring (they are scalars).")
-
-    def __rmatmul__(self, other):
-        self.__rmatmul__(other)
-
-    def __truediv__(self, other):
-        other = self._wrap(Number, PolyConst, other)
-        if not isinstance(other, PolyConst | PolyParam):
-            raise AlgebraicError("Cannot divide by variables in polynomial ring.")
-
-        return PolyMul(self, other)
-
-    def __rtruediv__(self, other):
-        raise AlgebraicError("Cannot perform right division in polynomial ring.")
-
-    def __pow__(self, other):
-        other = self._wrap(Number, PolyConst, other)
-        if not isinstance(other, PolyConst | PolyParam):
-             raise AlgebraicError(f"Cannot raise to powers of type {type(other)} in "
-                                  "polynomial ring. Only constants and parameters are allowed.")
-        return PolyExp(left=self, right=other)
-
-    # -- Other mathematical operations ---
-    
-    def diff(self, wrt: PolyVar) -> PolyPartialDiff:
-        return PolyPartialDiff(self, wrt)
-
-
-@dataclassabc(frozen=True)
-class PolyVar(Var, PolyRingExpr):
-    """ Variable TODO: desc """
-    name: str # overloads Leaf.name
-    shape: Shape = Shape.scalar() # ovearloads PolyRingExpr.shape
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        r = repr_type(self.shape)
-        idx = PolyVarIndex.from_var(self, state), # important comma!
-        r.set(MatrixIndex.scalar(), PolyIndex(idx), 1)
-        return r, state
-
-
-@dataclassabc(frozen=True)
-class PolyConst(Const, PolyRingExpr):
-    """ Constant TODO: desc """
-    value: Number # overloads Const.value
-    name: str = "" # overloads Leaf.name
-    shape: Shape = Shape.scalar() # ovearloads PolyRingExpr.shape
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        r = repr_type(self.shape)
-        r.set(MatrixIndex.scalar(), PolyIndex.constant(), self.value)
-        return r, state
-
-
-@dataclassabc(frozen=True)
-class PolyParam(Param, PolyRingExpr):
-    """ Polynomial parameter TODO: desc """
-    name: str # overloads Leaf.name
-    shape: Shape = Shape.scalar() # overloads PolyRingExpr.shape
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        if self not in state.parameters:
-            raise MissingParameters("Cannot construct representation because "
-                                    f"value for parameter {self} was not given.")
-
-        return PolyConst(state.parameters[self]).to_repr(repr_type, state) # type: ignore[call-arg]
-
-
-@dataclass(eq=False)
-class PolyAdd(BinaryOp, PolyRingExpr):
-    """ Addition operator between scalar polynomials. """
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
-        # Make a new empty representation
-        r = repr_type(self.shape)
-
-        # Make representations for existing stuff
-        for node in self.children():
-            nrepr, state = node.to_repr(repr_type, state)
-
-            # Process non-zero terms
-            for entry in nrepr.entries():
-                for term in nrepr.terms(entry):
-                    s = r.at(entry, term) + nrepr.at(entry, term)
-                    r.set(entry, term, s)
-
-        return r, state
-
-    def __str__(self) -> str:
-        return f"({self.left} + {self.right})"
-
-
-@dataclass(eq=False)
-class PolySub(BinaryOp, PolyRingExpr, Reducible):
-    """ Subtraction operator between scalar polynomials. """
-
-    def reduce(self) -> Expr:
-        """ See :py:meth:`mdpoly.algebra.Reducible.reduce`. """
-        return self.left + (-1 * self.right)
-
-    def __str__(self) -> str:
-        return f"({self.left} - {self.right})"
-
-
-@dataclass(eq=False)
-class PolyMul(BinaryOp, PolyRingExpr):
-    """ Multiplication operator between scalar polynomials. """
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
-        r = repr_type(self.shape)
-
-        lrepr, state = self.left.to_repr(repr_type, state)
-        rrepr, state = self.right.to_repr(repr_type, state)
-
-        # Non zero entries are the intersection since if either is zero the
-        # result is zero
-        nonzero_entries = set(lrepr.entries()) & set(rrepr.entries())
-        for entry in nonzero_entries:
-            # Compute polynomial product between non-zero entries
-            for lterm, rterm in product(lrepr.terms(entry), rrepr.terms(entry)):
-                # Compute where the results should go
-                term = PolyIndex.product(lterm, rterm)
-
-                # Compute product
-                p = r.at(entry, term) + lrepr.at(entry, lterm) * rrepr.at(entry, rterm)
-                r.set(entry, term, p)
-
-        return r, state
-
-    def __str__(self) -> str:
-        return f"({self.left} * {self.right})"
-
-
-@dataclass(eq=False)
-class PolyExp(BinaryOp, PolyRingExpr, Reducible):
-    """ Exponentiation operator between scalar polynomials. """
-
-    @property # type: ignore[override]
-    def right(self) -> Const: # type: ignore[override]
-        if not isinstance(super().right, Const):
-            raise AlgebraicError(f"Cannot raise {self.left} to {self.right} because"
-                                 f"{self.right} is not a constant.")
-
-        return cast(Const, super().right)
-
-    @right.setter
-    def right(self, right: Expr) -> None:
-        # Workaround necessary because of a bug
-        # https://bugs.python.org/issue14965
-        BinaryOp.right.fset(self, right) # type: ignore
-
-    def reduce(self) -> Expr:
-        """ See :py:meth:`mdpoly.algebra.Reducible.reduce`. """
-        var = self.left
-        if not isinstance(self.right.value, int):
-            raise NotImplementedError("Cannot raise to non-integer powers (yet).")
-
-        ntimes = self.right.value - 1
-        return reduce(operator.mul, (var for _ in range(ntimes)), var)
-
-    def __str__(self) -> str:
-        return f"({self.left} ** {self.right})"
-
-
-@dataclassabc(eq=False)
-class PolyPartialDiff(UnaryOp, PolyRingExpr):
-    """ Partial differentiation of scalar polynomials. """
-    wrt: PolyVar
-
-    @property
-    def shape(self) -> Shape:
-        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
-        return self.inner.shape
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
-        r = repr_type(self.shape)
-        lrepr, state = self.inner.to_repr(repr_type, state)
-
-        entry = MatrixIndex.scalar()
-        wrt = state.index(self.wrt)
-
-        for term in lrepr.terms(entry):
-            # sadly, walrus operator does not support unpacking
-            # https://bugs.python.org/issue43143
-            # if ((newterm, fac) := PolyIndex.differentiate(term, wrt)):
-            if result := PolyIndex.differentiate(term, wrt):
-                newterm, fac = result
-                r.set(entry, newterm, fac * lrepr.at(entry, term))
-
-        return r, state
-
-    def __str__(self) -> str:
-        return f"(∂_{self.wrt} {self.inner})"
-
-    def replace(self, old: PolyRingExpr, new: PolyRingExpr) -> Self:
-        """ Overloads :py:meth:`mdpoly.abc.Expr.replace` """
-        if self.wrt == old:
-            if not isinstance(new, PolyVar):
-                # FIXME: implement chain rule
-                raise AlgebraicError(f"Cannot take a derivative with respect to {new}.")
-
-            self.wrt = cast(PolyVar, new)
-
-        return cast(Self, Expr.replace(self, old, new))
-
-
-# ┏┳┓┏━┓╺┳╸┏━┓╻╻ ╻   ┏━┓╻  ┏━╸┏━╸┏┓ ┏━┓┏━┓
-# ┃┃┃┣━┫ ┃ ┣┳┛┃┏╋┛   ┣━┫┃  ┃╺┓┣╸ ┣┻┓┣┳┛┣━┫
-# ╹ ╹╹ ╹ ╹ ╹┗╸╹╹ ╹   ╹ ╹┗━╸┗━┛┗━╸┗━┛╹┗╸╹ ╹
-
-class MatrixExpr(Expr):
-    r""" Expression with the algebraic properties of a matrix ring and / or
-    module (depending on the shape).
-
-    We denote with :math:`R` a polynomial ring.
-
-    - If the shape is square, like :math:`(n, n)` then this is :math:`M_n(R)`
-      the ring of matrices over :math:`R`.
-
-    - If the shape is something else like a row or column (:math:`(m, p)`,
-      :math:`(1, n)` or :math:`(n, 1)`) this is a module, i.e. an algebra with
-      addition and scalar multiplication, where the "scalars" come from
-      :math:`R`.
-    
-    Furthermore some operators that are usually expected from matrices are
-    already included (eg. transposition).
-    """
-
-    @property
-    def algebra(self) -> Algebra:
-        return Algebra.matrix_ring
-
-    def at(self, entry: MatrixIndex, scalar_type: type) -> Expr:
-        raise NotImplementedError
-
-    def __add__(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, MatConst, other)
-        return MatAdd(self, other)
-
-    def __sub__(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, MatConst, other)
-        return MatSub(self, other)
-
-    def __rsub__(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, MatConst, other)
-        return MatSub(other, self)
-
-    def __mul__(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, MatConst, other)
-        return MatScalarMul(other, self) 
-
-    def __rmul__(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, MatConst, other)
-        return MatScalarMul(other, self)
-
-    def __matmul__(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, MatConst, other)
-        return MatMul(self, other)
-
-    def __rmatmul(self, other):
-        self._assert_same_algebra(self, other)
-        other = self._wrap(Number, MatConst, other)
-        return MatMul(other, self)
-
-    def __truediv__(self, scalar):
-        scalar = self._wrap_if_constant(scalar)
-        raise NotImplementedError
-
-    def transpose(self) -> MatrixExpr:
-        """ Matrix transposition. """
-        raise NotImplementedError
-        return MatTranspose(self)
-
-    @property
-    def T(self) -> MatrixExpr:
-        """ Shorthand for :py:meth:`mdpoly.algebra.MatrixExpr.transpose`. """
-        return self.transpose()
-
-    def to_scalar(self, scalar_type: type):
-        """ Convert to a scalar expression. """
-        raise NotImplementedError
-
-
-@dataclassabc(frozen=True)
-class MatConst(Const, MatrixExpr):
-    """ Matrix constant. TODO: desc. """
-    value: Sequence[Sequence[Number]] # Row major, overloads Const.value
-    shape: Shape # overloads Expr.shape
-    name: str = "" # overloads Leaf.name
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        r = repr_type(self.shape)
-
-        for i, row in enumerate(self.value):
-            for j, val in enumerate(row):
-                r.set(MatrixIndex(row=i, col=j), PolyIndex.constant(), val)
-
-        return r, state
-
-
-    def __str__(self) -> str:
-        if not self.name:
-            return repr(self.value)
-
-        return self.name
-
-
-
-T = TypeVar("T", bound=Var)
-
-@dataclassabc(frozen=True)
-class MatVar(Var, MatrixExpr):
-    """ Matrix of polynomial variables. TODO: desc """
-    name: str # overloads Leaf.name
-    shape: Shape # overloads Expr.shape
-
-    # TODO: review this API, can be moved elsewhere?
-    def to_scalars(self, scalar_var_type: Type[T]) -> Iterable[tuple[MatrixIndex, T]]:
-        for row in range(self.shape.rows):
-            for col in range(self.shape.cols):
-                var = scalar_var_type(name=f"{self.name}_[{row},{col}]") # type: ignore[call-arg]
-                entry = MatrixIndex(row, col)
-
-                yield entry, var
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        r = repr_type(self.shape)
-        
-        # FIXME: do not hardcode scalar type
-        for entry, var in self.to_scalars(PolyVar):
-            idx = PolyVarIndex.from_var(var, state), # important comma!
-            r.set(entry, PolyIndex(idx), 1)
-
-        return r, state
-
-    def __str__(self) -> str:
-        return self.name
-
-
-@dataclassabc(frozen=True)
-class MatParam(Param, MatrixExpr):
-    """ Matrix parameter. TODO: desc. """
-    name: str
-    shape: Shape
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        if self not in state.parameters:
-            raise MissingParameters("Cannot construct representation because "
-                                    f"value for parameter {self} was not given.")
-
-        # FIXME: add conversion to scalar variables
-        # Ignore typecheck because dataclassabc has not type stub
-        return MatConst(state.parameters[self]).to_repr(repr_type, state) # type: ignore[call-arg]
-
-    def __str__(self) -> str:
-        return self.name
-
-
-@dataclassabc
-class MatAdd(BinaryOp, MatrixExpr):
-    """ Addition operator between matrices. """
-
-    @property
-    def shape(self) -> Shape:
-        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
-        if not self.left.shape == self.right.shape:
-            raise AlgebraicError(f"Cannot add matrices {self.left} and {self.right} "
-                                 "with different shapes.")
-        return self.left.shape
-
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
-        # Make a new empty representation
-        r = repr_type(self.shape)
-
-        # Make representations for existing stuff
-        for node in self.children():
-            nrepr, state = node.to_repr(repr_type, state)
-
-            # Process non-zero terms
-            for entry in nrepr.entries():
-                for term in nrepr.terms(entry):
-                    s = r.at(entry, term) + nrepr.at(entry, term)
-                    r.set(entry, term, s)
-
-        return r, state
-
-    def __str__(self) -> str:
-        return f"({self.left} + {self.right})"
-
-
-@dataclassabc
-class MatSub(BinaryOp, MatrixExpr, Reducible):
-    """ Subtraction operator between matrices. """
-
-    @property
-    def shape(self) -> Shape:
-        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
-        if not self.left.shape == self.right.shape:
-            raise AlgebraicError(f"Cannot subtract matrices {self.left} and {self.right} "
-                                 f"with different shapes, {self.left.shape} and {self.right.shape}.")
-        return self.left.shape
-
-    def reduce(self) -> Expr:
-        """ See :py:meth:`mdpoly.algebra.Reducible.reduce`. """
-        return self.left + (-1 * self.right)
-
-    def __str__(self) -> str:
-        return f"({self.left} - {self.right})"
-
-
-@dataclassabc
-class MatElemMul(BinaryOp, MatrixExpr):
-    """ Elementwise Matrix Multiplication. """
-
-    @property
-    def shape(self) -> Shape:
-        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
-        if not self.left.shape == self.right.shape:
-            raise AlgebraicError("Cannot perform element-wise multiplication of matrices "
-                                 f"{self.left} and {self.right} with different shapes, "
-                                 f"{self.left.shape} and {self.right.shape}")
-        return self.left.shape
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
-        r = repr_type(self.shape)
-
-        lrepr, state = self.left.to_repr(repr_type, state)
-        rrepr, state = self.right.to_repr(repr_type, state)
-
-        # Non zero entries are the intersection since if either is zero the
-        # result is zero
-        nonzero_entries = set(lrepr.entries()) & set(rrepr.entries())
-        for entry in nonzero_entries:
-            # Compute polynomial product between non-zero entries
-            for lterm, rterm in product(lrepr.terms(entry), rrepr.terms(entry)):
-                # Compute where the results should go
-                term = PolyIndex.product(lterm, rterm)
-
-                # Compute product
-                p = r.at(entry, term) + lrepr.at(entry, lterm) * rrepr.at(entry, rterm)
-                r.set(entry, term, p)
-
-        return r, state
-
-    def __str__(self) -> str:
-        return f"({self.left} .* {self.right})"
-
-
-@dataclassabc
-class MatScalarMul(BinaryOp, MatrixExpr):
-    """ Matrix-Scalar Multiplication. Assumes scalar is on the left and matrix
-    on the right. """
-
-    @property
-    def shape(self) -> Shape:
-        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
-        if not self.left.shape == Shape.scalar():
-            raise InvalidShape(f"Matrix-scalar product assumes that left argumet {self.left} "
-                               f"but it has shape {self.left.shape}")
-
-        return self.right.shape
-
-
-    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
-        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
-        r = repr_type(self.shape)
-
-        scalar_repr, state = self.left.to_repr(repr_type, state)
-        mat_repr, state = self.right.to_repr(repr_type, state)
-
-        for entry in mat_repr.entries():
-            scalar_terms = scalar_repr.terms(MatrixIndex.scalar())
-            mat_terms = mat_repr.terms(entry)
-            for scalar_term, mat_term in product(scalar_terms, mat_terms):
-                term = PolyIndex.product(scalar_term, mat_term)
-
-                p = r.at(entry, term) + scalar_repr.at(entry, scalar_term) + mat_repr.at(entry, mat_term)
-                r.set(entry, term, p)
-
-        return r, state
-
-    def __str__(self) -> str:
-        return f"({self.left} * {self.right})"
-
-
-@dataclass(eq=False)
-class MatDotProd(BinaryOp, MatrixExpr):
-    """ Dot product. """
-
-    @property
-    def shape(self) -> Shape:
-        if not self.left.shape.is_row():
-            raise AlgebraicError(f"Left operand {self.left} must be a row!")
-
-        if not self.right.shape.is_col():
-            raise AlgebraicError(f"Right operand {self.right} must be a column!")
-
-        if self.left.shape.cols != self.right.shape.rows:
-            raise AlgebraicError(f"Rows of {self.right} and columns {self.left} do not match!")
-
-        return Shape.scalar()
-
-
-@dataclass(eq=False)
-class MatMul(BinaryOp, MatrixExpr):
-    """ Matrix Multiplication. """
-
-    @property
-    def shape(self) -> Shape:
-        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
-        if not self.left.shape.rows == self.right.shape.cols:
-            raise AlgebraicError("Cannot perform matrix multiplication between "
-                                 f"{self.left} and {self.right} (shapes {self.left.shape} and {self.right.shape})")
-
-        return Shape(self.left.shape.rows, self.right.shape.cols)
-
-
-    def __str__(self) -> str:
-        return f"({self.left} @ {self.right})"
-
-
-@dataclass(eq=False)
-class MatTranspose(UnaryOp, MatrixExpr):
-    """ Matrix transposition """
-
-    @property
-    def shape(self) -> Shape:
-        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
-        return Shape(self.inner.shape.cols, self.inner.shape.rows)
diff --git a/mdpoly/errors.py b/mdpoly/errors.py
index 85432b8..262f27e 100644
--- a/mdpoly/errors.py
+++ b/mdpoly/errors.py
@@ -1,7 +1,6 @@
-
 class AlgebraicError(Exception):
     """ This is raised whenever an invalid mathematical operation is performed.
-    See also :py:mod:`mdpoly.algebra`.
+    See also :py:mod:`mdpoly.operations`.
     """
 
 class InvalidShape(Exception):
diff --git a/mdpoly/expressions/__init__.py b/mdpoly/expressions/__init__.py
new file mode 100644
index 0000000..2efe342
--- /dev/null
+++ b/mdpoly/expressions/__init__.py
@@ -0,0 +1,66 @@
+from __future__ import annotations
+from typing import TYPE_CHECKING
+
+from abc import abstractmethod
+from typing import Type
+from dataclasses import field
+from dataclassabc import dataclassabc
+
+from ..abc import Expr, Nothing
+
+if TYPE_CHECKING:
+    from ..abc import ReprT
+    from ..state import State
+
+
+@dataclassabc(eq=False)
+class BinaryOp(Expr):
+    left: Expr = field()
+    right: Expr = field()
+
+
+@dataclassabc(eq=False)
+class UnaryOp(Expr):
+    right: Expr = field()
+
+    @property
+    def inner(self) -> Expr:
+        """ Inner expression on which the operator is acting, alias for right. """
+        return self.right
+
+    @property
+    def left(self) -> Expr:
+        return Nothing()
+
+    @left.setter
+    def left(self, left) -> None:
+        if not isinstance(left, Nothing):
+            raise ValueError("Cannot set left of left-acting unary operator "
+                             "to something that is not of type Nothing.")
+
+
+class Reducible(Expr):
+    """ Reducible Expression.
+
+    Algebraic expression that can be written in terms of other (existing)
+    expression objects, i.e. can be reduced to another expression made of
+    simpler operations. For example subtraction can be written in term of
+    addition and multiplication.
+    """
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
+        return self.reduce().to_repr(repr_type, state)
+
+    @abstractmethod
+    def reduce(self) -> Expr:
+        """ Reduce the expression to its basic elements """
+
+
+# TODO: review this idea before implementing
+# class Simplifiable(Expr):
+#     """ Simplifiable Expression. """
+# 
+#     @abstractmethod
+#     def simplify(self) -> Expr:
+#         """ Simplify the expression """
diff --git a/mdpoly/expressions/matrix.py b/mdpoly/expressions/matrix.py
new file mode 100644
index 0000000..65f86c2
--- /dev/null
+++ b/mdpoly/expressions/matrix.py
@@ -0,0 +1,175 @@
+from __future__ import annotations
+from typing import TYPE_CHECKING
+
+from typing import Type, TypeVar, Iterable, Sequence
+from dataclassabc import dataclassabc
+
+from .poly import PolyVar, PolyConst
+from ..abc import Expr, Var, Const, Param, Algebra
+from ..index import MatrixIndex, PolyIndex, PolyVarIndex
+from ..errors import MissingParameters
+
+from .. import operations
+
+if TYPE_CHECKING:
+    from ..abc import ReprT
+    from ..index import Shape, Number
+    from ..state import State
+
+
+class MatrixExpr(Expr):
+    r""" Expression with the algebraic properties of a matrix ring and / or
+    module (depending on the shape).
+
+    We denote with :math:`R` a polynomial ring.
+
+    - If the shape is square, like :math:`(n, n)` then this is :math:`M_n(R)`
+      the ring of matrices over :math:`R`.
+
+    - If the shape is something else like a row or column (:math:`(m, p)`,
+      :math:`(1, n)` or :math:`(n, 1)`) this is a module, i.e. an algebra with
+      addition and scalar multiplication, where the "scalars" come from
+      :math:`R`.
+    
+    Furthermore some operators that are usually expected from matrices are
+    already included (eg. transposition).
+    """
+
+    @property
+    def algebra(self) -> Algebra:
+        return Algebra.matrix_ring
+
+    def __add__(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, MatConst, other)
+        return operations.add.MatAdd(self, other)
+
+    def __sub__(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, MatConst, other)
+        return operations.add.MatSub(self, other)
+
+    def __rsub__(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, MatConst, other)
+        return operations.add.MatSub(other, self)
+
+    def __neg__(self):
+        # FIXME: Create PolyNeg?
+        return operations.mul.MatScalarMul(PolyConst(-1), self)
+
+    def __mul__(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, MatConst, other)
+
+        # TODO: case distiction based on shapes
+        return operations.mul.MatScalarMul(other, self) 
+
+
+    def __rmul__(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, MatConst, other)
+        return operations.mul.MatScalarMul(other, self)
+
+    def __matmul__(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, MatConst, other)
+        return operations.MatMul(self, other)
+
+    def __rmatmul(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, MatConst, other)
+        return operations.mul.MatMul(other, self)
+
+    def __truediv__(self, scalar):
+        scalar = self._wrap_if_constant(scalar)
+        raise NotImplementedError
+
+    def transpose(self) -> MatrixExpr:
+        """ Matrix transposition. """
+        raise NotImplementedError
+        return operations.transpose.MatTranspose(self)
+
+    @property
+    def T(self) -> MatrixExpr:
+        """ Shorthand for :py:meth:`mdpoly.expressions.matrix.MatrixExpr.transpose`. """
+        return self.transpose()
+
+    def to_scalar(self, scalar_type: type):
+        """ Convert to a scalar expression. """
+        raise NotImplementedError
+
+
+@dataclassabc(frozen=True)
+class MatConst(Const, MatrixExpr):
+    """ Matrix constant. TODO: desc. """
+    value: Sequence[Sequence[Number]] # Row major, overloads Const.value
+    shape: Shape # overloads Expr.shape
+    name: str = "" # overloads Leaf.name
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        r = repr_type(self.shape)
+
+        for i, row in enumerate(self.value):
+            for j, val in enumerate(row):
+                r.set(MatrixIndex(row=i, col=j), PolyIndex.constant(), val)
+
+        return r, state
+
+
+    def __str__(self) -> str:
+        if not self.name:
+            return repr(self.value)
+
+        return self.name
+
+
+
+T = TypeVar("T", bound=Var)
+
+@dataclassabc(frozen=True)
+class MatVar(Var, MatrixExpr):
+    """ Matrix of polynomial variables. TODO: desc """
+    name: str # overloads Leaf.name
+    shape: Shape # overloads Expr.shape
+
+    # TODO: review this API, can be moved elsewhere?
+    def to_scalars(self, scalar_var_type: Type[T]) -> Iterable[tuple[MatrixIndex, T]]:
+        for row in range(self.shape.rows):
+            for col in range(self.shape.cols):
+                var = scalar_var_type(name=f"{self.name}_[{row},{col}]") # type: ignore[call-arg]
+                entry = MatrixIndex(row, col)
+
+                yield entry, var
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        r = repr_type(self.shape)
+        
+        # FIXME: do not hardcode scalar type
+        for entry, var in self.to_scalars(PolyVar):
+            idx = PolyVarIndex.from_var(var, state), # important comma!
+            r.set(entry, PolyIndex(idx), 1)
+
+        return r, state
+
+    def __str__(self) -> str:
+        return self.name
+
+
+@dataclassabc(frozen=True)
+class MatParam(Param, MatrixExpr):
+    """ Matrix parameter. TODO: desc. """
+    name: str
+    shape: Shape
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        if self not in state.parameters:
+            raise MissingParameters("Cannot construct representation because "
+                                    f"value for parameter {self} was not given.")
+
+        # FIXME: add conversion to scalar variables
+        # Ignore typecheck because dataclassabc has not type stub
+        return MatConst(state.parameters[self]).to_repr(repr_type, state) # type: ignore[call-arg]
+
+    def __str__(self) -> str:
+        return self.name
diff --git a/mdpoly/expressions/poly.py b/mdpoly/expressions/poly.py
new file mode 100644
index 0000000..3fe83ca
--- /dev/null
+++ b/mdpoly/expressions/poly.py
@@ -0,0 +1,169 @@
+from __future__ import annotations
+from typing import TYPE_CHECKING
+
+from typing import Type, Iterable
+from functools import reduce
+from itertools import chain, combinations_with_replacement
+from dataclassabc import dataclassabc
+
+from ..abc import Expr, Var, Const, Param, Algebra, ReprT
+from ..index import Shape, MatrixIndex, PolyIndex, PolyVarIndex, Number
+from ..errors import AlgebraicError, MissingParameters, InvalidShape
+
+from .. import operations
+
+if TYPE_CHECKING:
+    from ..index import Number
+    from ..state import State
+
+
+class PolyRingExpr(Expr):
+    r""" Expression with the algebraic behaviour of a polynomial ring.
+
+    This is the algebra of :math:`\mathbb{R}[x_1, \ldots, x_n]`. Note that the
+    polynomials are scalars.
+    """
+
+    # -- Properties ---
+
+    @property
+    def algebra(self) -> Algebra:
+        return Algebra.poly_ring
+
+    @property
+    def shape(self) -> Shape:
+        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
+        if self.left.shape != self.right.shape:
+            raise InvalidShape(f"Cannot perform operation {repr(self)} with "
+                                f"shapes {self.left.shape} and {self.right.shape}.")
+        return self.left.shape
+
+    # --- Helper methods for construction ---
+
+    @classmethod
+    def from_powers(cls, variable: PolyRingExpr, exponents: Iterable[int]) -> PolyRingExpr:
+        """ Given a variable, say :math:`x`, and a list of exponents ``[0, 3, 5]``
+        returns a polynomial :math:`x^5 + x^3 + 1`.
+        """
+        nonzero_exponents = filter(lambda e: e != 0, exponents)
+        # Ignore typecheck because dataclassabc has no typing stub
+        constant_term = PolyConst(1 if 0 in exponents else 0) # type: ignore[call-arg]
+        # FIXME: remove annoying 0
+        return sum((variable ** e for e in nonzero_exponents), constant_term)
+
+
+    @classmethod
+    def make_combinations(cls, variables: Iterable[PolyVar], max_degree: int) -> Iterable[PolyRingExpr]:
+        """ Make combinations of terms.
+
+        For example given :math:`x, y` and *max_degree* of 2 generates
+        :math:`x^2, xy, x, y^2, y, 1`.
+        """
+        # Ignore typecheck because dataclassabc has no typing stub
+        vars_and_const = chain(variables, (PolyConst(1),)) # type: ignore[call-arg]
+        for comb in combinations_with_replacement(vars_and_const, max_degree):
+            yield cast(PolyRingExpr, reduce(operator.mul, comb))
+
+    # --- Operator Overloading ---
+
+    def __add__(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, PolyConst, other)
+        return operations.add.PolyAdd(self, other)
+
+    def __radd__(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, PolyConst, other)
+        return operations.add.PolyAdd(other, self)
+
+    def __sub__(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, PolyConst, other)
+        return operations.add.PolySub(self, other)
+
+    def __rsub__(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, PolyConst, other)
+        return operations.add.PolyAdd(other, self)
+
+    def __neg__(self):
+        # FIXME: Create PolyNeg?
+        return operations.mul.PolyMul(PolyConst(-1), self)
+
+    def __mul__(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, PolyConst, other)
+        return operations.mul.PolyMul(self, other)
+
+    def __rmul__(self, other):
+        self._assert_same_algebra(self, other)
+        other = self._wrap(Number, PolyConst, other)
+        return operations.mul.PolyMul(other, self)
+
+    def __matmul__(self, other):
+        raise AlgebraicError("Cannot perform matrix multiplication in polynomial ring (they are scalars).")
+
+    def __rmatmul__(self, other):
+        self.__rmatmul__(other)
+
+    def __truediv__(self, other):
+        other = self._wrap(Number, PolyConst, other)
+        if not isinstance(other, PolyConst | PolyParam):
+            raise AlgebraicError("Cannot divide by variables in polynomial ring.")
+
+        return operations.mul.PolyMul(self, other)
+
+    def __rtruediv__(self, other):
+        raise AlgebraicError("Cannot perform right division in polynomial ring.")
+
+    def __pow__(self, other):
+        other = self._wrap(Number, PolyConst, other)
+        if not isinstance(other, PolyConst | PolyParam):
+             raise AlgebraicError(f"Cannot raise to powers of type {type(other)} in "
+                                  "polynomial ring. Only constants and parameters are allowed.")
+        return operations.exp.PolyExp(left=self, right=other)
+
+    # -- Other mathematical operations ---
+    
+    def diff(self, wrt: PolyVar) -> operations.derivative.PolyPartialDiff:
+        return operations.derivative.PolyPartialDiff(right=self, wrt=wrt)
+
+
+@dataclassabc(frozen=True)
+class PolyVar(Var, PolyRingExpr):
+    """ Variable TODO: desc """
+    name: str # overloads Leaf.name
+    shape: Shape = Shape.scalar() # ovearloads PolyRingExpr.shape
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        r = repr_type(self.shape)
+        idx = PolyVarIndex.from_var(self, state), # important comma!
+        r.set(MatrixIndex.scalar(), PolyIndex(idx), 1)
+        return r, state
+
+
+@dataclassabc(frozen=True)
+class PolyConst(Const, PolyRingExpr):
+    """ Constant TODO: desc """
+    value: Number # overloads Const.value
+    name: str = "" # overloads Leaf.name
+    shape: Shape = Shape.scalar() # ovearloads PolyRingExpr.shape
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        r = repr_type(self.shape)
+        r.set(MatrixIndex.scalar(), PolyIndex.constant(), self.value)
+        return r, state
+
+
+@dataclassabc(frozen=True)
+class PolyParam(Param, PolyRingExpr):
+    """ Polynomial parameter TODO: desc """
+    name: str # overloads Leaf.name
+    shape: Shape = Shape.scalar() # overloads PolyRingExpr.shape
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        if self not in state.parameters:
+            raise MissingParameters("Cannot construct representation because "
+                                    f"value for parameter {self} was not given.")
+
+        return PolyConst(state.parameters[self]).to_repr(repr_type, state) # type: ignore[call-arg]
diff --git a/mdpoly/index.py b/mdpoly/index.py
index 5211245..070412c 100644
--- a/mdpoly/index.py
+++ b/mdpoly/index.py
@@ -1,15 +1,17 @@
 from __future__ import annotations
+from typing import TYPE_CHECKING
+
+from typing import Self, NamedTuple, Optional, Sequence
+from itertools import filterfalse
 
 from .errors import InvalidShape
 from .util import partition, isclose
 
-from itertools import filterfalse
-from typing import Self, NamedTuple, Optional, Sequence, TYPE_CHECKING
-
 if TYPE_CHECKING:
     from .state import State
     from .abc import Var
 
+# FIXME: this is not really an appropriate place for this -> constants?
 Number = int | float
 
 
@@ -206,6 +208,11 @@ class PolyIndex(tuple[PolyVarIndex]):
         return cls(PolyVarIndex(k, v) for k, v in d.items())
 
     @classmethod
+    def from_expr(cls, product_of_vars) -> Self:
+        # TODO: implement product of variables
+        raise NotImplementedError
+
+    @classmethod
     def constant(cls) -> Self:
         """ Index of the constant term. """
         return cls((PolyVarIndex.constant(),))
diff --git a/mdpoly/operations/__init__.py b/mdpoly/operations/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/mdpoly/operations/__init__.py
diff --git a/mdpoly/operations/add.py b/mdpoly/operations/add.py
new file mode 100644
index 0000000..bf0c5a0
--- /dev/null
+++ b/mdpoly/operations/add.py
@@ -0,0 +1,81 @@
+from __future__ import annotations
+from typing import TYPE_CHECKING
+
+from typing import Type
+from dataclassabc import dataclassabc
+
+from ..abc import Expr
+from ..errors import AlgebraicError
+
+from ..expressions import BinaryOp, Reducible
+from ..expressions.matrix import MatrixExpr
+from ..expressions.poly import PolyRingExpr
+
+if TYPE_CHECKING:
+    from ..abc import ReprT
+    from ..index import Shape
+    from ..state import State
+
+
+@dataclassabc(eq=False)
+class MatAdd(BinaryOp, MatrixExpr):
+    """ Addition operator between matrices. """
+
+    @property
+    def shape(self) -> Shape:
+        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
+        if not self.left.shape == self.right.shape:
+            raise AlgebraicError(f"Cannot add matrices {self.left} and {self.right} "
+                                 "with different shapes.")
+        return self.left.shape
+
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
+        # Make a new empty representation
+        r = repr_type(self.shape)
+
+        # Make representations for existing stuff
+        for node in self.children():
+            nrepr, state = node.to_repr(repr_type, state)
+
+            # Process non-zero terms
+            for entry in nrepr.entries():
+                for term in nrepr.terms(entry):
+                    s = r.at(entry, term) + nrepr.at(entry, term)
+                    r.set(entry, term, s)
+
+        return r, state
+
+    def __str__(self) -> str:
+        return f"({self.left} + {self.right})"
+
+
+@dataclassabc(eq=False)
+class MatSub(BinaryOp, MatrixExpr, Reducible):
+    """ Subtraction operator between matrices. """
+
+    @property
+    def shape(self) -> Shape:
+        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
+        if not self.left.shape == self.right.shape:
+            raise AlgebraicError(f"Cannot subtract matrices {self.left} and {self.right} "
+                                 f"with different shapes, {self.left.shape} and {self.right.shape}.")
+        return self.left.shape
+
+    def reduce(self) -> Expr:
+        """ See :py:meth:`mdpoly.expressions.Reducible.reduce`. """
+        return self.left + (-1 * self.right)
+
+    def __str__(self) -> str:
+        return f"({self.left} - {self.right})"
+
+
+@dataclassabc(eq=False)
+class PolyAdd(MatAdd, PolyRingExpr):
+    ...
+
+
+@dataclassabc(eq=False)
+class PolySub(MatSub, PolyRingExpr):
+    ...
diff --git a/mdpoly/operations/derivative.py b/mdpoly/operations/derivative.py
new file mode 100644
index 0000000..6e02b0e
--- /dev/null
+++ b/mdpoly/operations/derivative.py
@@ -0,0 +1,62 @@
+from __future__ import annotations
+from typing import TYPE_CHECKING
+
+from typing import Type, Self, cast
+from dataclassabc import dataclassabc
+
+from ..abc import Expr
+from ..errors import AlgebraicError
+from ..index import Shape, MatrixIndex, PolyIndex
+
+from ..expressions import UnaryOp
+from ..expressions.poly import PolyRingExpr, PolyVar
+
+if TYPE_CHECKING:
+    from ..abc import ReprT
+    from ..state import State
+
+
+# TODO: implement vector derivatives
+
+
+@dataclassabc(eq=False)
+class PolyPartialDiff(UnaryOp, PolyRingExpr):
+    """ Partial differentiation of scalar polynomials. """
+    wrt: PolyVar
+
+    @property
+    def shape(self) -> Shape:
+        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
+        return self.inner.shape
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
+        r = repr_type(self.shape)
+        lrepr, state = self.inner.to_repr(repr_type, state)
+
+        entry = MatrixIndex.scalar()
+        wrt = state.index(self.wrt)
+
+        for term in lrepr.terms(entry):
+            # sadly, walrus operator does not support unpacking
+            # https://bugs.python.org/issue43143
+            # if ((newterm, fac) := PolyIndex.differentiate(term, wrt)):
+            if result := PolyIndex.differentiate(term, wrt):
+                newterm, fac = result
+                r.set(entry, newterm, fac * lrepr.at(entry, term))
+
+        return r, state
+
+    def __str__(self) -> str:
+        return f"(∂_{self.wrt} {self.inner})"
+
+    def replace(self, old: PolyRingExpr, new: PolyRingExpr) -> Self:
+        """ Overloads :py:meth:`mdpoly.abc.Expr.replace` """
+        if self.wrt == old:
+            if not isinstance(new, PolyVar):
+                # FIXME: implement chain rule
+                raise AlgebraicError(f"Cannot take a derivative with respect to {new}.")
+
+            self.wrt = cast(PolyVar, new)
+
+        return cast(Self, Expr.replace(self, old, new))
diff --git a/mdpoly/operations/exp.py b/mdpoly/operations/exp.py
new file mode 100644
index 0000000..fbc67f5
--- /dev/null
+++ b/mdpoly/operations/exp.py
@@ -0,0 +1,52 @@
+from __future__ import annotations
+from typing import TYPE_CHECKING
+
+from operator import mul as opmul
+from typing import cast
+from functools import reduce
+from dataclasses import dataclass
+
+from ..abc import Expr, Const
+from ..errors import AlgebraicError
+
+from ..expressions import BinaryOp, Reducible
+from ..expressions.matrix import MatrixExpr
+from ..expression.poly import PolyRingExpr
+
+
+# TODO: implement matrix exponential, use caley-hamilton thm magic
+@dataclass(eq=False)
+class MatExp(BinaryOp, MatrixExpr, Reducible):
+    def __init__(self):
+        raise NotImplementedError
+
+
+@dataclass(eq=False)
+class PolyExp(BinaryOp, PolyRingExpr, Reducible):
+    """ Exponentiation operator between scalar polynomials. """
+
+    @property # type: ignore[override]
+    def right(self) -> Const: # type: ignore[override]
+        if not isinstance(super().right, Const):
+            raise AlgebraicError(f"Cannot raise {self.left} to {self.right} because"
+                                 f"{self.right} is not a constant.")
+
+        return cast(Const, super().right)
+
+    @right.setter
+    def right(self, right: Expr) -> None:
+        # Workaround necessary because of a bug
+        # https://bugs.python.org/issue14965
+        BinaryOp.right.fset(self, right) # type: ignore
+
+    def reduce(self) -> Expr:
+        """ See :py:meth:`mdpoly.expressions.Reducible.reduce`. """
+        var = self.left
+        if not isinstance(self.right.value, int):
+            raise NotImplementedError("Cannot raise to non-integer powers (yet).")
+
+        ntimes = self.right.value - 1
+        return reduce(opmul, (var for _ in range(ntimes)), var)
+
+    def __str__(self) -> str:
+        return f"({self.left} ** {self.right})"
diff --git a/mdpoly/operations/mul.py b/mdpoly/operations/mul.py
new file mode 100644
index 0000000..a2a3cc4
--- /dev/null
+++ b/mdpoly/operations/mul.py
@@ -0,0 +1,186 @@
+from __future__ import annotations
+from typing import TYPE_CHECKING
+
+from typing import Type
+from itertools import product
+from dataclasses import dataclass
+from dataclassabc import dataclassabc
+
+from ..index import Shape
+from ..errors import AlgebraicError, InvalidShape
+from ..index import MatrixIndex, PolyIndex
+
+from ..expressions import BinaryOp, Reducible
+from ..expressions.matrix import MatrixExpr
+from ..expression.poly import PolyRingExpr
+
+if TYPE_CHECKING:
+    from ..abc import ReprT
+    from ..state import State
+
+
+# ┏┳┓┏━┓╺┳╸┏━┓╻╻ ╻   ┏━┓┏━┓┏━┓╺┳┓╻ ╻┏━╸╺┳╸┏━┓
+# ┃┃┃┣━┫ ┃ ┣┳┛┃┏╋┛   ┣━┛┣┳┛┃ ┃ ┃┃┃ ┃┃   ┃ ┗━┓
+# ╹ ╹╹ ╹ ╹ ╹┗╸╹╹ ╹   ╹  ╹┗╸┗━┛╺┻┛┗━┛┗━╸ ╹ ┗━┛
+
+
+@dataclassabc
+class MatElemMul(BinaryOp, MatrixExpr):
+    """ Elementwise Matrix Multiplication. """
+
+    @property
+    def shape(self) -> Shape:
+        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
+        if not self.left.shape == self.right.shape:
+            raise AlgebraicError("Cannot perform element-wise multiplication of matrices "
+                                 f"{self.left} and {self.right} with different shapes, "
+                                 f"{self.left.shape} and {self.right.shape}")
+        return self.left.shape
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
+        r = repr_type(self.shape)
+
+        lrepr, state = self.left.to_repr(repr_type, state)
+        rrepr, state = self.right.to_repr(repr_type, state)
+
+        # Non zero entries are the intersection since if either is zero the
+        # result is zero
+        nonzero_entries = set(lrepr.entries()) & set(rrepr.entries())
+        for entry in nonzero_entries:
+            # Compute polynomial product between non-zero entries
+            for lterm, rterm in product(lrepr.terms(entry), rrepr.terms(entry)):
+                # Compute where the results should go
+                term = PolyIndex.product(lterm, rterm)
+
+                # Compute product
+                p = r.at(entry, term) + lrepr.at(entry, lterm) * rrepr.at(entry, rterm)
+                r.set(entry, term, p)
+
+        return r, state
+
+    def __str__(self) -> str:
+        return f"({self.left} .* {self.right})"
+
+
+@dataclassabc
+class MatScalarMul(BinaryOp, MatrixExpr):
+    """ Matrix-Scalar Multiplication. Assumes scalar is on the left and matrix
+    on the right. """
+
+    @property
+    def shape(self) -> Shape:
+        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
+        if not self.left.shape == Shape.scalar():
+            raise InvalidShape(f"Matrix-scalar product assumes that left argumet {self.left} "
+                               f"but it has shape {self.left.shape}")
+
+        return self.right.shape
+
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
+        r = repr_type(self.shape)
+
+        scalar_repr, state = self.left.to_repr(repr_type, state)
+        mat_repr, state = self.right.to_repr(repr_type, state)
+
+        for entry in mat_repr.entries():
+            scalar_terms = scalar_repr.terms(MatrixIndex.scalar())
+            mat_terms = mat_repr.terms(entry)
+            for scalar_term, mat_term in product(scalar_terms, mat_terms):
+                term = PolyIndex.product(scalar_term, mat_term)
+
+                p = r.at(entry, term) + scalar_repr.at(entry, scalar_term) + mat_repr.at(entry, mat_term)
+                r.set(entry, term, p)
+
+        return r, state
+
+    def __str__(self) -> str:
+        return f"({self.left} * {self.right})"
+
+
+@dataclass(eq=False)
+class MatMul(BinaryOp, MatrixExpr):
+    """ Matrix Multiplication. """
+
+    @property
+    def shape(self) -> Shape:
+        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
+        if not self.left.shape.rows == self.right.shape.cols:
+            raise AlgebraicError("Cannot perform matrix multiplication between "
+                                 f"{self.left} and {self.right} (shapes {self.left.shape} and {self.right.shape})")
+
+        return Shape(self.left.shape.rows, self.right.shape.cols)
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
+        r = repr_type(self.shape)
+
+        lrepr, state = self.left.to_repr(repr_type, state)
+        rrepr, state = self.right.to_repr(repr_type, state)
+
+        # Compute matrix product
+        for row in range(self.left.shape.rows):
+            for col in range(self.right.shape.cols):
+                for i in range(self.left.shape.cols):
+                    raise NotImplementedError
+
+        return r, state
+
+
+    def __str__(self) -> str:
+        return f"({self.left} @ {self.right})"
+
+
+@dataclass(eq=False)
+class MatDotProd(BinaryOp, MatrixExpr, Reducible):
+    """ Dot product. """
+
+    @property
+    def shape(self) -> Shape:
+        if not self.left.shape.is_row():
+            raise AlgebraicError(f"Left operand {self.left} must be a row!")
+
+        if not self.right.shape.is_col():
+            raise AlgebraicError(f"Right operand {self.right} must be a column!")
+
+        if self.left.shape.cols != self.right.shape.rows:
+            raise AlgebraicError(f"Rows of {self.right} and columns {self.left} do not match!")
+
+        return Shape.scalar()
+
+
+# ┏━┓┏━┓╻  ╻ ╻┏┓╻┏━┓┏┳┓╻┏━┓╻     ┏━┓┏━┓┏━┓╺┳┓╻ ╻┏━╸╺┳╸
+# ┣━┛┃ ┃┃  ┗┳┛┃┗┫┃ ┃┃┃┃┃┣━┫┃     ┣━┛┣┳┛┃ ┃ ┃┃┃ ┃┃   ┃ 
+# ╹  ┗━┛┗━╸ ╹ ╹ ╹┗━┛╹ ╹╹╹ ╹┗━╸   ╹  ╹┗╸┗━┛╺┻┛┗━┛┗━╸ ╹ 
+
+
+@dataclass(eq=False)
+class PolyMul(BinaryOp, PolyRingExpr):
+    """ Multiplication operator between scalar polynomials. """
+
+    def to_repr(self, repr_type: Type[ReprT], state: State) -> tuple[ReprT, State]:
+        """ See :py:meth:`mdpoly.abc.Expr.to_repr`. """
+        r = repr_type(self.shape)
+
+        lrepr, state = self.left.to_repr(repr_type, state)
+        rrepr, state = self.right.to_repr(repr_type, state)
+
+        # Non zero entries are the intersection since if either is zero the
+        # result is zero
+        nonzero_entries = set(lrepr.entries()) & set(rrepr.entries())
+        for entry in nonzero_entries:
+            # Compute polynomial product between non-zero entries
+            for lterm, rterm in product(lrepr.terms(entry), rrepr.terms(entry)):
+                # Compute where the results should go
+                term = PolyIndex.product(lterm, rterm)
+
+                # Compute product
+                p = r.at(entry, term) + lrepr.at(entry, lterm) * rrepr.at(entry, rterm)
+                r.set(entry, term, p)
+
+        return r, state
+
+    def __str__(self) -> str:
+        return f"({self.left} * {self.right})"
diff --git a/mdpoly/operations/transpose.py b/mdpoly/operations/transpose.py
new file mode 100644
index 0000000..750d0c0
--- /dev/null
+++ b/mdpoly/operations/transpose.py
@@ -0,0 +1,20 @@
+from __future__ import annotations
+from typing import TYPE_CHECKING
+
+from dataclasses import dataclass
+
+from ..expressions import UnaryOp
+from ..expressions.matrix import MatrixExpr
+
+if TYPE_CHECKING:
+    from ..index import Shape
+
+
+@dataclass(eq=False)
+class MatTranspose(UnaryOp, MatrixExpr):
+    """ Matrix transposition """
+
+    @property
+    def shape(self) -> Shape:
+        """ See :py:meth:`mdpoly.abc.Expr.shape`. """
+        return Shape(self.inner.shape.cols, self.inner.shape.rows)
diff --git a/mdpoly/representations.py b/mdpoly/representations.py
index 8268c42..32c75d7 100644
--- a/mdpoly/representations.py
+++ b/mdpoly/representations.py
@@ -1,10 +1,14 @@
-from .abc import Repr 
-from .index import Number, Shape, MatrixIndex, PolyIndex
+from __future__ import annotations
+from typing import TYPE_CHECKING
 
 from typing import Iterable
-
 import numpy as np
-import numpy.typing as npt
+
+if TYPE_CHECKING:
+    import numpy.typing as npt
+
+    from .abc import Repr 
+    from .index import Number, Shape, MatrixIndex, PolyIndex
 
 
 class SparseRepr(Repr):
diff --git a/mdpoly/state.py b/mdpoly/state.py
index 0095916..092de84 100644
--- a/mdpoly/state.py
+++ b/mdpoly/state.py
@@ -1,5 +1,6 @@
 from __future__ import annotations
 from typing import TYPE_CHECKING
+
 from dataclasses import dataclass, field
 
 from .index import PolyVarIndex
diff --git a/mdpoly/util.py b/mdpoly/util.py
index ae726ef..52cbb87 100644
--- a/mdpoly/util.py
+++ b/mdpoly/util.py
@@ -1,7 +1,12 @@
-from .constants import NUMERICS_EPS
-
 from itertools import tee, filterfalse
+from typing import TYPE_CHECKING
+
 from typing import Iterable
+from .constants import NUMERICS_EPS
+
+if TYPE_CHECKING:
+    pass
+
 
 def partition(pred, iterable):
     """Partition entries into false entries and true entries.