Matematické Fórum

/* Rychle spatlaná třída, obsahuje spoustu bugů - nejsou ošetřeny velikosti matic, apod. Taky jsem možná měl indexovad (row,col) místo (col,row), ale zatím mi to stačí */
#ifndef _GF_MATRIX_H
#define _GF_MATRIX_H

#include <iostream>

namespace Gf
{
    template<class Elem>
    class Matrix
    {
        protected:
            size_t _cols;
            size_t _rows;
            Elem *_d;
        public:
            Matrix(const Matrix& matrix)
            {
                _cols = matrix._cols;
                _rows = matrix._rows;
                _d = new Elem[_cols*_rows];

                for (size_t i = 0; i < _cols*_rows; ++i)
                {
                    _d[i] = matrix._d[i];
                }
            }
            
            Matrix(Matrix&& matrix)
            {
                _cols = matrix._cols;
                _rows = matrix._rows;
                _d = matrix._d;
                
                matrix._d = 0;
                matrix._cols = 0;
                matrix._rows = 0;
            }
            
            Matrix(size_t cols, size_t rows)
            {
                _cols = cols;
                _rows = rows;
                _d = new Elem[cols*rows];
            }
            
            Matrix(size_t cols, size_t rows, const Elem& val)
            {
                _cols = cols;
                _rows = rows;
                _d = new Elem[cols*rows];

                for (size_t i = 0; i < _cols*_rows; ++i)
                {
                    _d[i] = val;
                }
            }

            ~Matrix() { delete []_d; }

            Matrix& operator=(const Matrix& matrix)
            {
                if (&matrix == this) return *this;
                if (_cols != matrix._cols || _rows != matrix._rows)
                {
                    delete []_d;
                    _cols = matrix._cols;
                    _rows = matrix._rows;
                    _d = new Elem[_cols*_rows];
                }
                
                for (size_t i = 0; i < _cols*_rows; ++i)
                {
                    _d[i] = matrix._d[i];
                }
                return *this;
            }
            
            Matrix& operator=(Matrix&& matrix)
            {
                if (&matrix == this) return *this;
                delete []_d;
                _cols = matrix._cols;
                _rows = matrix._rows;
                _d = matrix._d;
                
                matrix._d = 0;
                matrix._cols = 0;
                matrix._rows = 0;
                return *this;
            }

            const Elem& operator()(size_t col, size_t row) const { return _d[_cols*row + col]; }
            Elem& operator()(size_t col, size_t row)             { return _d[_cols*row + col]; }

            size_t cols() const          { return _cols; }
            size_t rows() const          { return _rows; }
            Elem*  row(size_t row) const { return _d+_cols*row; }

            Matrix& operator/=(const Elem& elem)
            {
                for (size_t i = 0; i < _cols*_rows; ++i) _d[i] /= elem;
                return *this;
            }

            Matrix& operator*=(const Elem& elem)
            {
                for (size_t i = 0; i < _cols*_rows; ++i) _d[i] /= elem;
                return *this;
            }

            Matrix transpose() const
            {
                Matrix result(_rows, _cols);

                for (size_t i = 0; i < _cols; ++i)
                    for (size_t j = 0; j < _rows; ++j)
                        result(j, i) = (*this)(i, j);

                return result;
            }

            void setIdentity()
            {
                for (size_t i = 0; i < _cols; ++i)
                    for (size_t j = 0; j < _rows; ++j)
                        (*this)(j, i) = ((j == i) ? 1 : 0);
            }

            Matrix submatrix(size_t excl_col, size_t excl_row) const
            {
                Matrix result(_cols-1, _rows - 1);
                size_t res_col = 0, res_row = 0;

                for (size_t i = 0; i < _rows; ++i)
                {
                    if (i == excl_row) continue;
                    
                    res_col = 0;
                    for (size_t j = 0; j < _cols; ++j)
                        if (j != excl_col) result(res_col++, res_row) = (*this)(j, i);
                        
                    ++res_row;
                }
                
                return result;
            }

            /*
             * Returns Cofactor Matrix (square matrix only)
             */
            Matrix cofactor() const
            {
                Matrix result(_cols, _rows);
                Elem sign;

                if (_cols == 1)
                {
                    result(0, 0) = 1;
                }
                else for (size_t row = 0; row < _rows; ++row)
                {
                    for (size_t col = 0; col < _cols; ++col)
                    {
                        sign = ((col+row) % 2 == 0 ? 1 : -1);
                        result(col, row) = sign * submatrix(col, row).determinant();
                    }
                }
                return result;
            }
            
