diff --git a/src/Math/Matrix3.h b/src/Math/Matrix3.h
index a71854a7c..0e55d5169 100644
--- a/src/Math/Matrix3.h
+++ b/src/Math/Matrix3.h
@@ -86,8 +86,8 @@ template<class T> class Matrix3: public Matrix<3, T> {
         /**
          * @brief 2D reflection matrix
          * @param normal    Normal of the line through which to reflect
-         *      (normalized)
          *
+         * Expects that the normal is normalized.
          * @see Matrix4::reflection()
          */
         static Matrix3<T> reflection(const Vector2<T>& normal) {
diff --git a/src/Math/Matrix4.h b/src/Math/Matrix4.h
index d7f28f1a4..ce34891e0 100644
--- a/src/Math/Matrix4.h
+++ b/src/Math/Matrix4.h
@@ -72,12 +72,10 @@ template<class T> class Matrix4: public Matrix<4, T> {
          * @param angle             Rotation angle (counterclockwise, in radians)
          * @param normalizedAxis    Normalized rotation axis
          *
-         * If possible, use faster alternatives like rotationX(), rotationY()
-         * and rotationZ().
+         * Expects that the rotation axis is normalized. If possible, use
+         * faster alternatives like rotationX(), rotationY() and rotationZ().
          * @see rotation() const, Matrix3::rotation(T), Vector3::xAxis(),
          *      Vector3::yAxis(), Vector3::zAxis(), deg(), rad()
-         * @attention Assertion fails on non-normalized rotation vector and
-         *      identity matrix is returned.
          */
         static Matrix4<T> rotation(T angle, const Vector3<T>& normalizedAxis) {
             CORRADE_ASSERT(MathTypeTraits<T>::equals(normalizedAxis.dot(), T(1)),
@@ -174,8 +172,8 @@ template<class T> class Matrix4: public Matrix<4, T> {
         /**
          * @brief 3D reflection matrix
          * @param normal    Normal of the plane through which to reflect
-         *      (normalized)
          *
+         * Expects that the normal is normalized.
          * @see Matrix3::reflection()
          */
         static Matrix4<T> reflection(const Vector3<T>& normal) {
@@ -293,7 +291,7 @@ template<class T> class Matrix4: public Matrix<4, T> {
         /**
          * @brief Inverted Euclidean transformation matrix
          *
-         * Assumes that the matrix represents Euclidean transformation (i.e.
+         * Expects that the matrix represents Euclidean transformation (i.e.
          * only rotation and translation, no scaling) and creates inverted
          * matrix from transposed rotation part and negated translation part.
          * Significantly faster than the general algorithm in inverted().
diff --git a/src/Math/Quaternion.h b/src/Math/Quaternion.h
index da9ccf642..79e95d496 100644
--- a/src/Math/Quaternion.h
+++ b/src/Math/Quaternion.h
@@ -38,7 +38,9 @@ template<class T> class Quaternion {
          * @param angle             Rotation angle (counterclockwise, in radians)
          * @param normalizedAxis    Normalized rotation axis
          *
-         * Assumes that the rotation axis is normalized.
+         * Expects that the rotation axis is normalized. @f[
+         *      q = [\boldsymbol a \cdot sin \frac \theta 2, cos \frac \theta 2]
+         * @f]
          */
         inline static Quaternion<T> fromRotation(T angle, const Vector3<T>& normalizedAxis) {
             CORRADE_ASSERT(MathTypeTraits<T>::equals(normalizedAxis.dot(), T(1)),
@@ -72,7 +74,10 @@ template<class T> class Quaternion {
         /**
          * @brief Rotation angle of unit quaternion
          *
-         * Assumes that the quaternion is normalized.
+         * Expects that the quaternion is normalized. @f[
+         *      \theta = 2 \cdot acos q_S
+         * @f]
+         * @see rotationAxis(), fromRotation()
          */
         inline T rotationAngle() const {
             CORRADE_ASSERT(MathTypeTraits<T>::equals(lengthSquared(), T(1)),
@@ -84,7 +89,10 @@ template<class T> class Quaternion {
         /**
          * @brief Rotation axis of unit quaternion
          *
-         * Assumes that the quaternion is normalized.
+         * Expects that the quaternion is normalized. @f[
+         *      \boldsymbol a = \frac{\boldsymbol q_V}{\sqrt{1 - q_S^2}}
+         * @f]
+         * @see rotationAngle(), fromRotation()
          */
         inline Vector3<T> rotationAxis() const {
             CORRADE_ASSERT(MathTypeTraits<T>::equals(lengthSquared(), T(1)),
@@ -126,7 +134,9 @@ template<class T> class Quaternion {
         /**
          * @brief Multiply with scalar and assign
          *
-         * The computation is done in-place.
+         * The computation is done in-place. @f[
+         *      q \cdot a = [\boldsymbol q_V \cdot a, q_S \cdot a]
+         * @f]
          */
         inline Quaternion<T>& operator*=(T number) {
             _vector *= number;
@@ -135,23 +145,25 @@ template<class T> class Quaternion {
         }
 
         /**
-         * @brief Divide with scalar and assign
+         * @brief Multiply with scalar
          *
-         * The computation is done in-place.
+         * @see operator*=(T)
          */
-        inline Quaternion<T>& operator/=(T number) {
-            _vector /= number;
-            _scalar /= number;
-            return *this;
+        inline Quaternion<T> operator*(T scalar) const {
+            return Quaternion<T>(*this)*=scalar;
         }
 
         /**
-         * @brief Multiply with scalar
+         * @brief Divide with scalar and assign
          *
-         * @see operator*=(T)
+         * The computation is done in-place. @f[
+         *      \frac q a = [\frac {\boldsymbol q_V} a, \frac {q_S} a]
+         * @f]
          */
-        inline Quaternion<T> operator*(T scalar) const {
-            return Quaternion<T>(*this)*=scalar;
+        inline Quaternion<T>& operator/=(T number) {
+            _vector /= number;
+            _scalar /= number;
+            return *this;
         }
 
         /**
@@ -166,7 +178,10 @@ template<class T> class Quaternion {
         /**
          * @brief Multiply with quaternion
          *
-         * The computation is *not* done in-place.
+         * @f[
+         *      p q = [p_S \boldsymbol q_V + q_S \boldsymbol p_V + \boldsymbol p_V \times \boldsymbol q_V,
+         *             p_S q_S - \boldsymbol p_V \cdot \boldsymbol q_V]
+         * @f]
          */
         inline Quaternion<T> operator*(const Quaternion<T>& other) const {
             return {_scalar*other._vector + other._scalar*_vector + Vector3<T>::cross(_vector, other._vector),
@@ -174,30 +189,38 @@ template<class T> class Quaternion {
         }
 
         /**
-         * @brief Multiply with quaternion and assign
+         * @brief %Quaternion length squared
          *
-         * @see operator*(const Quaternion<T>&) const
+         * Should be used instead of length() for comparing quaternion length
+         * with other values, because it doesn't compute the square root. @f[
+         *      |q|^2 = \boldsymbol q_V \cdot \boldsymbol q_V + q_S q_S
+         * @f]
          */
-        inline Quaternion<T>& operator*=(const Quaternion<T>& other) {
-            return (*this = *this * other);
-        }
-
-        /** @brief %Quaternion length squared */
         inline T lengthSquared() const {
             return _vector.dot() + _scalar*_scalar;
         }
 
-        /** @brief %Quaternion length */
+        /**
+         * @brief %Quaternion length
+         *
+         * @see lengthSquared()
+         */
         inline T length() const {
             return std::sqrt(lengthSquared());
         }
 
-        /** @brief Normalized quaternion */
+        /** @brief Normalized quaternion (of length 1) */
         inline Quaternion<T> normalized() const {
             return (*this)/length();
         }
 
-        /** @brief Conjugated quaternion */
+        /**
+         * @brief Conjugated quaternion
+         *
+         * @f[
+         *      q^* = [-\boldsymbol q_V, q_S]
+         * @f]
+         */
         inline Quaternion<T> conjugated() const {
             return {-_vector, _scalar};
         }
@@ -206,7 +229,9 @@ template<class T> class Quaternion {
          * @brief Inverted quaternion
          *
          * See invertedNormalized() which is faster for normalized
-         * quaternions.
+         * quaternions. @f[
+         *      q^{-1} = \frac{q^*}{|q|^2} = \frac{[-\boldsymbol q_V, q_S]}{|q|^2}
+         * @f]
          */
         inline Quaternion<T> inverted() const {
             return conjugated()/lengthSquared();
@@ -215,8 +240,10 @@ template<class T> class Quaternion {
         /**
          * @brief Inverted normalized quaternion
          *
-         * Equivalent to conjugated(). Assumes that the quaternion is
-         * normalized.
+         * Equivalent to conjugated(). Expects that the quaternion is
+         * normalized. @f[
+         *      q^{-1} = q^* = [-\boldsymbol q_V, q_S] ~~~~~ |q| = 1
+         * @f]
          */
         inline Quaternion<T> invertedNormalized() const {
             CORRADE_ASSERT(MathTypeTraits<T>::equals(lengthSquared(), T(1)),
diff --git a/src/Math/RectangularMatrix.h b/src/Math/RectangularMatrix.h
index e87e8c4bb..691174766 100644
--- a/src/Math/RectangularMatrix.h
+++ b/src/Math/RectangularMatrix.h
@@ -63,6 +63,10 @@ template<std::size_t size, class T> class Vector;
 
 See @ref matrix-vector for brief introduction. See also Matrix (square) and
 Vector.
+
+The data are stored in column-major order, to reflect that, all indices in
+math formulas are in reverse order (i.e. @f$ \boldsymbol A_{ji} @f$ instead
+of @f$ \boldsymbol A_{ij} @f$).
 */
 template<std::size_t cols, std::size_t rows, class T> class RectangularMatrix {
     static_assert(cols != 0 && rows != 0, "Matrix cannot have zero elements");
@@ -199,6 +203,20 @@ template<std::size_t cols, std::size_t rows, class T> class RectangularMatrix {
             return !operator==(other);
         }
 
+        /**
+         * @brief Add and assign matrix
+         *
+         * The computation is done in-place. @f[
+         *      \boldsymbol A_{ji} = \boldsymbol A_{ji} + \boldsymbol B_{ji}
+         * @f]
+         */
+        RectangularMatrix<cols, rows, T>& operator+=(const RectangularMatrix<cols, rows, T>& other) {
+            for(std::size_t i = 0; i != cols*rows; ++i)
+                _data[i] += other._data[i];
+
+            return *this;
+        }
+
         /**
          * @brief Add matrix
          *
@@ -209,19 +227,12 @@ template<std::size_t cols, std::size_t rows, class T> class RectangularMatrix {
         }
 
         /**
-         * @brief Add and assign matrix
+         * @brief Negated matrix
          *
-         * More efficient than operator+(), because it does the computation
-         * in-place.
+         * The computation is done in-place. @f[
+         *      \boldsymbol A_{ji} = -\boldsymbol A_{ji}
+         * @f]
          */
-        RectangularMatrix<cols, rows, T>& operator+=(const RectangularMatrix<cols, rows, T>& other) {
-            for(std::size_t i = 0; i != cols*rows; ++i)
-                _data[i] += other._data[i];
-
-            return *this;
-        }
-
-        /** @brief Negated matrix */
         RectangularMatrix<cols, rows, T> operator-() const {
             RectangularMatrix<cols, rows, T> out;
 
@@ -231,20 +242,12 @@ template<std::size_t cols, std::size_t rows, class T> class RectangularMatrix {
             return out;
         }
 
-        /**
-         * @brief Subtract matrix
-         *
-         * @see operator-=()
-         */
-        inline RectangularMatrix<cols, rows, T> operator-(const RectangularMatrix<cols, rows, T>& other) const {
-            return RectangularMatrix<cols, rows, T>(*this)-=other;
-        }
-
         /**
          * @brief Subtract and assign matrix
          *
-         * More efficient than operator-(), because it does the computation
-         * in-place.
+         * The computation is done in-place. @f[
+         *      \boldsymbol A_{ji} = \boldsymbol A_{ji} - \boldsymbol B_{ji}
+         * @f]
          */
         RectangularMatrix<cols, rows, T>& operator-=(const RectangularMatrix<cols, rows, T>& other) {
             for(std::size_t i = 0; i != cols*rows; ++i)
@@ -253,6 +256,15 @@ template<std::size_t cols, std::size_t rows, class T> class RectangularMatrix {
             return *this;
         }
 
+        /**
+         * @brief Subtract matrix
+         *
+         * @see operator-=()
+         */
+        inline RectangularMatrix<cols, rows, T> operator-(const RectangularMatrix<cols, rows, T>& other) const {
+            return RectangularMatrix<cols, rows, T>(*this)-=other;
+        }
+
         /**
          * @brief Multiply matrix with number
          *
@@ -269,8 +281,9 @@ template<std::size_t cols, std::size_t rows, class T> class RectangularMatrix {
         /**
          * @brief Multiply matrix with number and assign
          *
-         * More efficient than operator*(U), because it does the computation
-         * in-place.
+         * The computation is done in-place. @f[
+         *      \boldsymbol A_{ji} = a \boldsymbol A_{ji}
+         * @f]
          */
         #ifndef DOXYGEN_GENERATING_OUTPUT
         template<class U> inline typename std::enable_if<std::is_arithmetic<U>::value, RectangularMatrix<cols, rows, T>&>::type operator*=(U number) {
@@ -299,8 +312,9 @@ template<std::size_t cols, std::size_t rows, class T> class RectangularMatrix {
         /**
          * @brief Divide matrix with number and assign
          *
-         * More efficient than operator/(), because it does the computation
-         * in-place.
+         * The computation is done in-place. @f[
+         *      \boldsymbol A_{ji} = \frac{\boldsymbol A_{ji}} a
+         * @f]
          */
         #ifndef DOXYGEN_GENERATING_OUTPUT
         template<class U> inline typename std::enable_if<std::is_arithmetic<U>::value, RectangularMatrix<cols, rows, T>&>::type operator/=(U number) {
@@ -313,7 +327,13 @@ template<std::size_t cols, std::size_t rows, class T> class RectangularMatrix {
             return *this;
         }
 
-        /** @brief Multiply matrix */
+        /**
+         * @brief Multiply matrix
+         *
+         * @f[
+         *      (\boldsymbol {AB})_{ji} = \sum_{k=0}^{m-1} \boldsymbol A_{ki} \boldsymbol B_{jk}
+         * @f]
+         */
         template<std::size_t size> RectangularMatrix<size, rows, T> operator*(const RectangularMatrix<size, cols, T>& other) const {
             RectangularMatrix<size, rows, T> out;
 
@@ -329,7 +349,9 @@ template<std::size_t cols, std::size_t rows, class T> class RectangularMatrix {
          * @brief Multiply vector
          *
          * Internally the same as multiplying with one-column matrix, but
-         * returns vector.
+         * returns vector. @f[
+         *      (\boldsymbol {Aa})_i = \sum_{k=0}^{m-1} \boldsymbol A_{ki} \boldsymbol a_k
+         * @f]
          */
         Vector<rows, T> operator*(const Vector<rows, T>& other) const {
             return operator*(static_cast<RectangularMatrix<1, rows, T>>(other));
diff --git a/src/Math/Vector.h b/src/Math/Vector.h
index c8ad4ec9d..838b75c25 100644
--- a/src/Math/Vector.h
+++ b/src/Math/Vector.h
@@ -57,10 +57,9 @@ template<std::size_t size, class T> class Vector: public RectangularMatrix<1, si
         /**
          * @brief Angle between normalized vectors (in radians)
          *
-         * @f[
-         * \phi = acos \left(\frac{a \cdot b}{|a| \cdot |b|} \right)
+         * Expects that both vectors are normalized. @f[
+         * \theta = acos \left(\frac{a \cdot b}{|a| \cdot |b|} \right)
          * @f]
-         * @attention Both vectors must be normalized.
          */
         inline static T angle(const Vector<size, T>& normalizedA, const Vector<size, T>& normalizedB) {
             CORRADE_ASSERT(MathTypeTraits<T>::equals(normalizedA.dot(), T(1)) && MathTypeTraits<T>::equals(normalizedB.dot(), T(1)),
@@ -200,10 +199,9 @@ template<std::size_t size, class T> class Vector: public RectangularMatrix<1, si
          * @brief Dot product of the vector
          *
          * Should be used instead of length() for comparing vector length with
-         * other values, because it doesn't compute the square root, just the
-         * dot product: @f$ a \cdot a < length \cdot length @f$ is faster
-         * than @f$ \sqrt{a \cdot a} < length @f$.
-         *
+         * other values, because it doesn't compute the square root. @f[
+         *      |a|^2 = a \cdot a
+         * @f]
          * @see dot(const Vector<size, T>&, const Vector<size, T>&)
          */
         inline T dot() const {
@@ -225,9 +223,11 @@ template<std::size_t size, class T> class Vector: public RectangularMatrix<1, si
         }
 
         /**
-         * @brief Vector projected onto another
+         * @brief %Vector projected onto another
          *
-         * Returns vector projected onto line defined by @p other.
+         * Returns vector projected onto line defined by @p other. @f[
+         *      \boldsymbol a_1 = \frac{\boldsymbol a \cdot \boldsymbol b}{\boldsymbol b \cdot \boldsymbol b} \boldsymbol b
+         * @f]
          */
         inline Vector<size, T> projected(const Vector<size, T>& other) const {
             return other*dot(*this, other)/other.dot();
diff --git a/src/Math/Vector3.h b/src/Math/Vector3.h
index 1e81d7f8e..ddec15a4a 100644
--- a/src/Math/Vector3.h
+++ b/src/Math/Vector3.h
@@ -94,6 +94,7 @@ template<class T> class Vector3: public Vector<3, T> {
          * @brief Cross product
          *
          * @f[
+         * \boldsymbol a \times \boldsymbol b =
          * \begin{pmatrix} c_0 \\ c_1 \\ c_2 \end{pmatrix} =
          * \begin{pmatrix}a_1b_2 - a_2b_1 \\ a_2b_0 - a_0b_2 \\ a_0b_1 - a_1b_0 \end{pmatrix}
          * @f]