units.cpp

template<int cols>
class Matrix {
	typedef double Row[cols];
	Row*const mat;
	const int r;
	Matrix(const Matrix&); // need copy constructor, but I'm cheating
public:
	int rows() const { return r; }
	Matrix(int rows):mat(new Row[rows]), r(rows) {}
	double& operator()(int row, int col) {
		if(row < 0 || row >= r || col <0 || col >= cols) throw "bug!";
		return mat[row][col];
	}
};

template<int rows>
struct TransposedMatrix {
	int cols() const { return mat.rows(); }
	TransposedMatrix(Matrix<rows>& m) : mat(m) {}
	double operator()(int row, int col) { return mat(col,row); }
private:
	Matrix<rows>& mat;
};

template<int n>
TransposedMatrix<n> transpose(Matrix<n>& m) { return m; }

template<int rows, int cols>
struct SmallMatrix {
	double mat[rows][cols];
	double& operator()(int row, int col) {
		if(row < 0 || row >= rows || col <0 || col >= cols) throw "bug!";
		return mat[row][col];
	}
	SmallMatrix(){
		for(int i=0; i<rows; i++)
			for(int j=0; j<cols; j++)
				mat[i][j] = 0;
	}
};

template<int rows, int cols>
SmallMatrix<rows,cols> operator*(TransposedMatrix<rows>& a, Matrix<cols>& b) {
	if(a.cols() != b.rows()) throw "bug!";
	return matmul<rows,cols>(a, b, a.cols());
}

template<int rows, int cols, int n>
SmallMatrix<rows,cols> operator*(SmallMatrix<rows, n>& a, SmallMatrix<n, cols>& b) {
	return matmul<rows,cols>(a, b, n);
}

template<int rows, int cols, typename M1, typename M2>
SmallMatrix<rows,cols> matmul(M1& a, M2& b, int n){
	SmallMatrix<rows, cols> c;
	for(int i=0; i<rows; i++)
		for(int j=0; j<cols; j++)
			for(int k=0; k<n; k++)
				c(i,j) += a(i,k) * b(k,j);
	return c;
}

template<int rows, int cols> // Gaussian elemination alg simplified from Wikipedia, with different variable names
void row_reduce(SmallMatrix<rows, cols>& mat) {
	int n = rows < cols ? rows : cols;
	for(int i=0; i<n; i++) {
		int p = i;
		while(p<n && mat(p,i)==0) p++;
		for(int j=0; j<cols; j++){
			double temp = mat(i,j);
			mat(i,j) = mat(p,j);
			mat(p,j) = temp;
		}
		double scale = mat(i,i);
		for(int j=i; j<cols; j++)
			mat(i,j) /= scale;
		for(int j=0; j<rows; j++) {
			if(j==i) continue;
			double factor = mat(j,i);
			for(int k=i; k<cols; k++)
				mat(j,k) -= mat(i,k) * factor;
		}
	}
}

template<int n>
SmallMatrix<n, n> inverse(SmallMatrix<n, n> mat) {
	SmallMatrix<n, 2*n> temp;
	for(int i=0; i<n; i++) // glue on the identity matrix
		for(int j=0; j<n; j++)
			temp(i,j) = mat(i,j),
			temp(i,j+n) = i == j;
	row_reduce(temp);
	SmallMatrix<n,n> ans;
	for(int i=0; i<n; i++)
		for(int j=0; j<n; j++)
			ans(i,j) = temp(i,j+n);
	return ans;
}

struct Line { double slope, intercept; };

Line least_squares_fit(const double x[], const double y[], int n) {
	Matrix<1> Y(n);
	for(int i=0; i<n; i++) Y(i,0) = y[i];
	Matrix<2> X(n);
	for(int i=0; i<n; i++)
		X(i,0) = x[i],
		X(i,1) = 1;
	TransposedMatrix<2> Xt = transpose(X);
	SmallMatrix<2,1> ans = inverse(Xt*X)*(Xt*Y);
	Line ret;
	ret.slope = ans(0,0);
	ret.intercept = ans(1,0);
	return ret;
}

#include<cstdio>
int main() {
	double x[] = {1,2,3}, y[] = {7,8,9}; // input
	Line line = least_squares_fit(x, y, 3);
	printf("y=%gx+%g\n", line.slope, line.intercept);
}