polynom.hpp

Go to the documentation of this file.

/* /////////////////////////////////////////////////////////////////////////////
//
//  File:       polynom.hpp
//  Created:    2005/9/30    8:38
//
//  Author: Andrey Somsikov
*/

#ifndef __ARAGELI_polynom_hpp__
#define __ARAGELI_polynom_hpp__

#include <cstddef>
#include <iostream>
#include <vector>
#include <iomanip>
#include <algorithm>
#include <cmath>

/*
#include <stdlib.h>
#include <ios>
#include <vector>
#include "arageli_defs.h"
#include "matrices.h"
*/

#include "config.hpp"
#include "refcntr.hpp"
#include "iteradapt.hpp"
#include "factory.hpp"
#include "powerest.hpp"
#include "exception.hpp"
#include "_utility.hpp"

namespace Arageli
{

template <typename F>
class polynom
{
public:
    typedef F coef_type;
    typedef std::pair<polynom<F>, polynom<F> > div_result_type;
protected:
    template <typename F>
    void fixLeadingZeroes()
    {
        std::size_t k = c.size()-1;
        while (c[k] == factory<coef_type>::null() && k > 0)
            k--;
        c.resize(k+1);
    }

public:
    std::size_t size() const { return c.size(); }
    std::size_t deg() const { return size()-1; }

    void setLe(const std::vector<F> &coefficients)
    {
        std::size_t k = coefficients.size()-1;
        while (coefficients[k] == factory<coef_type>::null() && k > 0)
            k--;
        c.reserve(k + 1);
        for (std::size_t i = 0; i <= k; i++)
            c.push_back(coefficients[i]);
    }

    void setLe(const vector<F> &coefficients)
    {
        std::size_t k = coefficients.get_n()-1;
        while (coefficients[k] == factory<coef_type>::null() && k > 0)
            k--;
        c.reserve(k + 1);
        for (std::size_t i = 0; i <= k; i++)
            c.push_back(coefficients[i]);
    }

    void setBe(const std::vector<F> &coefficients)
    {

        int k = 0;
        while (coefficients[k] == factory<coef_type>::null() && k < (int)coefficients.size()-1)
            k++;
        c.reserve((int)coefficients.size() - k);
        for (int i = (int)coefficients.size() - 1; i >= k; i--)
            c.push_back(coefficients[i]);
    }

    void setBe(const vector<F> &coefficients)
    {

        int k = 0;
        while (coefficients[k] == factory<coef_type>::null() && k < (int)coefficients.get_n()-1)
            k++;
        c.reserve((int)coefficients.get_n() - k);
        for (int i = (int)coefficients.get_n() - 1; i >= k; i--)
            c.push_back(coefficients[i]);
    }

    const std::vector<F> &getData() const
    {
        return c;
    }

    div_result_type div (const polynom<F>& a) const
    {
        if ((*this).size() < a.size())
            return div_result_type(polynom<F>(), *this);
        polynom<F> q;
        q.c.resize((*this).size() - a.size() + 1, factory<coef_type>::null());
        polynom<F> r;
        r.c = (*this).c;
        std::size_t rSize = r.size();

        while (rSize >= a.size())
        {
            F qValue = r.c[rSize - 1] / a[a.size()-1];
            q.c[rSize - a.size()] = qValue;
            for (std::size_t i = 0, j = rSize - a.size(); i < a.size(); i++, j++)
                r.c[j] -= qValue * a[i];
            if (factory<coef_type>::null() != r.c[rSize - 1])
                throw "can't divide";
            rSize--;
        }
        q.fixLeadingZeroes<F>();
        r.fixLeadingZeroes<F>();

        return div_result_type(q, r);
    }

    polynom<F>& operator = (const polynom<F>& copy) 
    { 
        c = copy.c;
        return *this; 
    }

    F& operator [] (std::size_t i)
    {
        return c[i];
    }

    const F& operator [] (std::size_t i) const
    {
        return c[i];
    }

    polynom<F> operator - ( ) const
    {
        polynom<F> res(*this);
        for (std::size_t i = 0; i < res.size(); i++)
            res[i] = -res[i];
        return res;
    }

    template <typename F>
        bool operator < (const polynom<F>& a) const
    {
        if ((*this).size() < a.size())
            return true;
        if ((*this).size() == a.size())
            return (*this)[(*this).size() - 1] < a[a.size()-1];
        return false;
    }

    template <typename F>
        bool operator <= (const polynom<F>& a) const
    {
        if ((*this).size() < a.size())
            return true;
        if ((*this).size() == a.size())
            return (*this)[(*this).size() - 1] <= a[a.size()-1];
        return false;
    }

    template <typename F>
        bool operator > (const polynom<F>& a) const
    {
        if ((*this).size() > a.size())
            return true;
        if ((*this).size() == a.size())
            return (*this)[(*this).size() - 1] > a[a.size()-1];
        return false;
    }

    template <typename F>
        bool operator >= (const polynom<F>& a) const
    {
        if ((*this).size() > a.size())
            return true;
        if ((*this).size() == a.size())
            return (*this)[(*this).size() - 1] >= a[a.size()-1];
        return false;
    }

