CP-Algorithms Library

This documentation is automatically generated by competitive-verifier/competitive-verifier

View the Project on GitHub cp-algorithms/cp-algorithms-aux

:heavy_check_mark: cp-algo/linalg/vector.hpp

Depends on

Required by

Verified with

Code

#ifndef CP_ALGO_LINALG_VECTOR_HPP
#define CP_ALGO_LINALG_VECTOR_HPP
#include "../random/rng.hpp"
#include "../math/modint.hpp"
#include <functional>
#include <algorithm>
#include <valarray>
#include <iostream>
#include <iterator>
#include <cassert>
namespace cp_algo::math::linalg {
    template<class vec, typename base>
    struct valarray_base: std::valarray<base> {
        using Base = std::valarray<base>;
        using Base::Base;

        valarray_base(base const& t): Base(t, 1) {}

        auto begin() {return std::begin(to_valarray());}
        auto begin() const {return std::begin(to_valarray());}
        auto end() {return std::end(to_valarray());}
        auto end() const {return std::end(to_valarray());}

        bool operator == (vec const& t) const {return std::ranges::equal(*this, t);}
        bool operator != (vec const& t) const {return !(*this == t);}

        vec operator-() const {return Base::operator-();}

        static vec from_range(auto const& R) {
            vec res(std::ranges::distance(R));
            std::ranges::copy(R, res.begin());
            return res;
        }
        Base& to_valarray() {return static_cast<Base&>(*this);}
        Base const& to_valarray() const {return static_cast<Base const&>(*this);}
    };

    template<class vec, typename base>
    vec operator+(valarray_base<vec, base> const& a, valarray_base<vec, base> const& b) {
        return a.to_valarray() + b.to_valarray();
    }
    template<class vec, typename base>
    vec operator-(valarray_base<vec, base> const& a, valarray_base<vec, base> const& b) {
        return a.to_valarray() - b.to_valarray();
    }

    template<class vec, typename base>
    struct vec_base: valarray_base<vec, base> {
        using Base = valarray_base<vec, base>;
        using Base::Base;

        static vec ei(size_t n, size_t i) {
            vec res(n);
            res[i] = 1;
            return res;
        }

        virtual void add_scaled(vec const& b, base scale, size_t i = 0) {
            for(; i < size(*this); i++) {
                (*this)[i] += scale * b[i];
            }
        }
        virtual vec const& normalize() {
            return static_cast<vec&>(*this);
        }
        virtual base normalize(size_t i) {
            return (*this)[i];
        }
        void read() {
            for(auto &it: *this) {
                std::cin >> it;
            }
        }
        void print() const {
            std::ranges::copy(*this, std::ostream_iterator<base>(std::cout, " "));
            std::cout << "\n";
        }
        static vec random(size_t n) {
            vec res(n);
            std::ranges::generate(res, random::rng);
            return res;
        }
        // Concatenate vectors
        vec operator |(vec const& t) const {
            vec res(size(*this) + size(t));
            res[std::slice(0, size(*this), 1)] = *this;
            res[std::slice(size(*this), size(t), 1)] = t;
            return res;
        }

        // Generally, vec shouldn't be modified
        // after its pivot index is set
        std::pair<size_t, base> find_pivot() {
            if(pivot == size_t(-1)) {
                pivot = 0;
                while(pivot < size(*this) && normalize(pivot) == base(0)) {
                    pivot++;
                }
                if(pivot < size(*this)) {
                    pivot_inv = base(1) / (*this)[pivot];
                }
            }
            return {pivot, pivot_inv};
        }
        void reduce_by(vec &t) {
            auto [pivot, pinv] = t.find_pivot();
            if(pivot < size(*this)) {
                add_scaled(t, -normalize(pivot) * pinv, pivot);
            }
        }
    private:
        size_t pivot = -1;
        base pivot_inv;
    };

    template<typename base>
    struct vec: vec_base<vec<base>, base> {
        using Base = vec_base<vec<base>, base>;
        using Base::Base;
    };

    template<modint_type base>
    struct vec<base>: vec_base<vec<base>, base> {
        using Base = vec_base<vec<base>, base>;
        using Base::Base;

        void add_scaled(vec const& b, base scale, size_t i = 0) override {
            for(; i < size(*this); i++) {
                (*this)[i].add_unsafe(scale.getr() * b[i].getr());
            }
            if(++counter == 8) {
                for(auto &it: *this) {
                    it.pseudonormalize();
                }
                counter = 0;
            }
        }
        vec const& normalize() override {
            for(auto &it: *this) {
                it.normalize();
            }
            return *this;
        }
        base normalize(size_t i) override {
            return (*this)[i].normalize();
        }
    private:
        size_t counter = 0;
    };
}
#endif // CP_ALGO_LINALG_VECTOR_HPP

