From 1d3b7c95cbe111fe0191cea62374fb3630f2e202 Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Sun, 2 Jun 2019 14:40:29 +0300
Subject: [PATCH 01/16] Fix a bug in the prefix version of vector compose
 function

---
 source/linal/vector.h | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/source/linal/vector.h b/source/linal/vector.h
index c9b2a17..767cc88 100644
--- a/source/linal/vector.h
+++ b/source/linal/vector.h
@@ -174,9 +174,9 @@ template<typename T, unsigned N>
 inline Vector<T, N+1> compose(T s, const Vector<T, N> &v)
 {
 	Vector<T, N+1> r;
+	r[0] = s;
 	for(unsigned i=0; i<N; ++i)
-		r[i] = v[i];
-	r[N] = s;
+		r[i+1] = v[i];
 	return r;
 }
 
-- 
2.43.0


From 0a8781509bafe347b9c32f1106891eac18354965 Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Sun, 2 Jun 2019 15:28:34 +0300
Subject: [PATCH 02/16] Improve Vector constructor for C++11

The new overload allows vectors of any size to be constructed by passing
the component values as constructor arguments.
---
 source/linal/vector.h | 18 ++++++++++++++++++
 tests/vector.cpp      | 25 +++++++++++++++++++++++++
 2 files changed, 43 insertions(+)

diff --git a/source/linal/vector.h b/source/linal/vector.h
index 767cc88..31597d4 100644
--- a/source/linal/vector.h
+++ b/source/linal/vector.h
@@ -84,9 +84,14 @@ public:
 	use by Matrix row accessor. */
 	Vector(const T *, unsigned);
 
+#if __cplusplus >= 201103L
+	template<typename... Args>
+	Vector(T, Args...);
+#else
 	Vector(T, T);
 	Vector(T, T, T);
 	Vector(T, T, T, T);
+#endif
 	template<typename U>
 	Vector(const Vector<U, N> &);
 
@@ -126,6 +131,18 @@ inline Vector<T, N>::Vector(const T *d, unsigned stride)
 		(*this)[i] = d[i*stride];
 }
 
+#if __cplusplus >= 201103L
+template<typename T, unsigned N>
+template<typename... Args>
+inline Vector<T, N>::Vector(T x_, Args... v)
+{
+	static_assert(1+sizeof...(v)==N, "Incorrect number of arguments in Vector constructor");
+	(*this)[0] = x_;
+	unsigned i = 1;
+	for(auto c: std::initializer_list<T> { static_cast<T>(v)... })
+		(*this)[i++] = c;
+}
+#else
 /* The compiler won't instantiate these unless they are used.  Trying to use
 them on the wrong class results in an error. */
 template<typename T, unsigned N>
@@ -151,6 +168,7 @@ inline Vector<T, N>::Vector(T x_, T y_, T z_, T w_)
 	this->VectorComponents<T, 4>::z = z_;
 	this->VectorComponents<T, 4>::w = w_;
 }
+#endif
 
 template<typename T, unsigned N>
 template<typename U>
diff --git a/tests/vector.cpp b/tests/vector.cpp
index 14aea87..cf07329 100644
--- a/tests/vector.cpp
+++ b/tests/vector.cpp
@@ -18,6 +18,7 @@ public:
 	static const char *get_name() { return "Vector"; }
 
 private:
+	void constructors();
 	void component_aliases();
 	void composition();
 	void slice();
@@ -30,6 +31,7 @@ private:
 
 VectorTests::VectorTests()
 {
+	add(&VectorTests::constructors, "Constructors");
 	add(&VectorTests::component_aliases, "Component aliases");
 	add(&VectorTests::composition, "Compose");
 	add(&VectorTests::slice, "Slice");
@@ -40,6 +42,29 @@ VectorTests::VectorTests()
 	add(&VectorTests::unit_vec, "Normalize");
 }
 
+void VectorTests::constructors()
+{
+	static int data[] = { 1, 2, 3 };
+
+	Vector2i v1(1, 2);
+	EXPECT_EQUAL(v1[0], 1);
+	EXPECT_EQUAL(v1[1], 2);
+
+	Vector3i v2(data);
+	EXPECT_EQUAL(v2[0], 1);
+	EXPECT_EQUAL(v2[1], 2);
+	EXPECT_EQUAL(v2[2], 3);
+
+#if __cplusplus >= 201103L
+	LinAl::Vector<int, 1> v3(1);
+	EXPECT_EQUAL(v3[0], 1);
+
+	LinAl::Vector<int, 5> v4(1, 2, 3, 4, 5);
+	EXPECT_EQUAL(v4[0], 1);
+	EXPECT_EQUAL(v4[4], 5);
+#endif
+}
+
 void VectorTests::component_aliases()
 {
 	Vector3i v;
-- 
2.43.0


From 8d7271ce74a1550c1987e94e1770cd9fe0aeb865 Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Sun, 2 Jun 2019 16:21:40 +0300
Subject: [PATCH 03/16] Improve readability of Matrix test cases with typedefs

---
 tests/matrix.cpp | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/tests/matrix.cpp b/tests/matrix.cpp
index 1b01cbc..6f57888 100644
--- a/tests/matrix.cpp
+++ b/tests/matrix.cpp
@@ -8,6 +8,12 @@ using namespace Msp;
 
 class MatrixTests: public Test::RegisteredTest<MatrixTests>
 {
+private:
+	typedef LinAl::SquareMatrix<double, 2> Matrix2d;
+	typedef LinAl::Matrix<double, 3, 2> Matrix3x2d;
+	typedef LinAl::Matrix<double, 2, 3> Matrix2x3d;
+	typedef LinAl::SquareMatrix<double, 3> Matrix3d;
+
 public:
 	MatrixTests();
 
@@ -22,7 +28,7 @@ private:
 MatrixTests::MatrixTests()
 {
 	add(&MatrixTests::multiply, "Multiplication");
-	add(&MatrixTests::invert,   "Inversion");
+	add(&MatrixTests::invert, "Inversion");
 }
 
 void MatrixTests::multiply()
@@ -35,10 +41,10 @@ void MatrixTests::multiply()
 		9, 3, 9, 6
 	};
 
-	LinAl::Matrix<double, 3, 2> a(data);
-	LinAl::Matrix<double, 2, 3> b(data+6);
-	EXPECT_EQUAL(a*b, (LinAl::Matrix<double, 3, 3>(data+12)));
-	EXPECT_EQUAL(b*a, (LinAl::Matrix<double, 2, 2>(data+21)));
+	Matrix3x2d a(data);
+	Matrix2x3d b(data+6);
+	EXPECT_EQUAL(a*b, Matrix3d(data+12));
+	EXPECT_EQUAL(b*a, Matrix2d(data+21));
 }
 
 template<typename T, unsigned N>