    template <typename F>
        polynom<F> operator + (const polynom<F>& a) const
    {
        polynom<F> res;
        const polynom<F> *op;

        if (a.size() > (*this).size())
        {
            res = polynom<F>(a);
            op = this;
        }
        else
        {
            res = polynom<F>(*this);
            op = &a;
        }

        for (std::size_t i = 0; i < op->size(); i++)
            res[i] += (*op)[i];

        res.fixLeadingZeroes<F>();
        return res;
    }

    template <typename F>
        polynom<F> operator - (const polynom<F>& a) const
    {
        return *this + (-a);
    }

    template <typename F>
        polynom<F> operator * (const polynom<F>& a) const
    {
        std::size_t i, j;
        std::vector<F> res(a.size() + (*this).size() - 1, factory<coef_type>::null());

        for (i = 0; i < a.size(); i++)
            for (j = 0; j < (*this).size(); j++)
            {
                res[i+j] += a[i] * (*this)[j];
            }

            std::size_t k = res.size() - 1;
            while (res[k] == factory<coef_type>::null() && k > 0)
                k--;
            res.resize(k + 1);

            for (i = 0; i < res.size()/2; i++)
                std::swap(res[i], res[res.size() - i - 1]);

            return polynom<F>(res);
    }

    template <typename F>
        polynom<F> operator / (const polynom<F>& a) const
    {
        return div(a).first;
    }

    template <typename F>
        polynom<F> operator % (const polynom<F>& a) const
    {
        return div(a).second;
    }

    polynom<F>& operator += (const polynom<F>& a)
    {
        *this = *this + a;
        return *this;
    }

    polynom<F>& operator -= (const polynom<F>& a)
    {
        *this = *this - a;
        return *this;
    }

    polynom<F>& operator *= (const polynom<F>& a)
    {
        *this = *this * a;
        return *this;
    }

    polynom<F>& operator /= (const polynom<F>& a)
    {
        *this = *this / a;
        return *this;
    }

    int operator == (const polynom<F>& a) const
    {
        return (*this).c == a.c;
    }

    int operator != (const polynom<F>& a) const
    {
        return !((*this) == a);
    }

protected:
    std::vector<F> c;
public:
    polynom<F>()
    {
        c.push_back(factory<coef_type>::null());
    }

    polynom<F>(const F& value, std::size_t deg = 0)
    {
        c.resize(deg+1, factory<coef_type>::null());
        c[deg] = value;
    }

    polynom<F>(const polynom<F>& copy)
    {
        c = copy.c;
    }

    polynom<F>(const std::vector<F> &coeff, bool le = false)
    {
        if (le)
            setLe(coeff);
        else
            setBe(coeff);
    }

    polynom<F>(const vector<F> &coefficients, bool le = false)
    {
        if (le)
            setLe(coefficients);
        else
            setBe(coefficients);
    }

    polynom<F>(const F *coefficients, std::size_t count, bool le = false)
    {
        std::vector<F> vCoeff;
        vCoeff.reserve(count);
        for (std::size_t i = 0; i < count; i++)
            vCoeff.push_back(coefficients[i]);
        if (le)
            setLe(vCoeff);
        else
            setBe(vCoeff);
    }

/*
    polynom<T> (const char* s)
    {
        std::istringstream buf(s);
        // WARNING. Correct if there are no virtual function.
        buf >> *this;
    }
*/

    ~polynom<F>()
    {
    }
public:
    static polynom<F> constant(const F& value)
    {
        return polynom<F>(value, 0);
    }
};

/*

template <typename F, typename FT>
std::istream& input_list
(
    std::istream& in,
    polynom<F>& x,
    const char* first_bracket = "(",
    const char* second_bracket = ")",
    const char* separator = ","
);

template <typename T, typename FT>
std::istream& operator>> (std::istream& out, polynom<T>& x)
{ return input_list(out, x); }
*/

template <typename F, typename FT>
std::ostream& operator << (std::ostream & s, const polynom<F> &v)
{
    std::size_t i = v.deg() + 1;
    s << "(";
    do {
        i--;
        s << v[i] << ", ";              
    } while(i != 0);
    s << ")";
    /*
    if (s.fail ())
    matrix_error (st_output_error, "polynom::output error");
    */
    return s;
}

template <typename F, typename FT>
polynom<F> abs(const polynom<F> &maxVal)
{
    polynom<F> res(factory<coef_type>::unit(), maxVal.deg());
    for (std::size_t i = 0; i < res.size(); i++)
        res[i] = abs(maxVal[i]);
    return res;
}

template <typename F, typename FT>
polynom<F> normalize(const polynom<F> &val)
{
    polynom<F> res(val);
    F d = val[val.size()-1];
    for (std::size_t i = 0; i < val.size(); i++)
        res[i] = val[i] / d;
    return res;
}


} // namespace Arageli


#ifdef ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE
#define ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_POLYNOM
#include "polynom.cpp"
#undef  ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_POLYNOM
#endif

#endif /*__ARAGELI_polynom_hpp__*/
/* End of file polynom.hpp */

---