#line 1 "cp-algo/linalg/vector.hpp"


#line 1 "cp-algo/random/rng.hpp"


#include <chrono>
#include <random>
namespace cp_algo::random {
    uint64_t rng() {
        static std::mt19937_64 rng(
            std::chrono::steady_clock::now().time_since_epoch().count()
        );
        return rng();
    }
}

#line 1 "cp-algo/math/modint.hpp"


#line 1 "cp-algo/math/common.hpp"


#include <functional>
#include <cstdint>
namespace cp_algo::math {
#ifdef CP_ALGO_MAXN
    const int maxn = CP_ALGO_MAXN;
#else
    const int maxn = 1 << 19;
#endif
    const int magic = 128; // threshold for sizes to run the naive algo

    auto bpow(auto const& x, int64_t n, auto const& one, auto op) {
        if(n == 0) {
            return one;
        } else {
            auto t = bpow(x, n / 2, one, op);
            t = op(t, t);
            if(n % 2) {
                t = op(t, x);
            }
            return t;
        }
    }
    auto bpow(auto x, int64_t n, auto ans) {
        return bpow(x, n, ans, std::multiplies{});
    }
    template<typename T>
    T bpow(T const& x, int64_t n) {
        return bpow(x, n, T(1));
    }
}

#line 4 "cp-algo/math/modint.hpp"
#include <iostream>
namespace cp_algo::math {
    template<typename modint>
    struct modint_base {
        static int64_t mod() {
            return modint::mod();
        }
        modint_base(): r(0) {}
        modint_base(int64_t rr): r(rr % mod()) {
            r = std::min(r, r + mod());
        }
        modint inv() const {
            return bpow(to_modint(), mod() - 2);
        }
        modint operator - () const {return std::min(-r, mod() - r);}
        modint& operator /= (const modint &t) {
            return to_modint() *= t.inv();
        }
        modint& operator *= (const modint &t) {
            if(mod() <= uint32_t(-1)) {
                r = r * t.r % mod();
            } else {
                r = __int128(r) * t.r % mod();
            }
            return to_modint();
        }
        modint& operator += (const modint &t) {
            r += t.r; r = std::min(r, r - mod());
            return to_modint();
        }
        modint& operator -= (const modint &t) {
            r -= t.r; r = std::min(r, r + mod());
            return to_modint();
        }
        modint operator + (const modint &t) const {return modint(to_modint()) += t;}
        modint operator - (const modint &t) const {return modint(to_modint()) -= t;}
        modint operator * (const modint &t) const {return modint(to_modint()) *= t;}
        modint operator / (const modint &t) const {return modint(to_modint()) /= t;}
        auto operator <=> (const modint_base &t) const = default;
        int64_t rem() const {return 2 * r > (uint64_t)mod() ? r - mod() : r;}

        // Only use if you really know what you're doing!
        uint64_t modmod() const {return 8ULL * mod() * mod();};
        void add_unsafe(uint64_t t) {r += t;}
        void pseudonormalize() {r = std::min(r, r - modmod());}
        modint const& normalize() {
            if(r >= (uint64_t)mod()) {
                r %= mod();
            }
            return to_modint();
        }
        uint64_t& setr() {return r;}
        uint64_t getr() const {return r;}
    private:
        uint64_t r;
        modint& to_modint() {return static_cast<modint&>(*this);}
        modint const& to_modint() const {return static_cast<modint const&>(*this);}
    };
    template<typename modint>
    std::istream& operator >> (std::istream &in, modint_base<modint> &x) {
        return in >> x.setr();
    }
    template<typename modint>
    std::ostream& operator << (std::ostream &out, modint_base<modint> const& x) {
        return out << x.getr();
    }

    template<typename modint>
    concept modint_type = std::is_base_of_v<modint_base<modint>, modint>;

    template<int64_t m>
    struct modint: modint_base<modint<m>> {
        static constexpr int64_t mod() {return m;}
        using Base = modint_base<modint<m>>;
        using Base::Base;
    };

    struct dynamic_modint: modint_base<dynamic_modint> {
        static int64_t mod() {return m;}
        static void switch_mod(int64_t nm) {m = nm;}
        using Base = modint_base<dynamic_modint>;
        using Base::Base;