@@ -59,7 +65,7 @@ void MatrixTests::invert()
 		1, 4, 2, 4, 8, 6, 2, 2, 4
 	};
 
-	LinAl::SquareMatrix<double, 3> m(data);
-	LinAl::SquareMatrix<double, 3> i = LinAl::invert(m);
+	Matrix3d m(data);
+	Matrix3d i = LinAl::invert(m);
 	EXPECT(is_identity(i*m));
 }
-- 
2.43.0


From 44a1fb77d9199ce76fec21ee8aaa6df38e141bf0 Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Sun, 2 Jun 2019 16:22:24 +0300
Subject: [PATCH 04/16] Copy Matrix test cases for DynamicMatrix

---
 tests/dynamicmatrix.cpp | 67 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 67 insertions(+)
 create mode 100644 tests/dynamicmatrix.cpp

diff --git a/tests/dynamicmatrix.cpp b/tests/dynamicmatrix.cpp
new file mode 100644
index 0000000..769945d
--- /dev/null
+++ b/tests/dynamicmatrix.cpp
@@ -0,0 +1,67 @@
+#include <limits>
+#include <msp/linal/dynamicmatrix.h>
+#include <msp/test/test.h>
+
+using namespace std;
+using namespace Msp;
+
+class DynamicMatrixTests: public Test::RegisteredTest<DynamicMatrixTests>
+{
+private:
+	typedef LinAl::DynamicMatrix<double> Matrixd;
+
+public:
+	DynamicMatrixTests();
+
+	static const char *get_name() { return "DynamicMatrix"; }
+
+private:
+	void multiply();
+	void invert();
+};
+
+
+DynamicMatrixTests::DynamicMatrixTests()
+{
+	add(&DynamicMatrixTests::multiply, "Multiplication");
+	add(&DynamicMatrixTests::invert, "Inversion");
+}
+
+void DynamicMatrixTests::multiply()
+{
+	static double data[] =
+	{
+		1, 2, 1, 2, 1, 2,
+		3, 0, 3, 0, 0, 3,
+		3, 6, 3, 3, 6, 3, 6, 3, 6,
+		9, 3, 9, 6
+	};
+
+	Matrixd a(3, 2, data);
+	Matrixd b(2, 3, data+6);
+	EXPECT_EQUAL(a*b, Matrixd(3, 3, data+12));
+	EXPECT_EQUAL(b*a, Matrixd(2, 2, data+21));
+}
+
+template<typename T>
+bool is_identity(const LinAl::DynamicMatrix<T> &m)
+{
+	static const T limit = numeric_limits<T>::epsilon()*4;
+	for(unsigned i=0; i<m.rows(); ++i)
+		for(unsigned j=0; j<m.columns(); ++j)
+			if(abs(T(i==j)-m(i, j))>limit)
+				return false;
+	return true;
+}
+
+void DynamicMatrixTests::invert()
+{
+	static double data[] =
+	{
+		1, 4, 2, 4, 8, 6, 2, 2, 4
+	};
+
+	Matrixd m(3, 3, data);
+	Matrixd i = LinAl::invert(m);
+	EXPECT(is_identity(i*m));
+}
-- 
2.43.0


From 459c248db8804b2a51010c89b72d3a0d326267e9 Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Sun, 2 Jun 2019 16:23:08 +0300
Subject: [PATCH 05/16] Fix a memory access error in DynamicMatrix::invert

Gauss_jordan returns a reference to the result, which is a local
variable.  It needs to be assigned to *this to retain it.
---
 source/linal/dynamicmatrix.h | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/source/linal/dynamicmatrix.h b/source/linal/dynamicmatrix.h
index 6d45657..1b191a1 100644
--- a/source/linal/dynamicmatrix.h
+++ b/source/linal/dynamicmatrix.h
@@ -308,7 +308,8 @@ inline DynamicMatrix<T> &DynamicMatrix<T>::invert()
 	for(unsigned i=0; i<rows_; ++i)
 		r(i, i) = T(1);
 
-	return gauss_jordan(*this, r);
+	gauss_jordan(*this, r);
+	return *this = r;
 }
 
 template<typename T>
-- 
2.43.0


From 4b895388d12b160f486c8b49b3296b15eeed7cc2 Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Sun, 2 Jun 2019 19:13:14 +0300
Subject: [PATCH 06/16] Add an interpolation sub-library

Initially including polynomials and splines, with cubic Hermite splines
being the only implemented spline type.
---
 Build                              |   1 +
 source/interpolate/dummy.cpp       |   3 +
 source/interpolate/hermitespline.h | 109 +++++++++++++++
 source/interpolate/polynomial.h    | 168 +++++++++++++++++++++++
 source/interpolate/spline.h        | 206 +++++++++++++++++++++++++++++
 5 files changed, 487 insertions(+)
 create mode 100644 source/interpolate/dummy.cpp
 create mode 100644 source/interpolate/hermitespline.h
 create mode 100644 source/interpolate/polynomial.h
 create mode 100644 source/interpolate/spline.h

