linear-algebra-latex.cpp

#include <bits/stdc++.h>

using namespace std;

#include <iostream>
#include <vector>

using std::vector;

template <typename T>
string to_string(T t) {
    string res;
    bool f = t < 0;
    if (f) {
        t = -t;
    }
    while (t != 0) {
        res += static_cast<char>('0' + t % 10);
        t /= 10;
    }
    if (res.empty()) {
        res = "0";
    }
    if (f) res += "-";
    return string(res.rbegin(), res.rend());
}

template <typename T>
T myGcd(T a, T b) {
    while (b != T()) {
        a %= b;

        T tmp = a;
        a = b;
        b = tmp;
    }
    return a;
}

template <typename T>
class Polynomial {
    vector<T> data;

 public:
    void fix() {
        while (!data.empty() && data.back() == T()) {
            data.pop_back();
        }
    }

    void fix(size_t sz) {
        if (data.size() < sz) {
            data.resize(sz);
        }
    }

    size_t size() const {
        return data.size();
    }

    Polynomial(const T& t): data({t}) {
        fix();
    }

    Polynomial(const vector<T>& v): data(v) {
        fix();
    }

    Polynomial(): data() {}

    template <typename I>
    Polynomial(I first, I last): data(first, last) {
        fix();
    }

    // equality

    bool operator == (const Polynomial<T>& to) const {
        return this->data == to.data;
    }

    bool operator != (const Polynomial<T>& to) const {
        return !(*this == to);
    }

    // calculation

    T operator[] (size_t pos) const {
        if (pos >= size()) {
            return T();
        }

        return data[pos];
    }

    int Degree() const {
        return static_cast<int>(size()) - 1;
    }

    T operator() (T t) const {
        T cur(1);
        T res = T();

        for (const auto& it : data) {
            res += cur * it;
            cur *= t;
        }

        return res;
    }

    // iterators

    typename vector<T>::const_iterator begin() const {
        return data.begin();
    }

    typename vector<T>::const_iterator end() const {
        return data.end();
    }

    // plus

    Polynomial<T>& operator += (const Polynomial<T>& to) {
        fix(to.size());
        for (size_t i = 0; i != to.size(); ++i) {
            data[i] += to[i];
        }
        fix();
        return *this;
    }

    Polynomial<T> operator + (const Polynomial<T>& to) const {
        Polynomial<T> tmp = *this;
        tmp += to;
        return tmp;
    }

    // minus

    Polynomial<T>& operator -= (const Polynomial<T>& to) {
        fix(to.size());
        for (size_t i = 0; i != to.size(); ++i) {
            data[i] -= to[i];
        }
        fix();
        return *this;
    }

    Polynomial<T> operator - (const Polynomial<T>& to) const {
        Polynomial<T> tmp = *this;
        tmp -= to;
        return tmp;
    }

    // multiplication

    Polynomial<T> operator * (const Polynomial<T>& to) const {
        vector<T> res(size() + to.size());

        for (size_t i = 0; i != size(); ++i) {
            for (size_t j = 0; j != to.size(); ++j) {
                res[i + j] += data[i] * to[j];
            }
        }

        return Polynomial<T>(std::move(res));
    }

    Polynomial<T>& operator *= (const Polynomial<T>& to) {
        return *this = *this * to;
    }

    // composition

    Polynomial<T> operator & (const Polynomial<T>& to) const {
        Polynomial<T> res;
        Polynomial<T> cur = T(1);
        for (size_t i = 0; i < size(); ++i) {
            res += cur * data[i];
            cur *= to;
        }
        return res;
    }

    // shift

    Polynomial<T> operator << (size_t cnt) const {
        vector<T> tmp(cnt);
        tmp.insert(tmp.end(), begin(), end());
        return Polynomial<T>(tmp);
    }

    // division

