MRSparsePolynomial.h

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
#include <cassert>
#include <utility>
#include <vector>
 
namespace MR
{
 
 
template <typename C, typename D, D M>
class SparsePolynomial;

template <typename T, typename D, D M>
using SparsePolynomialProduct = SparsePolynomial<decltype( std::declval<T>() * std::declval<T>() ), D, M>;

template <typename T, typename C, typename D, D M>
SparsePolynomialProduct<T,D,M> mulAs( const SparsePolynomial<C,D,M>& a, const SparsePolynomial<C,D,M>& b );

template <typename C, typename D, D M>
SparsePolynomial<C,D,M> operator *( const SparsePolynomial<C,D,M>& a, const SparsePolynomial<C,D,M>& b ) { return mulAs<C>( a, b ); }

template <typename T, typename C, typename D, D M>
[[nodiscard]] int signOfProductsDiff( const SparsePolynomial<C,D,M>& a, const SparsePolynomial<C,D,M>& b,
    const SparsePolynomial<C,D,M>& c, const SparsePolynomial<C,D,M>& d );

template <typename C, typename D, D M>
class SparsePolynomial
{
    static_assert( M > 0 );
public:
    using Term = std::pair<D, C>;

    SparsePolynomial() = default;

    SparsePolynomial( std::vector<Term> && );

    SparsePolynomial( C c0, D d1, C c1 );

    SparsePolynomial( C c0, D d1, C c1, D d2, C c2 );

    [[nodiscard]] static SparsePolynomial fromUnsortedTerms( std::vector<Term> && terms );

    void setZeroCoeff( D d )
    {
        auto it = std::lower_bound( terms_.begin(), terms_.end(), d,
            []( const Term & t, D d ) { return t.first < d; } );
        if ( it != terms_.end() && it->first == d )
            terms_.erase( it );
    }

    [[nodiscard]] bool empty() const { return terms_.empty(); }

    [[nodiscard]] bool isPositive() const;

    [[nodiscard]] const std::vector<Term> & get() const { return terms_; }
 
    SparsePolynomial& operator +=( const SparsePolynomial& b );
    SparsePolynomial& operator -=( const SparsePolynomial& b );
    [[nodiscard]] friend SparsePolynomial operator +( SparsePolynomial a, const SparsePolynomial& b ) { a += b; return a; }
    [[nodiscard]] friend SparsePolynomial operator -( SparsePolynomial a, const SparsePolynomial& b ) { a -= b; return a; }
    template <typename T, typename C2, typename D2, D2 M2>
    friend SparsePolynomialProduct<T,D2,M2> mulAs( const SparsePolynomial<C2,D2,M2>& a, const SparsePolynomial<C2,D2,M2>& b );
 
private:
    void mergeTerms_();
 