            /*
             * Returns Adjugate Matrix (square matrix only)
             */
            Matrix adjugate() const
            {
                Matrix result(_rows, _cols);
                Elem sign;
                
                if (_cols == 1)
                {
                    result(0, 0) = 1;
                }
                else for (size_t row = 0; row < _rows; ++row)
                {
                    for (size_t col = 0; col < _cols; ++col)
                    {
                        sign = ((col+row) % 2 == 0 ? 1 : -1);
                        result(row, col) = sign * submatrix(col, row).determinant();
                    }
                }
                return result;
            }

            /*
             * Matrix determinant (square matrix only)
             */

            Elem determinant() const
            {
                if (_cols != _rows) return 0;

                if (_cols == 0) return 1;
                if (_cols == 1) return _d[0];
                if (_cols == 2) return _d[0] * _d[3] - _d[1] * _d[2];

                Elem det = 0;
                Elem sign;
                
                for (size_t i = 0; i < _cols; ++i)
                {
                    sign = (i % 2 == 0 ? 1 : -1);
                    det += sign * (*this)(i, 0) * submatrix(i, 0).determinant();
                }

                return det;
            }

            /*
             * Returns inverse matrix (square matrix only)
             */

            Matrix inverse() const
            {
                Elem det = determinant();
                Matrix result = adjugate();
                result /= det;

                return result;
            }

            /*
             * Matrix filling:
             *   Matrix<elem> a(3,3);
             *   a <<  5,  6,  7,
             *         8,  9, 10,
             *        11, 12, 13;
             */

            struct Filler
            {
                Elem *_d;

                Filler(Elem *d, const Elem& first) : _d(d) { *_d = first; ++_d; }
                Filler& operator,(const Elem& elem) { *_d = elem; ++_d; return *this; }
            };

            Filler operator<<(const Elem& elem) { return Filler(_d, elem); }

            /*
             * Write Matrix to std::ostream
             */
            friend std::ostream& operator<<(std::ostream& strm, const Matrix& m)
            {
                strm << "Matrix[";
                for (size_t i = 0; i < m._rows; ++i)
                {
                    strm << "[";
                    for (size_t j = 0; j < m._cols; ++j)
                    {
                        if (j != 0) strm << ", ";
                        strm << m(j, i);
                    }
                    strm << "]";
                }
                strm << "]";
                return strm;
            }
    };

    template<class Elem>
    Matrix<Elem> operator*(const Matrix<Elem>& m1, const Matrix<Elem>& m2)
    {
        Matrix<Elem> result(m2.cols(), m1.rows());

        size_t row, col, i;

        for (row = 0; row < result.rows(); ++row)
        {
            for (col = 0; col < result.cols(); ++col)
            {
                result(col, row) = 0;
                for (i = 0; i < m1.cols(); ++i)
                    result(col, row) += m1(i, row) * m2(col, i);
            }
        }
        
        return result;
    }
}

#endif

#ifndef _GF_POLYNOMAL_APPROX_H
#define _GF_POLYNOMAL_APPROX_H

#include <map>
#include <vector>
#include <string>

#include "Matrix.h"

namespace Gf
{
    Matrix<long double> PolynomalApprox(
        const std::map<long double, long double>& values,
        size_t degree);
}
#endif

#include "PolynomalApprox.h"

#include <iostream>
#include <cmath>
#include <sstream>

namespace Gf
{
    Matrix<long double> PolynomalApprox(const std::map<long double, long double>& values, size_t degree)
    {
        typedef Matrix<long double> MT;
        
        size_t i, j;
        MT y(1, values.size());
        MT A(degree + 1, values.size());

        i = 0;
        for (auto it = values.begin(); it != values.end(); ++it)
        {
            y(0, i) = it->second;
            for (j = 0; j < degree+1; ++j)
            {
                if (j == 0) A(j, i) = 1;
                else if (j == 1) A(j, i) = it->first;
                else A(j, i) = std::pow(it->first, (long double)j);
            }
            ++i;
        }

        MT AT = A.transpose();
        MT coefs = ((AT*A).inverse() * AT) * y;

        //std::cout << "coefs:" << std::endl << coefs;

        return coefs;
    }
}