    std::pair<Polynomial<T>, Polynomial<T>> divide(const Polynomial<T>& to) const {
        Polynomial<T> cur = *this;
        size_t sz = cur.size();
        vector<T> res(sz);
        for (size_t i1 = 0; i1 < sz; ++i1) {
            size_t i = sz - i1 - 1;
            if (i + 1 < to.size()) {
                break;
            }
            if (i >= cur.size()) {
                continue;
            }
            T cnt = cur[i] / to.data.back();
            size_t j = i + 1 - to.size();
            res[j] += cnt;
            cur -= (to << j) * cnt;
        }
        cur.fix();
        return std::make_pair(Polynomial<T>(res), cur);
    }

    Polynomial<T> operator / (const Polynomial<T>& to) const {
        return this->divide(to).first;
    }

    Polynomial<T> operator /= (const Polynomial<T>& to) {
        return *this = *this / to;
    }

    Polynomial<T> operator % (const Polynomial<T>& to) const {
        return this->divide(to).second;
    }

    Polynomial<T> operator %= (const Polynomial<T>& to) {
        return *this = *this % to;
    }

    Polynomial<T> operator , (const Polynomial<T>& b) const {
        Polynomial<T> tmp = myGcd(*this, b);
        if (tmp.size() != 0) {
            tmp /= tmp.data.back();
        }
        return tmp;
    }
};

template <typename T>
bool operator == (const T& a, const Polynomial<T>& b) {
    return Polynomial<T>(a) == b;
}

template <typename T>
bool operator != (const T& a, const Polynomial<T>& b) {
    return Polynomial<T>(a) != b;
}

template <typename T>
Polynomial<T> operator + (const T& a, const Polynomial<T>& b) {
    return Polynomial<T>(a) + b;
}

template <typename T>
Polynomial<T> operator - (const T& a, const Polynomial<T>& b) {
    return Polynomial<T>(a) - b;
}

template <typename T>
Polynomial<T> operator * (const T& a, const Polynomial<T>& b) {
    return Polynomial<T>(a) * b;
}

template <typename T>
Polynomial<T> operator / (const T& a, const Polynomial<T>& b) {
    return Polynomial<T>(a) / b;
}

template <typename T>
Polynomial<T> operator , (const T& a, const Polynomial<T>& b) {
    return Polynomial<T>(a) , b;
}

template <typename T>
string latex(const Polynomial<T>& p) {
    string res;
    bool printed = false;
    for (size_t i1 = 0; i1 != p.size(); ++i1) {
        size_t i = p.size() - i1 - 1;
        const T& t = p[i];
        if (t == T()) {
            continue;
        }

        if (t > T() && printed) {
            res += "+";
        }
        if ((t != T(1) && t != T(-1)) || i == 0) {
            res += to_string(t);
            if (i != 0) {
                res += "";
            }
        } else {
            if (t == T(-1)) {
                res += "-";
            }
        }

        if (i != 0) {
            res += "x";
            if (i != 1) {
                res += "^";
                res += "{" + to_string(i) + "}";
            }
        }

        printed = true;
    }

    if (!printed) {
        res += "0";
    }

    return res;
}
/*
namespace static_if_detail {

    struct identity {
        template<typename T>
        T operator()(T&& x) const {
            return std::forward<T>(x);
        }
    };

    template<bool Cond>
    struct statement {
        template<typename F>
        void then(const F& f){
            f(identity());
        }
    
        template<typename F>
        void else_(const F&){}
    };
    
    template<>
    struct statement<false> {
        template<typename F>
        void then(const F&){}
    
        template<typename F>
        void else_(const F& f){
            f(identity());
        }
    };
    
} //end of namespace static_if_detail

template<bool Cond, typename F>
static_if_detail::statement<Cond> static_if(F const& f){
    static_if_detail::statement<Cond> if_;
    if_.then(f);
    return if_;
}

*/

template <typename T, size_t n, size_t m>
class Matrix {
    array<array<T, m>, n> arr;

