resultant.cpp

Go to the documentation of this file.

/*****************************************************************************
    
    resultant.cpp -- See declarations in resultant.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_resultant)

#include "vector.hpp"

#include "resultant.hpp"


namespace Arageli
{


template <typename P1, typename P2, typename M>
M& sylvester_fill (const P1& p1, const P2& p2, M& res)
{
    typedef typename M::size_type size_type;
    typedef typename P1::monom_const_iterator P1I;
    typedef typename P2::monom_const_iterator P2I;

    size_type
        d1 = p1.degree(),
        d2 = p2.degree(),
        s = d1 + d2;

    res.resize(s, s);

    size_type i = 0;
    for(; i != d2; ++i)
    {
        size_type margin = i + d1;
        for
        (
            P1I p1i = p1.monoms_begin(), p1end = p1.monoms_end();
            p1i != p1end; ++p1i
        )
            res(i, margin - p1i->degree()) = p1i->coef();
    }

    for(; i != s; ++i)
        for
        (
            P2I p2i = p2.monoms_begin(), p2end = p2.monoms_end();
            p2i != p2end; ++p2i
        )
            res(i, i - p2i->degree()) = p2i->coef();

    return res;
}


template <typename P1, typename P2, typename M, typename MTFactory>
M& sylvester (const P1& p1, const P2& p2, M& res, MTFactory fctr)
{
    res.assign
    (
        p1.degree() + p2.degree(), p1.degree() + p2.degree(),
        res.is_empty() ? fctr.null() : fctr.null(res(0, 0))
    );

    return sylvester_fill(p1, p2, res);
}


template <typename P1, typename P2, typename SRChain, typename PCFactory>
void subresultants_nonmodular
(const P1& p1, const P2& p2, SRChain& s, PCFactory fctr)
{
    // Source: Buhberger and others, p. 168, alg. 4.1
    // We have a problem with sign of the result subresultants.
    
    typedef typename SRChain::difference_type index;
    typedef typename SRChain::element_type P;   // maybe there is better chose
    typedef typename P::coef_type Coef;
    typedef typename P::degree_type Degree;

    typename P1::degree_type m = p1.degree();
    typename P2::degree_type n = p2.degree();

    index j;
    if(m > n)j = m-1;
    else j = n;

    s.resize(j+2);

    s[j+1] = p1;
    s[j] = p2;

    vector<Coef, false> rr(j+2);

    rr[j+1] = fctr.unit();

    for(;;)
    {
        P& sj = s[j];
        index r = sj.degree();  // P should return -1 when polynomial is null.

        for(index k = j-1; k >= r+1; --k)
            s[k] = null(s[k]);

        if(j > r && r >= 0)
            s[r] = (sj*pow(sj.leading_coef(), j-r)) / pow(rr[j+1], j-r);

        if(r <= 0)break;

        {
            P pq; Coef mult;
            polynom_pseudodivide_simple(s[j+1], sj, pq, s[r-1], mult);
        }

        s[r-1] /= pow(-rr[j+1], j-r+2);
        j = r-1;
        rr[j+1] = s[j+1].leading_coef();
    }

    //output_aligned(std::cout << "\nSub Resultants =\n", s);
}


template <typename P1, typename P2, typename PCFactory>
typename P1::coef_type resultant_nonmodular
(const P1& p1, const P2& p2, PCFactory fctr)
{
    // WARNING! Need to choose from P1::coef_type and P2::coef_type
    // insted just P1::coef_type as a type of the return value.

    vector<P1, false> s;    // WARNING! Need to choose from P1 and P2 insted just P1.
    subresultants_nonmodular(p1, p2, s, fctr);
    
    ARAGELI_ASSERT_2(s[0].degree() <= 0);

    //ARAGELI_DEBUG_EXEC_2(if(!is_null(s[0])))
    //{
    //  ARAGELI_DEBUG_EXEC_2(typename P1::coef_type _dres = resultant_mod1(p1, p2, fctr));
    //  ARAGELI_ASSERT_2(!is_null(_dres));
    //  ARAGELI_ASSERT_2
    //  (
    //      s[0].leading_coef()/sign(s[0].leading_coef().leading_coef()) ==
    //      _dres/sign(_dres.leading_coef())
    //  );
    //}
    //ARAGELI_DEBUG_EXEC_2(else)
    //{
    //  ARAGELI_ASSERT_2(is_null(resultant_mod1(p1, p2, fctr)));
    //}

    return s[0].leading_coef_cpy();
}


}


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

---