diff --git a/Build b/Build
index 20875c4..ff363a4 100644
--- a/Build
+++ b/Build
@@ -4,6 +4,7 @@ package "mspmath"
 
 	library "mspmath"
 	{
+		source "source/interpolate";
 		source "source/linal";
 		source "source/geometry";
 		install true;
diff --git a/source/interpolate/dummy.cpp b/source/interpolate/dummy.cpp
new file mode 100644
index 0000000..ec10aec
--- /dev/null
+++ b/source/interpolate/dummy.cpp
@@ -0,0 +1,3 @@
+#include "hermitespline.h"
+#include "polynomial.h"
+#include "spline.h"
diff --git a/source/interpolate/hermitespline.h b/source/interpolate/hermitespline.h
new file mode 100644
index 0000000..057a9c4
--- /dev/null
+++ b/source/interpolate/hermitespline.h
@@ -0,0 +1,109 @@
+#ifndef MSP_INTERPOLATE_HERMITESPLINE_H_
+#define MSP_INTERPOLATE_HERMITESPLINE_H_
+
+#include "spline.h"
+
+namespace Msp {
+namespace Interpolate {
+
+template<typename T>
+struct Hermite
+{
+	static const Polynomial<T, 3> &h00();
+	static const Polynomial<T, 3> &h10();
+	static const Polynomial<T, 3> &h01();
+	static const Polynomial<T, 3> &h11();
+	static Polynomial<T, 3> make(T, T, T, T);
+};
+
+/**
+Implements a cubic Hermite spline.  It is parametrized by the value and slope
+at each knot.  The first derivative is guaranteed to be continuous, making it
+particularly useful as an interpolator.
+*/
+template<typename T, unsigned N = 1>
+class HermiteSpline: public Spline<T, 3, N>
+{
+public:
+	using typename Spline<T, 3, N>::Value;
+
+	struct Knot: Spline<T, 3, N>::Knot
+	{
+		Value dy;
+
+		Knot(T x_, Value y_, Value d): Spline<T, 3, N>::Knot(x_, y_), dy(d) { }
+	};
+
+	HermiteSpline(const std::vector<Knot> &);
+	HermiteSpline(Value, Value, Value, Value);
+};
+
+template<typename T>
+inline const Polynomial<T, 3> &Hermite<T>::h00()
+{
+	static Polynomial<T, 3> p(LinAl::Vector<T, 4>(2, -3, 0, 1));
+	return p;
+}
+
+template<typename T>
+inline const Polynomial<T, 3> &Hermite<T>::h10()
+{
+	static Polynomial<T, 3> p(LinAl::Vector<T, 4>(1, -2, 1, 0));
+	return p;
+}
+
+template<typename T>
+inline const Polynomial<T, 3> &Hermite<T>::h01()
+{
+	static Polynomial<T, 3> p(LinAl::Vector<T, 4>(-2, 3, 0, 0));
+	return p;
+}
+
+template<typename T>
+inline const Polynomial<T, 3> &Hermite<T>::h11()
+{
+	static Polynomial<T, 3> p(LinAl::Vector<T, 4>(1, -1, 0, 0));
+	return p;
+}
+
+template<typename T>
+inline Polynomial<T, 3> Hermite<T>::make(T p0, T m0, T p1, T m1)
+{
+	return h00()*p0 + h10()*m0 + h01()*p1 + h11()*m1;
+}
+
+
+template<typename T, unsigned N>
+inline HermiteSpline<T, N>::HermiteSpline(const std::vector<Knot> &k):
+	Spline<T, 3, N>(k.front())
+{
+	typedef SplineValue<T, N> SV;
+
+	this->reserve(k.size()-1);
+	for(unsigned i=1; i<k.size(); ++i)
+	{
+		T dx = k[i].x-k[i-1].x;
+		Polynomial<T, 1> t(LinAl::Vector<T, 2>(T(1)/dx, -k[i-1].x/dx));
+		Polynomial<T, 3> p[N];
+		for(unsigned j=0; j<N; ++j)
+			p[j] = Hermite<T>::make(SV::get(k[i-1].y, j), SV::get(k[i-1].dy, j)*dx, SV::get(k[i].y, j), SV::get(k[i].dy, j)*dx)(t);
+		this->add_segment(p, k[i].x);
+	}
+}
+
+template<typename T, unsigned N>
+inline HermiteSpline<T, N>::HermiteSpline(Value p0, Value m0, Value p1, Value m1):
+	Spline<T, 3, N>(Knot(0, p0, m0))
+{
+	typedef SplineValue<T, N> SV;
+
+	this->reserve(1);
+	Polynomial<T, 3> p[N];
+	for(unsigned i=0; i<N; ++i)
+		p[i] = Hermite<T>::make(SV::get(p0, i), SV::get(m0, i), SV::get(p1, i), SV::get(m1, i));
+	this->add_segment(p, T(1));
+}
+} // namespace Interpolate
+} // namespace Msp
+
+#endif
diff --git a/source/interpolate/polynomial.h b/source/interpolate/polynomial.h
new file mode 100644
index 0000000..859eea8
--- /dev/null
+++ b/source/interpolate/polynomial.h
@@ -0,0 +1,168 @@
+#ifndef MSP_INTERPOLATE_POLYNOMIAL_H_
+#define MSP_INTERPOLATE_POLYNOMIAL_H_
+
+#include <algorithm>
+#include <msp/linal/vector.h>
+
+namespace Msp {
+namespace Interpolate {
+
+/**
+Represents a polynomial of degree D with a single variable.  The coefficients
+are stored as a vector with the highest degree term first.
+*/
+template<typename T, unsigned D>
+class Polynomial
+{
+private:
+	LinAl::Vector<T, D+1> coeff;
+
+public:
+	Polynomial() { }
+	Polynomial(const LinAl::Vector<T, D+1> &c): coeff(c) { }
+	template<typename U, unsigned G>
+	Polynomial(const Polynomial<U, G> &);
+
+	unsigned degree() const { return D; }
+
+	LinAl::Vector<T, D+1> coefficients() { return coeff; }
+	const LinAl::Vector<T, D+1> coefficients() const { return coeff; }
+	T &coefficient(unsigned i) { return coeff[i]; }
+	const T &coefficient(unsigned i) const { return coeff[i]; }
+
+	Polynomial &operator+=(const Polynomial &);
+	Polynomial &operator-=(const Polynomial &);
+	Polynomial &operator*=(T);
+	Polynomial &operator/=(T);
+
+	/// Calculates the value of the polynomial at a particular point.
+	T operator()(T) const;
+
+	/** Substitutes the variable of the polynomial with another polynomial.  The
+	result is a new polynomial. */
+	template<unsigned G>
+	Polynomial<T, D*G> operator()(const Polynomial<T, G> &) const;
+};
+
+template<typename T, unsigned D>
+template<typename U, unsigned G>
+inline Polynomial<T, D>::Polynomial(const Polynomial<U, G> &other)
+{
+	char source_cannot_be_larger[D>=G ? 1 : -1];
+	for(unsigned i=0; i<=G; ++i)
+		coeff[D-i] = other.coefficient(G-i);
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> &Polynomial<T, D>::operator+=(const Polynomial<T, D> &p)
+{
+	coeff += p.coeff;
+	return *this;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> operator+(const Polynomial<T, D> &p, const Polynomial<T, D> &q)
+{
+	Polynomial<T, D> r(p);
+	return r += q;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> &Polynomial<T, D>::operator-=(const Polynomial<T, D> &p)
+{
+	coeff -= p.coeff;
+	return *this;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> operator-(const Polynomial<T, D> &p, const Polynomial<T, D> &q)
+{
+	Polynomial<T, D> r(p);
+	return r -= q;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> &Polynomial<T, D>::operator*=(T s)
+{
+	coeff *= s;
+	return *this;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> operator*(const Polynomial<T, D> &p, T s)
+{
+	Polynomial<T, D> r(p);
+	return r *= s;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> operator*(T s, const Polynomial<T, D> &p)
+{
+	return p*s;
+}
+
+template<typename T, unsigned D, unsigned G>
+inline Polynomial<T, D+G> operator*(const Polynomial<T, D> &p, const Polynomial<T, G> &q)
+{
+	LinAl::Vector<T, D+G+1> coeff;
+	for(unsigned i=0; i<=D; ++i)
+		for(unsigned j=0; j<=G; ++j)
+			coeff[i+j] += p.coefficient(i)*q.coefficient(j);
+	return Polynomial<T, D+G>(coeff);
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> &Polynomial<T, D>::operator/=(T s)
+{
+	coeff /= s;
+	return *this;
+}
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> operator/(const Polynomial<T, D> &p, T s)
+{
+	Polynomial<T, D> r(p);
+	return r /= s;
+}
+
+template<typename T, unsigned D>
+inline T Polynomial<T, D>::operator()(T x) const
+{
+	T r = T();
+	T x_i = T(1);
+	for(int i=D; i>=0; --i, x_i*=x)
+		r += coeff[i]*x_i;
+	return r;
+}
+
+template<typename T, unsigned D>
+template<unsigned G>
+inline Polynomial<T, D*G> Polynomial<T, D>::operator()(const Polynomial<T, G> &p) const
+{
+	Polynomial<T, D*G> r;
+	r.coefficient(D*G) = coeff[D];
+
+	LinAl::Vector<T, D*G+1> cf;
+	cf[0] = T(1);
+
+	for(unsigned i=1; i<=D; ++i)
+	{
+		for(int j=i*G; j>=0; --j)
+		{
+			T c = T();
+			for(unsigned k=0; (k<=static_cast<unsigned>(j) && k<=G); ++k)
+				c += cf[j-k]*p.coefficient(G-k);
+			cf[j] = c;
+		}
+
+		for(unsigned j=0; j<=i*G; ++j)
+			r.coefficient(D*G-j) += coeff[D-i]*cf[j];
+	}
+
+	return r;
+}
+
+} // namespace Interpolate
+} // namespace Msp
+
+#endif
diff --git a/source/interpolate/spline.h b/source/interpolate/spline.h
new file mode 100644
index 0000000..67c8df0
--- /dev/null
+++ b/source/interpolate/spline.h
@@ -0,0 +1,206 @@
+#ifndef MSP_INTERPOLATE_SPLINE_H_
+#define MSP_INTERPOLATE_SPLINE_H_
+
+#include <algorithm>
+#include <vector>
+#include "polynomial.h"
+
+namespace Msp {
+namespace Interpolate {
+
+template<typename T, unsigned N>
+struct SplineValue
+{
+	typedef LinAl::Vector<T, N> Type;
+	static T get(const Type &v, unsigned i) { return v[i]; }
+	static Type make(const T *v) { return Type(v); }
+};
+
+template<typename T>
+struct SplineValue<T, 1>
+{
+	typedef T Type;
+	static T get(const Type &v, unsigned) { return v; }
+	static Type make(const T *v) { return *v; }
+};
+
+/**
+Stores a spline of degree D.  It is a piecewise polynomial function with value
+continuity.  Derivatives are not guaranteed to be continuous.
+
+The template parameter N controls the dimensionality of the spline's values.
+Multi-dimensional splines can be used to implement parametric curves.
+
+While this class contains everything needed to store and evaluate splines, it
+cannot be used to construct non-empty splines.  See also HermiteSpline.
+*/
+template<typename T, unsigned D, unsigned N = 1>
+class Spline
+{
+public:
+	typedef typename SplineValue<T, N>::Type Value;
+
+	struct Knot
+	{
+		T x;
+		Value y;
+
+		Knot(): x(T()) { }
+		Knot(T x_, const Value &y_): x(x_), y(y_) { }
+	};
+
+	struct Segment
+	{
+		Knot start_knot;
+		Polynomial<T, D> polynomials[N];
+		Knot end_knot;
+
+		Value operator()(T) const;
+	};
+
+private:
+	enum { STRIDE = sizeof(Segment)-sizeof(Knot) };
+
+	unsigned short n_segments;
+	unsigned short capacity;
+	char *segments;
+
+public:
+	Spline();
+protected:
+	Spline(const Knot &);
+public:
+	Spline(const Spline &);
+	Spline &operator=(const Spline &);
+	~Spline();
+
+private:
+	static unsigned data_size(unsigned n) { return sizeof(Knot)+n*STRIDE; }
+protected:
+	void reserve(unsigned);
+	void add_segment(const Polynomial<T, D> [N], T);
+
+public:
+	unsigned degree() const { return D; }
+	unsigned size() const { return n_segments; }
+	const Segment &segment(unsigned i) const;
+	const Knot &knot(unsigned i) const { return i==0 ? segment(0).start_knot : segment(i-1).end_knot; }
+	const Polynomial<T, D> &polynomial(unsigned i, unsigned j = 0) const { return segment(i).polynomials[j]; }
+
+	Value operator()(T) const;
+};
+
+template<typename T, unsigned D, unsigned N>
+inline Spline<T, D, N>::Spline():
+	n_segments(0),
+	capacity(0),
+	segments(0)
+{ }
+
+template<typename T, unsigned D, unsigned N>
+inline Spline<T, D, N>::Spline(const Knot &k):
+	n_segments(0),
+	capacity(0),
+	segments(new char[sizeof(Knot)])
+{
+	new(segments) Knot(k);
+}
+
+template<typename T, unsigned D, unsigned N>
+inline Spline<T, D, N>::Spline(const Spline &s):
+	n_segments(0),
+	capacity(0),
+	segments(0)
+{
+	*this = s;
+}
+
+template<typename T, unsigned D, unsigned N>
+inline Spline<T, D, N> &Spline<T, D, N>::operator=(const Spline &s)
+{
+	if(s.segments)
+	{
+		reserve(s.n_segments);
+		std::copy(s.segments, s.segments+data_size(s.n_segments), segments);
+		n_segments = s.n_segments;
+	}
+	else
+	{
+		delete[] segments;
+		n_segments = 0;
+		capacity = 0;
+		segments = 0;
+	}
+}
+
+template<typename T, unsigned D, unsigned N>
+inline Spline<T, D, N>::~Spline()
+{
+	delete[] segments;
+}
+
+template<typename T, unsigned D, unsigned N>
+inline void Spline<T, D, N>::reserve(unsigned n)
+{
+	if(n<capacity)
+		return;
+
+	char *new_segs = new char[data_size(n)];
+	if(segments)
+	{
+		std::copy(segments, segments+data_size(n_segments), new_segs);
+		delete[] segments;
+	}
+	segments = new_segs;
+	capacity = n;
+}
+
+template<typename T, unsigned D, unsigned N>
+inline void Spline<T, D, N>::add_segment(const Polynomial<T, D> polys[N], T end_x)
+{
+	if(!segments)
+		throw std::logic_error("Spline::add_segment");
+
+	if(n_segments==capacity)
+		reserve(n_segments+1);
+	++n_segments;
+
+	Segment &seg = *reinterpret_cast<Segment *>(segments+(n_segments-1)*STRIDE);
+	std::copy(polys, polys+N, seg.polynomials);
+	seg.end_knot.x = end_x;
+	seg.end_knot.y = seg(end_x);
+}
+
+template<typename T, unsigned D, unsigned N>
+inline const typename Spline<T, D, N>::Segment &Spline<T, D, N>::segment(unsigned i) const
+{
+	return *reinterpret_cast<const Segment *>(segments+i*STRIDE);
+}
+
+template<typename T, unsigned D, unsigned N>
+inline typename Spline<T, D, N>::Value Spline<T, D, N>::operator()(T x) const
+{
+	if(!n_segments)
+		return Value();
+
+	const Segment *seg = &segment(0);
+	for(unsigned i=1; (i<n_segments && x>seg->end_knot.x); ++i)
+		seg = &segment(i);
+
+	return (*seg)(x);
+}
+
+template<typename T, unsigned D, unsigned N>
+inline typename Spline<T, D, N>::Value Spline<T, D, N>::Segment::operator()(T x) const
+{
+	T v[N];
+	for(unsigned i=0; i<N; ++i)
+		v[i] = polynomials[i](x);
+
+	return SplineValue<T, N>::make(v);
+}
+
+} // namespace Interpolate
+} // namespace Msp
+
+#endif
-- 
2.43.0


From 23583f28679d722e7354c1aa9951731e3dc79726 Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Tue, 4 Jun 2019 16:37:36 +0300
Subject: [PATCH 07/16] Move the Knot struct out of Spline

It doesn't depend on the degree of the spline, so this will generate less
different types and make it easier to create splines of different degrees
from the same knots.
---
 source/interpolate/spline.h | 22 +++++++++++++---------
 1 file changed, 13 insertions(+), 9 deletions(-)

diff --git a/source/interpolate/spline.h b/source/interpolate/spline.h
index 67c8df0..c801188 100644
--- a/source/interpolate/spline.h
+++ b/source/interpolate/spline.h
@@ -3,6 +3,7 @@
 
 #include <algorithm>
 #include <vector>
+#include "knot.h"
 #include "polynomial.h"
 
 namespace Msp {
@@ -24,6 +25,17 @@ struct SplineValue<T, 1>
 	static Type make(const T *v) { return *v; }
 };
 
+template<typename T, unsigned N>
+struct SplineKnot
+{
+	typedef typename SplineValue<T, N>::Type Value;
+	T x;
+	Value y;
+
+	SplineKnot(): x(T()) { }
+	SplineKnot(T x_, const Value &y_): x(x_), y(y_) { }
+};
+
 /**
 Stores a spline of degree D.  It is a piecewise polynomial function with value
 continuity.  Derivatives are not guaranteed to be continuous.
@@ -39,15 +51,7 @@ class Spline
 {
 public:
 	typedef typename SplineValue<T, N>::Type Value;
-
-	struct Knot
-	{
-		T x;
-		Value y;
-
-		Knot(): x(T()) { }
-		Knot(T x_, const Value &y_): x(x_), y(y_) { }
-	};
+	typedef SplineKnot<T, N> Knot;
 
 	struct Segment
 	{
-- 
2.43.0


From 1d33efd177208ed449b6333539bf13e4ffc93b97 Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Wed, 5 Jun 2019 22:24:35 +0300
Subject: [PATCH 08/16] Move a #include to where it's actually used

---
 source/interpolate/hermitespline.h | 1 +
 source/interpolate/spline.h        | 1 -
 2 files changed, 1 insertion(+), 1 deletion(-)

diff --git a/source/interpolate/hermitespline.h b/source/interpolate/hermitespline.h
index 057a9c4..509aef1 100644
--- a/source/interpolate/hermitespline.h
+++ b/source/interpolate/hermitespline.h
@@ -1,6 +1,7 @@
 #ifndef MSP_INTERPOLATE_HERMITESPLINE_H_
 #define MSP_INTERPOLATE_HERMITESPLINE_H_
 
+#include <vector>
 #include "spline.h"
 
 namespace Msp {
diff --git a/source/interpolate/spline.h b/source/interpolate/spline.h
index c801188..88626d2 100644
--- a/source/interpolate/spline.h
+++ b/source/interpolate/spline.h
@@ -2,7 +2,6 @@
 #define MSP_INTERPOLATE_SPLINE_H_
 
 #include <algorithm>
-#include <vector>
 #include "knot.h"
 #include "polynomial.h"
 
-- 
2.43.0


From acdb569a0370162ffc282b2e43221f34e35e7c58 Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Wed, 5 Jun 2019 22:25:02 +0300
Subject: [PATCH 09/16] Fix a missing return in assignment operator

---
 source/interpolate/spline.h | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/source/interpolate/spline.h b/source/interpolate/spline.h
index 88626d2..d761979 100644
--- a/source/interpolate/spline.h
+++ b/source/interpolate/spline.h
@@ -134,6 +134,8 @@ inline Spline<T, D, N> &Spline<T, D, N>::operator=(const Spline &s)
 		capacity = 0;
 		segments = 0;
 	}
+
+	return *this;
 }
 
 template<typename T, unsigned D, unsigned N>
-- 
2.43.0


From 0d4dcdc1c4717c88e79dfd0b7d4ebe34ff7cc548 Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Thu, 6 Jun 2019 00:22:10 +0300
Subject: [PATCH 10/16] Add a simple linearly interpolated "spline"

---
 source/interpolate/linearspline.h | 42 +++++++++++++++++++++++++++++++
 1 file changed, 42 insertions(+)
 create mode 100644 source/interpolate/linearspline.h

diff --git a/source/interpolate/linearspline.h b/source/interpolate/linearspline.h
new file mode 100644
index 0000000..2c88ceb
--- /dev/null
+++ b/source/interpolate/linearspline.h
@@ -0,0 +1,42 @@
+#ifndef MSP_INTERPOLATE_LINEARSPLINE_H_
+#define MSP_INTERPOLATE_LINEARSPLINE_H_
+
+#include <vector>
+#include "spline.h"
+
+namespace Msp {
+namespace Interpolate {
+
+/**
+A very simple type of spline.  It interpolates the value linearly between
+knots.
+*/
+template<typename T, unsigned N = 1>
+class LinearSpline: public Spline<T, 1, N>
+{
+public:
+	using typename Spline<T, 1, N>::Value;
+	using typename Spline<T, 1, N>::Knot;
+
+	LinearSpline(const std::vector<Knot> &k):
+		Spline<T, 1, N>(k.front())
+	{
+		typedef SplineValue<T, N> SV;
+
+		this->reserve(k.size()-1);
+		for(unsigned i=1; i<k.size(); ++i)
+		{
+			T dx = k[i].x-k[i-1].x;
+			Value slope = (k[i].y-k[i-1].y)/dx;
+			Polynomial<T, 1> p[N];
+			for(unsigned j=0; j<N; ++j)
+				p[j] = Polynomial<T, 1>(LinAl::Vector<T, 2>(SV::get(slope, j), SV::get(k[i-1].y, j)-k[i-1].x*SV::get(slope, j)));
+			this->add_segment(p, k[i].x);
+		}
+	}
+};
+
+} // namespace Interpolate
+} // namespace Msp
+
+#endif
-- 
2.43.0


From 92480ec2f41502e3cbdfe1fe3ca703c4521ff747 Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Thu, 6 Jun 2019 16:01:39 +0300
Subject: [PATCH 11/16] Remove a stray include

I thought of moving SplineKnot to a separate header but decided against
it (for now) and this got left behind.
---
 source/interpolate/spline.h | 1 -
 1 file changed, 1 deletion(-)

diff --git a/source/interpolate/spline.h b/source/interpolate/spline.h
index d761979..bd606d9 100644
--- a/source/interpolate/spline.h
+++ b/source/interpolate/spline.h
@@ -2,7 +2,6 @@
 #define MSP_INTERPOLATE_SPLINE_H_
 
 #include <algorithm>
-#include "knot.h"
 #include "polynomial.h"
 
 namespace Msp {
-- 
2.43.0


From c854dd42c57b0430558baed5151449c83e24507f Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Sat, 8 Jun 2019 02:16:06 +0300
Subject: [PATCH 12/16] Add a bezier spline type

---
 source/interpolate/bezierspline.h | 88 +++++++++++++++++++++++++++++++
 source/interpolate/dummy.cpp      |  2 +
 2 files changed, 90 insertions(+)
 create mode 100644 source/interpolate/bezierspline.h

diff --git a/source/interpolate/bezierspline.h b/source/interpolate/bezierspline.h
new file mode 100644
index 0000000..c6e295a
--- /dev/null
+++ b/source/interpolate/bezierspline.h
@@ -0,0 +1,88 @@
+#ifndef MSP_INTERPOLATE_BEZIERSPLINE_H_
+#define MSP_INTERPOLATE_BEZIERSPLINE_H_
+
+#include <vector>
+#include "spline.h"
+
+namespace Msp {
+namespace Interpolate {
+
+template<typename T, unsigned D>
+struct Bernstein
+{
+	static Polynomial<T, D> b(unsigned);
+};
+
+/**
+A spline composed of bezier curves.  Each segment has D-1 control points
+between the knots which determine the shape of the curve.  Bezier splines
+with 2- or 3-dimensional values are commonly used in graphics.
+*/
+template<typename T, unsigned D, unsigned N = 1>
+class BezierSpline: public Spline<T, D, N>
+{
+public:
+	using typename Spline<T, D, N>::Knot;
+
+	BezierSpline(const std::vector<Knot> &);
+};
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> Bernstein<T, D>::b(unsigned k)
+{
+	if(k>D)
+		throw std::invalid_argument("Bernstein::b");
+
+	LinAl::Vector<T, D+1> coeff;
+
+	int c = ((D-k)%2 ? -1 : 1);
+	for(unsigned i=0; i<=D-k; ++i)
+	{
+		coeff[i] = c;
+		c = c*-1*static_cast<int>(D-k-i)/static_cast<int>(i+1);
+	}
+
+	unsigned n = 1;
+	for(unsigned i=0; i<k; ++i)
+		n *= D-i;
+	for(unsigned i=2; i<=k; ++i)
+		n /= i;
+
+	for(unsigned i=0; i<=D; ++i)
+		coeff[i] *= n;
+
+	return Polynomial<T, D>(coeff);
+}
+
+template<typename T, unsigned D, unsigned N>
+inline BezierSpline<T, D, N>::BezierSpline(const std::vector<Knot> &k):
+	Spline<T, D, N>(k.front())
+{
+	typedef SplineValue<T, N> SV;
+
+	if((k.size()-1)%D)
+		throw std::invalid_argument("BezierSpline::BezierSpline");
+
+	Polynomial<T, D> b[D+1];
+	for(unsigned i=0; i<=D; ++i)
+		b[i] = Bernstein<T, D>::b(i);
+
+	for(unsigned i=0; i+D<k.size(); i+=D)
+	{
+		T dx = k[i+D].x-k[i].x;
+		Polynomial<T, 1> t(LinAl::Vector<T, 2>(T(1)/dx, -k[i].x/dx));
+		Polynomial<T, D> p[N];
+		for(unsigned j=0; j<N; ++j)
+		{
+			for(unsigned c=0; c<=D; ++c)
+				p[j] += b[c]*SV::get(k[i+c].y, j);
+			p[j] = p[j](t);
+		}
+		this->add_segment(p, k[i+D].x);
+	}
+}
+
+} // namespace Interpolate
+} // namespace Msp
+
+#endif
diff --git a/source/interpolate/dummy.cpp b/source/interpolate/dummy.cpp
index ec10aec..cfa902e 100644
--- a/source/interpolate/dummy.cpp
+++ b/source/interpolate/dummy.cpp
@@ -1,3 +1,5 @@
+#include "bezierspline.h"
 #include "hermitespline.h"
+#include "linearspline.h"
 #include "polynomial.h"
 #include "spline.h"
-- 
2.43.0


From febdd901857c88bcd4cd1be7abc4e6ae766e8a4a Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Wed, 30 Dec 2020 17:14:55 +0200
Subject: [PATCH 13/16] Support Matrix multiplication operators on SquareMatrix

---
 source/linal/squarematrix.h | 1 +
 1 file changed, 1 insertion(+)

diff --git a/source/linal/squarematrix.h b/source/linal/squarematrix.h
index b8531aa..5b9415a 100644
--- a/source/linal/squarematrix.h
+++ b/source/linal/squarematrix.h
@@ -24,6 +24,7 @@ public:
 	static SquareMatrix identity();
 
 	SquareMatrix &operator*=(const SquareMatrix &);
+	using Matrix<T, S, S>::operator*=;
 
 	SquareMatrix &invert();
 };
-- 
2.43.0


From 0202ce51b7a16a4b8cbc9d3257f46ade114eba34 Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Sun, 13 Mar 2022 21:29:31 +0200
Subject: [PATCH 14/16] Assign a version number

---
 Build | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/Build b/Build
index ff363a4..f183ff5 100644
--- a/Build
+++ b/Build
@@ -1,5 +1,7 @@
 package "mspmath"
 {
+	version "1.0";
+
 	require "mspdatafile";
 
 	library "mspmath"
-- 
2.43.0


From 291c11cf66e7083dc21fffea1afbdeaaad96077d Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Sun, 13 Mar 2022 21:30:11 +0200
Subject: [PATCH 15/16] Require C++11 for building

---
 Build                 |  5 +++++
 source/linal/vector.h | 35 +----------------------------------
 2 files changed, 6 insertions(+), 34 deletions(-)

diff --git a/Build b/Build
index f183ff5..7fa961c 100644
--- a/Build
+++ b/Build
@@ -4,6 +4,11 @@ package "mspmath"
 
 	require "mspdatafile";
 
+	build_info
+	{
+		standard CXX "c++11";
+	};
+
 	library "mspmath"
 	{
 		source "source/interpolate";
diff --git a/source/linal/vector.h b/source/linal/vector.h
index 31597d4..610ac21 100644
--- a/source/linal/vector.h
+++ b/source/linal/vector.h
@@ -84,14 +84,9 @@ public:
 	use by Matrix row accessor. */
 	Vector(const T *, unsigned);
 
-#if __cplusplus >= 201103L
 	template<typename... Args>
 	Vector(T, Args...);
-#else
-	Vector(T, T);
-	Vector(T, T, T);
-	Vector(T, T, T, T);
-#endif
+
 	template<typename U>
 	Vector(const Vector<U, N> &);
 
@@ -131,7 +126,6 @@ inline Vector<T, N>::Vector(const T *d, unsigned stride)
 		(*this)[i] = d[i*stride];
 }
 
-#if __cplusplus >= 201103L
 template<typename T, unsigned N>
 template<typename... Args>
 inline Vector<T, N>::Vector(T x_, Args... v)
@@ -142,33 +136,6 @@ inline Vector<T, N>::Vector(T x_, Args... v)
 	for(auto c: std::initializer_list<T> { static_cast<T>(v)... })
 		(*this)[i++] = c;
 }
-#else
-/* The compiler won't instantiate these unless they are used.  Trying to use
-them on the wrong class results in an error. */
-template<typename T, unsigned N>
-inline Vector<T, N>::Vector(T x_, T y_)
-{
-	this->VectorComponents<T, 2>::x = x_;
-	this->VectorComponents<T, 2>::y = y_;
-}
-
-template<typename T, unsigned N>
-inline Vector<T, N>::Vector(T x_, T y_, T z_)
-{
-	this->VectorComponents<T, 3>::x = x_;
-	this->VectorComponents<T, 3>::y = y_;
-	this->VectorComponents<T, 3>::z = z_;
-}
-
-template<typename T, unsigned N>
-inline Vector<T, N>::Vector(T x_, T y_, T z_, T w_)
-{
-	this->VectorComponents<T, 4>::x = x_;
-	this->VectorComponents<T, 4>::y = y_;
-	this->VectorComponents<T, 4>::z = z_;
-	this->VectorComponents<T, 4>::w = w_;
-}
-#endif
 
 template<typename T, unsigned N>
 template<typename U>
-- 
2.43.0


From ff8619c1c440fd1a7c8deaca763a33b8e1d53b5e Mon Sep 17 00:00:00 2001
From: Mikko Rasa <tdb@tdb.fi>
Date: Sun, 13 Mar 2022 21:32:47 +0200
Subject: [PATCH 16/16] Make gauss_jordan operate on columns instead of rows

Since matrices are stored in column-major order, this makes the algorithm
more amenable to possible simdification.
---
 source/linal/dynamicmatrix.h | 30 +++++++++++++++---------------
 source/linal/matrix.h        | 24 ++++++++++++------------
 source/linal/matrixops.h     | 28 ++++++++++++++--------------
 3 files changed, 41 insertions(+), 41 deletions(-)

diff --git a/source/linal/dynamicmatrix.h b/source/linal/dynamicmatrix.h
index 1b191a1..de4511c 100644
--- a/source/linal/dynamicmatrix.h
+++ b/source/linal/dynamicmatrix.h
@@ -45,9 +45,9 @@ public:
 	DynamicMatrix &operator+=(const DynamicMatrix &);
 	DynamicMatrix &operator-=(const DynamicMatrix &);
 
-	DynamicMatrix &exchange_rows(unsigned, unsigned);
-	DynamicMatrix &multiply_row(unsigned, T);
-	DynamicMatrix &add_row(unsigned, unsigned, T);
+	DynamicMatrix &exchange_columns(unsigned, unsigned);
+	DynamicMatrix &multiply_column(unsigned, T);
+	DynamicMatrix &add_column(unsigned, unsigned, T);
 
 	DynamicMatrix &invert();
 };
@@ -252,38 +252,38 @@ inline bool operator==(const DynamicMatrix<T> &m1, const DynamicMatrix<T> &m2)
 }
 
 template<typename T>
-inline DynamicMatrix<T> &DynamicMatrix<T>::exchange_rows(unsigned i, unsigned j)
+inline DynamicMatrix<T> &DynamicMatrix<T>::exchange_columns(unsigned i, unsigned j)
 {
-	if(i>=rows_ || j>=rows_)
-		throw std::out_of_range("DynamicMatrix::exchange_rows");
+	if(i>=columns_ || j>=columns_)
+		throw std::out_of_range("DynamicMatrix::exchange_columns");
 
 	using std::swap;
 	for(unsigned k=0; k<columns_; ++k)
-		swap(element(i, k), element(j, k));
+		swap(element(k, i), element(k, j));
 
 	return *this;
 }
 
 template<typename T>
-inline DynamicMatrix<T> &DynamicMatrix<T>::multiply_row(unsigned i, T s)
+inline DynamicMatrix<T> &DynamicMatrix<T>::multiply_column(unsigned i, T s)
 {
-	if(i>=rows_)
-		throw std::out_of_range("DynamicMatrix::multiply_row");
+	if(i>=columns_)
+		throw std::out_of_range("DynamicMatrix::multiply_column");
 
 	for(unsigned k=0; k<columns_; ++k)
-		element(i, k) *= s;
+		element(k, i) *= s;
 
 	return *this;
 }
 
 template<typename T>
-inline DynamicMatrix<T> &DynamicMatrix<T>::add_row(unsigned i, unsigned j, T s)
+inline DynamicMatrix<T> &DynamicMatrix<T>::add_column(unsigned i, unsigned j, T s)
 {
-	if(i>=rows_ || j>=rows_)
-		throw std::out_of_range("DynamicMatrix::exchange_rows");
+	if(i>=columns_ || j>=columns_)
+		throw std::out_of_range("DynamicMatrix::exchange_columns");
 
 	for(unsigned k=0; k<columns_; ++k)
-		element(j, k) += element(i, k)*s;
+		element(k, j) += element(k, i)*s;
 
 	return *this;
 }
diff --git a/source/linal/matrix.h b/source/linal/matrix.h
index 2711d40..bf21184 100644
--- a/source/linal/matrix.h
+++ b/source/linal/matrix.h
@@ -51,9 +51,9 @@ public:
 	Matrix &operator+=(const Matrix &);
 	Matrix &operator-=(const Matrix &);
 
-	Matrix &exchange_rows(unsigned, unsigned);
-	Matrix &multiply_row(unsigned, T);
-	Matrix &add_row(unsigned, unsigned, T);
+	Matrix &exchange_columns(unsigned, unsigned);
+	Matrix &multiply_column(unsigned, T);
+	Matrix &add_column(unsigned, unsigned, T);
 };
 
 template<typename T, unsigned M, unsigned N>
@@ -227,27 +227,27 @@ inline bool operator==(const Matrix<T, M, N> &a, const Matrix<T, M, N> &b)
 }
 
 template<typename T, unsigned M, unsigned N>
-inline Matrix<T, M, N> &Matrix<T, M, N>::exchange_rows(unsigned i, unsigned j)
+inline Matrix<T, M, N> &Matrix<T, M, N>::exchange_columns(unsigned i, unsigned j)
 {
 	using std::swap;
-	for(unsigned k=0; k<N; ++k)
-		swap(element(i, k), element(j, k));
+	for(unsigned k=0; k<M; ++k)
+		swap(element(k, i), element(k, j));
 	return *this;
 }
 
 template<typename T, unsigned M, unsigned N>
-inline Matrix<T, M, N> &Matrix<T, M, N>::multiply_row(unsigned i, T s)
+inline Matrix<T, M, N> &Matrix<T, M, N>::multiply_column(unsigned i, T s)
 {
-	for(unsigned k=0; k<N; ++k)
-		element(i, k) *= s;
+	for(unsigned k=0; k<M; ++k)
+		element(k, i) *= s;
 	return *this;
 }
 
 template<typename T, unsigned M, unsigned N>
-inline Matrix<T, M, N> &Matrix<T, M, N>::add_row(unsigned i, unsigned j, T s)
+inline Matrix<T, M, N> &Matrix<T, M, N>::add_column(unsigned i, unsigned j, T s)
 {
-	for(unsigned k=0; k<N; ++k)
-		element(j, k) += element(i, k)*s;
+	for(unsigned k=0; k<M; ++k)
+		element(k, j) += element(k, i)*s;
 	return *this;
 }
 
diff --git a/source/linal/matrixops.h b/source/linal/matrixops.h
index 84c8012..0988acc 100644
--- a/source/linal/matrixops.h
+++ b/source/linal/matrixops.h
@@ -22,37 +22,37 @@ inline T &gauss_jordan(T &m, T &r)
 	typedef typename T::ElementType V;
 	using std::abs;
 
-	for(unsigned i=0; i<m.rows(); ++i)
+	for(unsigned i=0; i<m.columns(); ++i)
 	{
 		unsigned pivot = i;
-		for(unsigned j=i+1; j<m.rows(); ++j)
-			if(abs(m.element(j, i))>abs(m.element(pivot, i)))
+		for(unsigned j=i+1; j<m.columns(); ++j)
+			if(abs(m.element(i, j))>abs(m.element(i, pivot)))
 				pivot = j;
 
-		if(m.element(pivot, i)==V(0))
+		if(m.element(i, pivot)==V(0))
 			throw not_invertible();
 
 		if(pivot!=i)
 		{
-			m.exchange_rows(i, pivot);
-			r.exchange_rows(i, pivot);
+			m.exchange_columns(i, pivot);
+			r.exchange_columns(i, pivot);
 		}
 
-		for(unsigned j=i+1; j<m.rows(); ++j)
+		for(unsigned j=i+1; j<m.columns(); ++j)
 		{
-			V a = -m.element(j, i)/m.element(i, i);
-			m.add_row(i, j, a);
-			r.add_row(i, j, a);
+			V a = -m.element(i, j)/m.element(i, i);
+			m.add_column(i, j, a);
+			r.add_column(i, j, a);
 		}
 
 		V a = V(1)/m.element(i, i);
-		m.multiply_row(i, a);
-		r.multiply_row(i, a);
+		m.multiply_column(i, a);
+		r.multiply_column(i, a);
 	}
 
-	for(unsigned i=m.rows(); i-->0; )
+	for(unsigned i=m.columns(); i-->0; )
 		for(unsigned j=i; j-->0; )
-			r.add_row(i, j, -m.element(j, i));
+			r.add_column(i, j, -m.element(i, j));
 
 	return r;
 }
-- 
2.43.0