 public:
    Matrix(initializer_list<initializer_list<T>> l) {
        auto it0 = l.begin();
        for (size_t i = 0; i < n; ++i) {
            assert(i < l.size());
            auto it1 = it0->begin();
            for (size_t j = 0; j < m; ++j) {
                assert(j < it0->size());
                arr[i][j] = *it1;
                ++it1;
            }
            ++it0;
        }
    }

    Matrix(initializer_list<T> l) {
        auto it = l.begin();
        for (size_t i = 0; i < n; ++i) {
            for (size_t j = 0; j < m; ++j) {
                assert(it != l.end());
                arr[i][j] = *(it++);
            }
        }
    }

    Matrix() {}

    const array<T, m>& operator[] (size_t i) const {
        return arr[i];
    }

    array<T, m>& operator[] (size_t i) {
        return arr[i];
    }
};

template <typename T, size_t n, size_t m>
Matrix<T, n, m> operator + (const Matrix<T, n, m>& a, const Matrix<T, n, m>& b) {
    Matrix<T, n, m> mt;
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < m; ++j) {
            mt[i][j] = a[i][j] + b[i][j];
        }
    }
    return mt;
}

template <typename T, size_t n, size_t m>
Matrix<T, n, m> operator - (const Matrix<T, n, m>& a, const Matrix<T, n, m>& b) {
    Matrix<T, n, m> mt;
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < m; ++j) {
            mt[i][j] = a[i][j] - b[i][j];
        }
    }
    return mt;
}

template <typename T, size_t n, size_t m>
bool operator == (const Matrix<T, n, m>& a, const Matrix<T, n, m>& b) {
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < m; ++j) {
            if (a[i][j] != b[i][j]) {
                return false;
            }
        }
    }
    return true;
}

template <typename T, size_t n, size_t m, size_t z>
Matrix<T, n, z> operator * (const Matrix<T, n, m>& a, const Matrix<T, m, z>& b) {
    Matrix<T, n, z> mt;
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < z; ++j) {
            mt[i][j] = 0;
            for (size_t k = 0; k < m; ++k) {
                mt[i][j] += a[i][k] * b[k][j];
            }
        }
    }
    return mt;
}

using M22 = Matrix<int, 2, 2>;
using M33 = Matrix<int, 3, 3>;
using M77 = Matrix<Polynomial<int>, 7, 7>;

template <typename T>
class Equation {
public:
    T t;
    Equation(T t): t(t) {}
};

template <typename A, size_t n, size_t m>
Equation<Matrix<A, n, m>> operator !(Matrix<A, n, m> mt) {
    return Equation<Matrix<A, n, m>>(mt);
}

template <typename A, typename B>
class Op {
public:
    A a;
    B b;
    char op;
    Op(A a, B b, char op): a(a), b(b), op(op) {}
};

template <typename A, typename B>
Op<A, B> operator + (A a, B b) {
    return Op<A, B>(a, b, '+');
}

template <typename A, typename B>
Op<A, B> operator - (A a, B b) {
    return Op<A, B>(a, b, '-');
}

template <typename A, typename B>
Op<A, B> operator * (A a, B b) {
    return Op<A, B>(a, b, '*');
}

template <typename T>
class Transponate {
public:
    T t;
    Transponate(T t): t(t) {}
};

template <typename T>
Transponate<T> tt(T t) {
    return Transponate<T>(t);
}

template <typename T>
class Square {
public:
    T t;
    Square(T t): t(t) {}
};

template <typename T>
Square<T> sq(T t) {
    return Square<T>(t);
}

template <typename T>
class Gather {
public:
    T t;
    Gather(T t): t(t) {}
};

template <typename T>
Gather<T> gather(T t) {
    return Gather<T>(t);
}

////////////

template <typename T>
string latexC(T t);

string latex(int i) {
    return to_string(i);
}

template <typename T, size_t n>
string latex(const array<T, n>& arr) {
    string res;
    for (size_t i = 0; i < n; ++i) {
        if (i != 0) {
            res += "& ";
        }
        res += latex(arr[i]);
    }
    return res;
}

