Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# /units.cpp

Created Apr 21, 2014
 template 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 struct TransposedMatrix { int cols() const { return mat.rows(); } TransposedMatrix(Matrix& m) : mat(m) {} double operator()(int row, int col) { return mat(col,row); } private: Matrix& mat; }; template TransposedMatrix transpose(Matrix& m) { return m; } template 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 SmallMatrix operator*(TransposedMatrix& a, Matrix& b) { if(a.cols() != b.rows()) throw "bug!"; return matmul(a, b, a.cols()); } template SmallMatrix operator*(SmallMatrix& a, SmallMatrix& b) { return matmul(a, b, n); } template SmallMatrix matmul(M1& a, M2& b, int n){ SmallMatrix c; for(int i=0; i // Gaussian elemination alg simplified from Wikipedia, with different variable names void row_reduce(SmallMatrix& mat) { int n = rows < cols ? rows : cols; for(int i=0; i SmallMatrix inverse(SmallMatrix mat) { SmallMatrix temp; for(int i=0; i ans; for(int i=0; i Y(n); for(int i=0; i X(n); for(int i=0; i 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 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); }
to join this conversation on GitHub. Already have an account? Sign in to comment