    std::vector<Term> terms_; 
};
 
template <typename C, typename D, D M>
SparsePolynomial<C,D,M>::SparsePolynomial( std::vector<Term> && terms ) : terms_( std::move( terms ) )
{
#ifndef NDEBUG
    for ( size_t i = 0; i < terms_.size(); ++i )
    {
        assert( terms_[i].first <= M );
        assert( terms_[i].second != 0 );
        assert( i == 0 || terms_[i - 1].first < terms_[i].first );
    }
#endif
}
 
template <typename C, typename D, D M>
SparsePolynomial<C,D,M>::SparsePolynomial( C c0, D d1, C c1 )
{
    assert( c1 != 0 );
    assert( d1 != 0 );
    if ( c0 != 0 )
        terms_.emplace_back( D(0), c0 );
    if ( d1 <= M )
        terms_.emplace_back( d1, c1 );
}
 
template <typename C, typename D, D M>
SparsePolynomial<C,D,M>::SparsePolynomial( C c0, D d1, C c1, D d2, C c2 )
{
    assert( c1 != 0 );
    assert( d1 != 0 );
    assert( c2 != 0 );
    assert( d2 != 0 );
    assert( d1 != d2 );
    if ( c0 != 0 )
        terms_.emplace_back( D(0), c0 );
    if ( d1 > d2 )
    {
        std::swap( d1, d2 );
        std::swap( c1, c2 );
    }
    if ( d1 <= M )
        terms_.emplace_back( d1, c1 );
    if ( d2 <= M )
        terms_.emplace_back( d2, c2 );
}
 
template <typename C, typename D, D M>
SparsePolynomial<C,D,M> SparsePolynomial<C,D,M>::fromUnsortedTerms( std::vector<Term> && terms )
{
    std::sort( terms.begin(), terms.end(), []( const Term & x, const Term & y ) { return x.first < y.first; } );
    assert( terms.empty() || terms.back().first <= M );
    SparsePolynomial res;
    res.terms_ = std::move( terms );
    res.mergeTerms_();
    return res;
}
 
template <typename C, typename D, D M>
bool SparsePolynomial<C,D,M>::isPositive() const
{
    if ( !terms_.empty() )
        return terms_.front().second > 0;
 
    assert (false);
    return false;
}
 
template <typename C, typename D, D M>
void SparsePolynomial<C,D,M>::mergeTerms_()
{
    size_t out = 0;
    for ( size_t i = 0; i < terms_.size(); )
    {
        auto deg = terms_[i].first;
        auto cf = std::move( terms_[i].second );
        for ( ++i; i < terms_.size() && terms_[i].first == deg; ++i )
            cf += terms_[i].second;
        if ( cf != 0 )
            terms_[out++] = { deg, std::move( cf ) };
    }
    terms_.resize( out );
}
 
template <typename C, typename D, D M>
SparsePolynomial<C,D,M>& SparsePolynomial<C,D,M>::operator +=( const SparsePolynomial& b )
{
    std::vector<Term> res;
    res.reserve( terms_.size() + b.terms_.size() );
    std::merge( std::make_move_iterator( terms_.begin() ), std::make_move_iterator( terms_.end() ),
        b.terms_.begin(), b.terms_.end(), std::back_inserter( res ),
        []( const Term & x, const Term & y ) { return x.first < y.first; } );
    terms_ = std::move( res );
    mergeTerms_();
    return * this;
}
 
template <typename C, typename D, D M>
SparsePolynomial<C,D,M>& SparsePolynomial<C,D,M>::operator -=( const SparsePolynomial& b )
{
    std::vector<Term> res;
    res.reserve( terms_.size() + b.terms_.size() );
    auto itA = terms_.begin();
    auto itB = b.terms_.begin();
    while ( itA != terms_.end() && itB != b.terms_.end() )
    {
        if ( itB->first < itA->first )
        {
            res.emplace_back( itB->first, -itB->second );
            ++itB;
        }
        else
        {
            res.push_back( std::move( *itA ) );
            ++itA;
        }
    }
    for ( ; itA != terms_.end(); ++itA )
        res.push_back( std::move( *itA ) );
    for ( ; itB != b.terms_.end(); ++itB )
        res.emplace_back( itB->first, -itB->second );
    terms_ = std::move( res );
    mergeTerms_();
    return * this;
}
 
template <typename T, typename C, typename D, D M>
[[nodiscard]] SparsePolynomialProduct<T,D,M> mulAs( const SparsePolynomial<C,D,M>& a, const SparsePolynomial<C,D,M>& b )
{
    using Res = SparsePolynomialProduct<T,D,M>;
    using Term = typename Res::Term;
    std::vector<Term> res;
    res.reserve( a.terms_.size() * b.terms_.size() );
    for ( const auto & [degA, cfA] : a.terms_ )
    {
        assert( cfA != 0 );
        for ( const auto & [degB, cfB] : b.terms_ )
        {
            assert( cfB != 0 );
            const auto deg = degA + degB;
            if ( deg > M )
                break;
            res.emplace_back( deg, T( cfA ) * T( cfB ) );
        }
    }
    return Res::fromUnsortedTerms( std::move( res ) );
}
 
template <typename T, typename C, typename D, D M>
int signOfProductsDiff( const SparsePolynomial<C,D,M>& a, const SparsePolynomial<C,D,M>& b,
    const SparsePolynomial<C,D,M>& c, const SparsePolynomial<C,D,M>& d )
{
    struct RowTerm
    {
        D deg;
        int i, j; 
        bool neg; 
    };
    auto greater = []( const RowTerm & x, const RowTerm & y ) { return x.deg > y.deg; };
    std::vector<RowTerm> heap;
    heap.reserve( a.get().size() + c.get().size() );
    if ( !b.get().empty() )
        for ( int i = 0; i < (int)a.get().size() && a.get()[i].first + b.get()[0].first <= M; ++i )
            heap.push_back( { a.get()[i].first + b.get()[0].first, i, 0, false } );
    if ( !d.get().empty() )
        for ( int i = 0; i < (int)c.get().size() && c.get()[i].first + d.get()[0].first <= M; ++i )
            heap.push_back( { c.get()[i].first + d.get()[0].first, i, 0, true } );
    std::make_heap( heap.begin(), heap.end(), greater );
 
    while ( !heap.empty() )
    {
        const auto deg = heap.front().deg;
        decltype( std::declval<T>() * std::declval<T>() ) coeff{};
        do
        {
            std::pop_heap( heap.begin(), heap.end(), greater );
            auto r = heap.back();
            heap.pop_back();
            const auto & f = ( r.neg ? c : a ).get();
            const auto & g = ( r.neg ? d : b ).get();
            const auto prod = T( f[r.i].second ) * T( g[r.j].second );
            if ( r.neg )
                coeff -= prod;
            else
                coeff += prod;
            if ( ++r.j < (int)g.size() && ( r.deg = f[r.i].first + g[r.j].first ) <= M )
            {
                heap.push_back( r );
                std::push_heap( heap.begin(), heap.end(), greater );
            }
        }
        while ( !heap.empty() && heap.front().deg == deg );
        if ( coeff > 0 )
            return 1;
        if ( coeff < 0 )
            return -1;
    }
    return 0;
}
 
}