template <typename T, size_t n, size_t m>
string latex(const Matrix<T, n, m>& mt) {
    string res = "\\begin{bmatrix}\n";
    for (size_t i = 0; i < n; ++i) {
        res += latex(mt[i]) + "\\\\\n";
    }
    res +=  "\\end{bmatrix}";
    return res;
}

template <typename T>
string latex(const Equation<T>& t) {
    return latex(t.t);
}

string latex(char c) {
    if (c == '*') return "\\cdot";
    return string(1, c);
}

template <typename A, typename B, typename C>
string latex(const Op<Op<A, B>, C>& p) {
    string l, r;
    if (p.op == '*' && p.a.op != '*') {
        l = latexC(p.a);
        r = latexC(p.b);
    } else {
        l = latex(p.a);
        r = latex(p.b);
    }
    return l + " " + latex(p.op) + " " + r;
}

template <typename A, typename B>
string latex(const Op<A, B>& p) {
    string l, r;
    if (p.op == '*') {
        l = latexC(p.a);
        r = latexC(p.b);
    } else {
        l = latex(p.a);
        r = latex(p.b);
    }
    return l + " " + latex(p.op) + " " + r;
}

template <typename T>
string latex(const Transponate<T>& t) {
    return latexC(t.t) + "^{\\mathrm{T}}";
}

template <typename T>
string latex(const Square<T>& t) {
    return latexC(t.t) + "^{2}";
}

template <typename T>
string latex(const Gather<T>& t) {
    return "\\begin{gather*}" + latex(t.t) + "\\end{gather*}";
}

////////////

template <typename T, size_t n, size_t m>
M22 calc(const Matrix<T, n, m>& mt) {
    return mt;
}

template <typename T>
M22 calc(const Equation<T>& t) {
    return calc(t.t);
}

template <typename A, typename B>
M22 calc(const Op<A, B>& p) {
    if (p.op == '+') {
        return calc(p.a) + calc(p.b);
    }
    if (p.op == '-') {
        return calc(p.a) - calc(p.b);
    }
    if (p.op == '*') {
        return calc(p.a) * calc(p.b);
    }
    assert(false);
}

template <typename T>
M22 calc(const Transponate<T>& t) {
    M22 mt = calc(t.t);
    swap(mt[0][1], mt[1][0]);
    return mt;
}

template <typename T>
M22 calc(const Square<T>& t) {
    return calc(t.t) * calc(t.t);
}

////////////

template <typename A, typename B>
string latexC(Op<A, B> t) {
    return "{\\left(" + latex(t) + "\\right)}";
}

template <typename T>
string latexC(T t) {
    return latex(t);
}

//////////

int f2(int x) {
    if (x == 0) return 0;
    if (x == 1) return 1;
    return f2(x - 2) * 2 + f2(x - 1);
}

int f2M(int n) {
    return (-pow(-1, n) + pow(2, n)) / 3;
}

int sgn(vector<int> v) {
    int t = 0;
    for (int i = 0; i < v.size(); i++) {
        for (int j = i + 1; j < v.size(); j++) {
            t += v[j] < v[i];
        }
    }
    if (t % 2 == 0) {
        return 1;
    }
    return -1;
}