        // Wrapper for temp switching
        auto static with_mod(int64_t tmp, auto callback) {
            struct scoped {
                int64_t prev = mod();
                ~scoped() {switch_mod(prev);}
            } _;
            switch_mod(tmp);
            return callback();
        }
    private:
        static int64_t m;
    };
    int64_t dynamic_modint::m = 0;
}

#line 6 "cp-algo/linalg/vector.hpp"
#include <algorithm>
#include <valarray>
#line 9 "cp-algo/linalg/vector.hpp"
#include <iterator>
#include <cassert>
namespace cp_algo::math::linalg {
    template<class vec, typename base>
    struct valarray_base: std::valarray<base> {
        using Base = std::valarray<base>;
        using Base::Base;

        valarray_base(base const& t): Base(t, 1) {}

        auto begin() {return std::begin(to_valarray());}
        auto begin() const {return std::begin(to_valarray());}
        auto end() {return std::end(to_valarray());}
        auto end() const {return std::end(to_valarray());}

        bool operator == (vec const& t) const {return std::ranges::equal(*this, t);}
        bool operator != (vec const& t) const {return !(*this == t);}

        vec operator-() const {return Base::operator-();}

        static vec from_range(auto const& R) {
            vec res(std::ranges::distance(R));
            std::ranges::copy(R, res.begin());
            return res;
        }
        Base& to_valarray() {return static_cast<Base&>(*this);}
        Base const& to_valarray() const {return static_cast<Base const&>(*this);}
    };

    template<class vec, typename base>
    vec operator+(valarray_base<vec, base> const& a, valarray_base<vec, base> const& b) {
        return a.to_valarray() + b.to_valarray();
    }
    template<class vec, typename base>
    vec operator-(valarray_base<vec, base> const& a, valarray_base<vec, base> const& b) {
        return a.to_valarray() - b.to_valarray();
    }

    template<class vec, typename base>
    struct vec_base: valarray_base<vec, base> {
        using Base = valarray_base<vec, base>;
        using Base::Base;

        static vec ei(size_t n, size_t i) {
            vec res(n);
            res[i] = 1;
            return res;
        }

        virtual void add_scaled(vec const& b, base scale, size_t i = 0) {
            for(; i < size(*this); i++) {
                (*this)[i] += scale * b[i];
            }
        }
        virtual vec const& normalize() {
            return static_cast<vec&>(*this);
        }
        virtual base normalize(size_t i) {
            return (*this)[i];
        }
        void read() {
            for(auto &it: *this) {
                std::cin >> it;
            }
        }
        void print() const {
            std::ranges::copy(*this, std::ostream_iterator<base>(std::cout, " "));
            std::cout << "\n";
        }
        static vec random(size_t n) {
            vec res(n);
            std::ranges::generate(res, random::rng);
            return res;
        }
        // Concatenate vectors
        vec operator |(vec const& t) const {
            vec res(size(*this) + size(t));
            res[std::slice(0, size(*this), 1)] = *this;
            res[std::slice(size(*this), size(t), 1)] = t;
            return res;
        }

        // Generally, vec shouldn't be modified
        // after its pivot index is set
        std::pair<size_t, base> find_pivot() {
            if(pivot == size_t(-1)) {
                pivot = 0;
                while(pivot < size(*this) && normalize(pivot) == base(0)) {
                    pivot++;
                }
                if(pivot < size(*this)) {
                    pivot_inv = base(1) / (*this)[pivot];
                }
            }
            return {pivot, pivot_inv};
        }
        void reduce_by(vec &t) {
            auto [pivot, pinv] = t.find_pivot();
            if(pivot < size(*this)) {
                add_scaled(t, -normalize(pivot) * pinv, pivot);
            }
        }
    private:
        size_t pivot = -1;
        base pivot_inv;
    };

    template<typename base>
    struct vec: vec_base<vec<base>, base> {
        using Base = vec_base<vec<base>, base>;
        using Base::Base;
    };

    template<modint_type base>
    struct vec<base>: vec_base<vec<base>, base> {
        using Base = vec_base<vec<base>, base>;
        using Base::Base;

        void add_scaled(vec const& b, base scale, size_t i = 0) override {
            for(; i < size(*this); i++) {
                (*this)[i].add_unsafe(scale.getr() * b[i].getr());
            }
            if(++counter == 8) {
                for(auto &it: *this) {
                    it.pseudonormalize();
                }
                counter = 0;
            }
        }
        vec const& normalize() override {
            for(auto &it: *this) {
                it.normalize();
            }
            return *this;
        }
        base normalize(size_t i) override {
            return (*this)[i].normalize();
        }
    private:
        size_t counter = 0;
    };
}
Back to top page