poly.h

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
#include <limits>
 
#include "../config.h"

namespace alfi::util::poly {
    template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
    void normalize(Container<Number>& p, const Number& epsilon = std::numeric_limits<Number>::epsilon()) {
        if (p.empty()) {
            return p.push_back(0);
        }
        auto p_start = std::find_if(p.begin(), p.end(), [&epsilon](Number v) { return std::abs(v) > epsilon; });
        if (p_start == p.end()) {
            --p_start;
        }
        if (p_start > p.begin()) {
            p.erase(p.begin(), p_start);
        }
    }

    template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
    Container<Number> mul(const Container<Number>& p1, const Container<Number>& p2) {
        if (p1.empty() || p2.empty()) {
            return {};
        }
        Container<Number> result(p1.size() + p2.size() - 1);
        for (SizeT i = 0; i < p1.size(); ++i) {
            for (SizeT j = 0; j < p2.size(); ++j) {
                result[i+j] += p1[i] * p2[j];
            }
        }
        return result;
    }

    template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
    std::pair<Container<Number>,Container<Number>> div(const Container<Number>& dividend, const Container<Number>& divisor, const Number& epsilon = std::numeric_limits<Number>::epsilon()) {
        const auto divisor_start = std::find_if(divisor.begin(), divisor.end(), [&epsilon](Number v) { return std::abs(v) > epsilon; });
 
        if (divisor_start == divisor.end() || dividend.size() < divisor.size()) {
            return {{}, dividend};
        }
 
        const auto divisor_start_idx = divisor_start - divisor.begin();
 
        const auto n = dividend.size();
        const auto m = divisor.size() - divisor_start_idx;
 
        Container<Number> quotient(n - m + 1, 0);
        Container<Number> remainder = dividend;
 
        for (SizeT i = 0; i <= n - m; ++i) {
            const Number factor = remainder[i] / divisor[divisor_start_idx];
            quotient[i] = factor;
            for (SizeT j = 0; j < m; ++j) {
                remainder[i+j] -= factor * divisor[divisor_start_idx+j];
            }
        }
 
        remainder.erase(remainder.begin(), remainder.end() - m + 1);
        return {quotient, remainder};
    }
}