int main() {
    int task = 4;
    if (task == 1) {
        auto a = !M22{3, 11, 4, 4};
        auto b = !M22{3, 8, 3, 4};
        auto c = !M22{3, 3, 8, 4};
        auto d = !M22{3, 10, 3, 5};
        auto e = !M22{6, 21, 7, 9};

        auto eq = a * b * c + d * tt(d * b) 
        - e * sq(c) 
        + a * tt(d * b)
         + e * tt(a * b) + d * b * b;

        auto res = calc(eq);

        cout << latex(gather(eq)) << endl;
    
        auto eq0 = a * b * c + d * tt(d * b) 
        + a * tt(d * b)
        - e * sq(c) 
         + e * tt(a * b) + d * b * b;

        assert(calc(eq0) == res);

        cout << latex(gather(eq0)) << endl;

        cout << "\nВынесем $" << latex(tt(d * b)) << "$ за скобки.\n";

        auto eq1 = a * b * c 
        + (d + a) * tt(d * b) 
        - e * sq(c) 
        + e * tt(a * b) 
        + d * b * b;

        assert(calc(eq1) == res);

        cout << latex(gather(eq1)) << "\n\n";

        cout << "Аналогичным образом вынесем $" << latex(c) << "$ за скобки." << endl;

        auto eq2 = (a * b - e * c) * c
        + (d + a) * tt(d * b) 
        + e * tt(a * b) 
        + d * b * b;

        assert(calc(eq2) == res);

        cout << latex(gather(eq2)) << endl;

        auto eq3 = (!calc(a * b) - !calc(e * c)) * c
        + !calc(d + a) * tt(!calc(d * b)) 
        + e * tt(!calc(a * b)) 
        + d * !calc(b * b);

        assert(calc(eq3) == res);

        cout << latex(gather(eq3)) << endl;

        auto eq4 = (!calc(!calc(a * b) - !calc(e * c))) * c
        + !calc(d + a) * !calc(tt(!calc(d * b)))
        + e * !calc(tt(!calc(a * b))) 
        + d * !calc(b * b);

        assert(calc(eq4) == res);

        cout << latex(gather(eq4)) << endl;

        auto eq5 = (!calc(!calc(a * b) - !calc(e * c))) * c
        + e * (
            !calc(tt(!calc(d * b)))
            + !calc(tt(!calc(a * b))) 
        )
        + d * !calc(b * b);

        assert(calc(eq5) == res);

        cout << "Вынесем $" << latex(e) << "$ за скобки" << endl;

        cout << latex(gather(eq5)) << endl;

        auto eq6 = (!calc(!calc(a * b) - !calc(e * c))) * c
        + e * !calc(
            !calc(tt(!calc(d * b)))
            + !calc(tt(!calc(a * b))) 
        )
        + d * !calc(b * b);

        assert(calc(eq6) == res);

        cout << latex(gather(eq6)) << endl;

        auto eq7 = !calc((!calc(!calc(a * b) - !calc(e * c))) * c)
        + !calc(e * !calc(
            !calc(tt(!calc(d * b)))
            + !calc(tt(!calc(a * b))) 
        ))
        + !calc(d * !calc(b * b));

        assert(calc(eq7) == res);

        cout << latex(gather(eq7)) << endl;

        cout << latex(gather(res)) << endl;
    }
    if (task == 2) {

        for (size_t i = 0; i < 15; ++i) {
            cout << f2M(i) << " ";
        }
        cout << endl;

        auto x = Polynomial<int>(vector<int>{0, 1});
        //Matrix<Polynomial<int>, 3, 3> A{0, x, 0, x, x, x, 0, x, 0};
        M33 A{0, 4, 0, 4, 4, 4, 0, 4, 0};
        auto cur = A;
        for (size_t i = 1; i <= 10; ++i) {
            cout << "A^{" << i << "} = " << latex(cur) << ", ";
            cur = cur * A;
        }
    }
    if (task == 4) {
        auto x = Polynomial<int>(vector<int>{0, 1});
        M77 A{
            {4, x, 10, 7, 1, 0, 9},
            {x, 8, 1, 9, 3, 2, x},
            {10, 1, 5, 2, 7, x, 4},
            {7, 9, 2, 10, x, 3, 1},
            {1, 3, 7, x, 9, 8, 6},
            {0, 2, x, 3, 8, 7, 5},
            {9, x, 4, 1, 6, 5, 3}
        };

        auto det = x - x;

        int n = 7;

        vector<int> v(n);
        for (int i = 0; i < n; i++) v[i] = i;

        do {
            int f = sgn(v);
            Polynomial<int> cur = f;
            for (int i = 0; i < n; i++) {
                cur *= A[i][v[i]];
            }
            det += cur;
        } while (next_permutation(v.begin(), v.end()));

        cout << latex(det) << endl;
        cout << latex(A) << endl;
    }
}