#include <vector>
#include <iostream>
#include <algorithm>
#include <functional>
#include <numeric>

using namespace std;

long gcd(long a, long b) {
    if (a < 0 || b < 0) return gcd(abs(a), abs(b));
    else if (min(a, b) == 0)
        return max(1l, max(a, b));
    else if (a == b)
        return a;
    else
        return gcd(min(a, b), max(a, b) - min(a, b));
}

struct vec {
    vec(vector<long> v) : v_(v) {}
    vector<long>& v() { return v_; }
    const vector<long>& v() const { return v_; }

    vec operator+(const vec&) const;
    vec operator*(long) const;
    long operator*(const vec&) const;

    void reduceByGCD();
    
    bool isZero() const { return v() == vector<long>(v().size()); }

private:
    vector<long> v_;
    long gcd() const;
};

long
vec::gcd() const
{
    if (v_.empty()) return 1;
    long res = v_[0];
    for (int i = 1; i < v_.size(); ++i)
        res = ::gcd(res, v_[i]);
    return res;
}

void
vec::reduceByGCD()
{
    long div = gcd();
    for (auto &a: v_) a /= div;
}

vec
vec::operator+(const vec& r) const
{
    vector<long> res(v().size());
    transform(this->v_.begin(), this->v_.end(), r.v_.begin(), res.begin(), plus<long>());
    return vec(res);
}

vec
vec::operator*(long k) const
{
    vector<long> res(v().size());
    transform(this->v_.begin(), this->v_.end(), res.begin(), bind(multiplies<long>(), placeholders::_1, k));
    return vec(res);
}

long
vec::operator*(const vec& r) const
{
    return inner_product(v_.begin(), v_.end(), r.v().begin(), 0l);
}

vec gschmidt_and_reduce(const vec& what, const vec& other)
{
    long sqNorm = other*other;
    vec res = what*sqNorm  + other*(what*other)*-1;
    res.reduceByGCD();
    return res;
}

vector<vec>
gschmidt_and_reduce(const vector<vec>& vs)
{
    vector<vec> res;
    for (size_t i = 0; i < vs.size(); ++i) {
        vec cur = vs[i];
        for (int k = 0; k < i; ++k)
            cur = gschmidt_and_reduce(cur, res[k]);
        res.push_back(cur);
    }
    return res;
}

int main() {
    vector<vec> A = {
        vec({ 3, 5, 4, -4 }),
        vec({ -2, -3, -2, 1 }),
        vec({ -5, -9, -6, 8 }),
        vec({ 1, 0, 0, 0 }),
        vec({ 0, 1, 0, 0 }),
        vec({ 0, 0, 1, 0 }),
        vec({ 0, 0, 0, 1 }),
    };
    vector<vec> orthogonilizedA = gschmidt_and_reduce(A);
    cout << "Orthogonalized a:"<< endl;
    for (int i = 0; i < orthogonilizedA.size(); ++i) {
        if (!orthogonilizedA[i].isZero()) {
            for (auto coordinate: orthogonilizedA[i].v())
                cout << coordinate << ", ";
            cout << endl;
        }
        if (i == 2)
            cout << "Orthogonal complement of a:" << endl;
    }

    vector<vec> B = {
        vec({ -4, -5, -3, 3 }),
        vec({ 2, 7, 0, -6 }),
        vec({ -1, 0, -1, 0 }),
        vec({ 1, 0, 0, 0 }),
        vec({ 0, 1, 0, 0 }),
        vec({ 0, 0, 1, 0 }),
        vec({ 0, 0, 0, 1 }),
    };
    vector<vec> orthogonilizedB = gschmidt_and_reduce(B);
    cout << "Orthogonalized b:"<< endl;
    for (int i = 0; i < orthogonilizedB.size(); ++i) {
        if (!orthogonilizedB[i].isZero()) {
            for (auto coordinate: orthogonilizedB[i].v())
                cout << coordinate << ", ";
            cout << endl;
        }
        if (i == 2)
            cout << "Orthogonal complement of b:" << endl;
    }
    return 0;
}