algebrslt.cpp

Go to the documentation of this file.

/*****************************************************************************
    
    algebrslt.cpp -- See declarations in algebrslt.hpp.

    This file is a part of the Arageli library.

    Copyright (C) Nikolai Yu. Zolotykh, 1999--2006
    Copyright (C) Sergey S. Lyalin, 2005--2006

*****************************************************************************/

#include "config.hpp"

#if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || \
    defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_algebrslt)

#include "exception.hpp"
#include "cmp.hpp"
#include "vector.hpp"
#include "sparse_polynom.hpp"
#include "polyalg.hpp"
#include "sturm.hpp"
#include "resultant.hpp"

#include "algebrslt.hpp"


namespace Arageli
{


template
<
    typename PolyVec,
    typename SegVec,
    typename Rslt,
    typename P,
    typename Seg,
    typename F
>
void algebraic_func_rslt
(
    const PolyVec& polyvec, // source polynomials
    SegVec segvec, // source segments
    Rslt rslt, // source resultant
    P& p, Seg& seg, // the result: polynomial and segment where root is located
    F f // the function that is performing on the source algebraic numbers
)
{
    // Implementation of this function is based on Loos, 1973. Source: Buhberger

    ARAGELI_ASSERT_0(polyvec.size() == segvec.size());

    //std::cout << "\n++++++++++ algebraic_func_rslt +++++++++++++\n";

    //typedef vector<Rslt, false> Rslt_factors;
    //typedef Rslt_factors::size_type size_type;

    //std::cout << "\nrslt = " << rslt;

    // Make squre-free resultant.   // WARNING! Is it needed? See strurm specification.
    Rslt gcd_rslt_df = gcd(rslt, diff(rslt));
    gcd_rslt_df /= gcd_rslt_df.leading_coef_cpy();  // WARNING! Is it needed?
    rslt /= gcd_rslt_df;

    //std::cout << "\nsf rslt = " << rslt;

    rslt = polynom_primpart(rslt);

    //std::cout << "\nprimpart rslt = " << rslt;

    vector<Rslt, false> ss;
    sturm_diff_sys(rslt, ss);
    typedef typename vector<Rslt, false>::size_type size_type;

    //output_aligned(std::cout << "\nrslt Sturm system =\n", ss);

    for(;;)
    {
        seg = f(segvec);

        //output_aligned(std::cout << "\nsegvec =\n", segvec);
        //std::cout << "seg = " << seg;

        size_type snr = sturm_number_roots(ss, seg);

        if(snr == 0)
            throw invalid_argument
            (/*
                "The operation cannot be performed on the specified "
                "arguments in algebraic number field"
            */);
        else if(snr == 1)   // Congratulations! The interval is found.
            break;
        else    // Make intervals more precise.
            for(size_type i = 0; i < polyvec.size(); ++i)
                if(!segvec[i].is_negligible())
                    interval_root_dichotomy
                    (
                        polyvec[i],
                        sign(polyvec[i].subs(segvec[i].first())),
                        segvec[i]
                    );
    }

    p = rslt;

    //std::cout << "\np = " << p;

    //std::cout << "\nprimpart p = " << p;
    //std::cout << "\n++++++++ END algebraic_func_rslt +++++++++++\n";
}


namespace _Internal
{

template <typename T>
struct vector_2_plus_algebraic_helper
{
    typename T::value_type operator() (const T& x) const
    {
        ARAGELI_ASSERT_0(x.size() == 2);
        return x.front() + x.back();
    }
};


template <typename T>
struct vector_2_minus_algebraic_helper
{
    typename T::value_type operator() (const T& x) const
    {
        ARAGELI_ASSERT_0(x.size() == 2);
        return x.front() - x.back();
    }
};


template <typename T>
struct vector_2_multiplies_algebraic_helper
{
    typename T::value_type operator() (const T& x) const
    {
        ARAGELI_ASSERT_0(x.size() == 2);
        return x.front() * x.back();
    }
};


template <typename T>
struct vector_2_divides_algebraic_helper
{
    typename T::value_type operator() (const T& x) const
    {
        ARAGELI_ASSERT_0(x.size() == 2);
        return x.front() / x.back();
    }
};


} // namespace _Internal


template
<
    typename P1, typename Seg1,
    typename P2, typename Seg2,
    typename P3, typename Seg3
>
void algebraic_plus
(
    const P1& p1, const Seg1& seg1,
    const P2& p2, const Seg2& seg2,
    P3& p3, Seg3& seg3
)
{
    typedef sparse_polynom<P1> Pxy;
    typedef vector<P1, false> Polyvec;
    typedef vector<Seg1, false> Segvec;

    //std::cout << "\np1 = " << p1 << "\nseg1 = " << seg1;
    //std::cout << "\np2 = " << p2 << "\nseg2 = " << seg2;

    Polyvec polyvec;
    polyvec.push_back(p1);
    polyvec.push_back(p2);
    
    Segvec segvec;
    segvec.push_back(seg1);
    segvec.push_back(seg2);
    
    P1 rslt = resultant(p1.subs(Pxy("((x)-x)")), p2);

    algebraic_func_rslt
    (
        polyvec, segvec, rslt, p3, seg3,
        _Internal::vector_2_plus_algebraic_helper<Segvec>()
    );
}


// WARNING! TEMPORARY.
template
<
    typename P1, typename Seg1,
    typename P2, typename Seg2,
    typename P3, typename Seg3
>
void algebraic_minus
(
    const P1& p1, const Seg1& seg1,
    const P2& p2, const Seg2& seg2,
    P3& p3, Seg3& seg3
)
{
    typedef sparse_polynom<P1> Pxy;
    typedef vector<P1, false> Polyvec;
    typedef vector<Seg1, false> Segvec;

    //std::cout << "\np1 = " << p1 << "\nseg1 = " << seg1;
    //std::cout << "\np2 = " << p2 << "\nseg2 = " << seg2;

    Polyvec polyvec;
    polyvec.push_back(p1);
    polyvec.push_back(p2);
    
    Segvec segvec;
    segvec.push_back(seg1);
    segvec.push_back(seg2);
    
    P1 rslt = resultant(p1.subs(Pxy("((x)+x)")), p2);

    algebraic_func_rslt
    (
        polyvec, segvec, rslt, p3, seg3,
        _Internal::vector_2_minus_algebraic_helper<Segvec>()
    );
}


// WARNING! TEMPORARY.
template
<
    typename P1, typename Seg1,
    typename P2, typename Seg2,
    typename P3, typename Seg3
>
void algebraic_multiplies
(
    const P1& p1, const Seg1& seg1,
    const P2& p2, const Seg2& seg2,
    P3& p3, Seg3& seg3
)
{
    typedef sparse_polynom<P1> Pxy;
    typedef vector<P1, false> Polyvec;
    typedef vector<Seg1, false> Segvec;

    //std::cout << "\np1 = " << p1 << "\nseg1 = " << seg1;
    //std::cout << "\np2 = " << p2 << "\nseg2 = " << seg2;

    Polyvec polyvec;
    polyvec.push_back(p1);
    polyvec.push_back(p2);
    
    Segvec segvec;
    segvec.push_back(seg1);
    segvec.push_back(seg2);
    
    Pxy p1n = reciprocal_poly(p1).subs(Pxy("x"));
    typedef typename P1::monom monom;
    for(typename Pxy::monom_iterator i = p1n.monoms_begin(); i != p1n.monoms_end(); ++i)
        i->coef() *= monom(unit<typename monom::coef_type>(), p1n.degree() - i->degree());
    //std::cout << "\np1n = " << p1n;
    
    P1 rslt = resultant(p1n, p2);

    algebraic_func_rslt
    (
        polyvec, segvec, rslt, p3, seg3,
        _Internal::vector_2_multiplies_algebraic_helper<Segvec>()
    );
}


// WARNING! TEMPORARY.
template
<
    typename P1, typename Seg1,
    typename P2, typename Seg2,
    typename P3, typename Seg3
>
void algebraic_divides
(
    const P1& p1, const Seg1& seg1,
    const P2& p2, Seg2 seg2,
    P3& p3, Seg3& seg3
)
{
    typedef sparse_polynom<P1> Pxy;
    typedef vector<P1, false> Polyvec;
    typedef vector<Seg1, false> Segvec;

    //std::cout << "\np1 = " << p1 << "\nseg1 = " << seg1;
    //std::cout << "\np2 = " << p2 << "\nseg2 = " << seg2;

    Polyvec polyvec;
    polyvec.push_back(p1);
    polyvec.push_back(p2);

    if(sign(seg2.first())*sign(seg2.second()) <= 0)
    {
        ARAGELI_ASSERT_0(!is_null(p2.subs(null<typename Seg2::value_type>())));
        int lsign = sign(p2.subs(seg2.first()));
        
        do
        {
            interval_root_dichotomy(p2, lsign, seg2);
        }while(sign(seg2.first())*sign(seg2.second()) <= 0);
    }
    
    Segvec segvec;
    segvec.push_back(seg1);
    segvec.push_back(seg2);
    
    P1 rslt = resultant(p1.subs(Pxy("((x)*x)")), p2);

    algebraic_func_rslt
    (
        polyvec, segvec, rslt, p3, seg3,
        _Internal::vector_2_divides_algebraic_helper<Segvec>()
    );
}



}


#endif  // #if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || ...

---