algebraic.hpp

Go to the documentation of this file.

/*****************************************************************************
    
    algebraic.hpp -- Algebraic number declaration.

    This file is a part of the Arageli library.

    Copyright (C) Nikolai Yu. Zolotykh, 1999--2006
    Copyright (C) Sergey S. Lyalin, 2006
    University of Nizhni Novgorod, Russia

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

#ifndef _ARAGELI_algebraic_hpp_
#define _ARAGELI_algebraic_hpp_

#include "config.hpp"

#include <iostream>
#include <sstream>

#include "_utility.hpp"
#include "factory.hpp"
#include "intalg.hpp"
#include "polyalg.hpp"
#include "sturm.hpp"
#include "algebrslt.hpp"
#include "sparse_polynom.hpp"
#include "rational.hpp"
#include "interval.hpp"

#include "std_import.hpp"


namespace Arageli
{


template <typename TP, typename TS, typename Poly, typename Seg>
struct algebraic_config_default
{
    static void init (const Poly& poly, const Seg& seg) { normalize(poly, seg); }
    static void get_value (const Poly& poly, const Seg& seg) {}
    static void before_oper (const Poly& poly, const Seg& seg) {}
    static void after_oper (const Poly& poly, const Seg& seg) { normalize(poly, seg); }
    
    static bool is_normal (const Poly& poly, const Seg& seg)
    {
        vector<Poly, false> ss;
        sturm_diff_sys(poly, ss);
        return sturm_number_roots(ss, seg) == 1;
        //if(seg.is_point())
        //  return is_null(poly.subs(seg.first()));
        //else
        //{
        //  vector<Seg, false> lims;
        //  sturm<TS>(poly, lims, seg);
        //  return lims.size() == 1;
        //}
    }
    
    static void normalize (const Poly& poly, const Seg& seg)
    {
        ARAGELI_ASSERT_0(is_normal(poly, seg));
    }
};


template
<
    typename TP = rational<>,
    typename TS = TP, 
    typename Poly = sparse_polynom<TP>,
    typename Seg = interval<TS>,
    typename Config = algebraic_config_default<TP, TS, Poly, Seg>
>
class algebraic
{
    Poly poly_m;    // polynomial with one root in seg segment that identifies the number
    Seg seg_m;  // segment where the number located
    Config config_m;    // configurator

    template <typename TP2, typename TS2, typename Poly2, typename Seg2, typename Cfg2>
    friend class algebraic;

public:

    typedef TP base_type;
    typedef TS interval_limit_type;
    typedef Poly polynom_type;
    typedef Seg interval_type;
    typedef Config config_type;

    algebraic () : poly_m(null<TP>()), seg_m(null<TS>()) {}

    algebraic (const TS& x) : seg_m(x)
    {
        poly_m += typename Poly::monom(unit<TP>(), 1);
        poly_m -= TP(x);
    }

    algebraic (const Poly& px) : poly_m(px)
    {
        vector<Seg, false> lims;
        sturm(px, lims);
        ARAGELI_ASSERT_0(!lims.is_empty());
        seg_m = lims.front();
    }

    // Creates number as the root of px in segment s.
    algebraic (const Poly& px, const Seg& s) : poly_m(px), seg_m(s)
    { config_m.init(poly_m, seg_m); }


    //algebraic (const char* s);


    operator double () const
    {
        // WARNING! TEMPORARY!
        interval<rational<> > fseg = seg_m;
        interval_root_precise(poly_m, fseg, rational<>("1/100000000000"));
        return fseg.first();
    }


    const Poly& polynom () const
    {
        config_m.get_value(poly_m, seg_m);
        return poly_m;
    }

    const Seg& segment () const { return seg_m; }

    Poly& polynom ()
    {
        config_m.get_value(poly_m, seg_m);
        return poly_m;
    }

    Seg& segment () { return seg_m; }


    bool is_normal () const { return config_m.is_normal(poly_m, seg_m); }

    void normalize () const { config_m.normalize(poly_m, seg_m); }


    algebraic operator- () const
    {
        if(is_null(poly_m) && is_null(seg_m))return *this;
        
        ARAGELI_ASSERT_0(is_normal());
        Poly poly_res; Seg seg_res;
        
        algebraic_minus
        (
            Poly(unit<TP>(), 1), seg(null<TS>()),
            poly_m, seg_m, poly_res, seg_res
        );

        return algebraic(poly_res, seg_res);
    }
    
    const algebraic& operator+ () const { return *this; }

    algebraic& operator++ ()
    {
        *this += algebraic(unit<TP>());
        return *this;
    }
    
    algebraic& operator-- ()
    {
        *this -= algebraic(unit<TP>());
        return *this;
    }
    
    algebraic operator++ (int) { algebraic t = *this; operator++(); return t; }
    algebraic operator-- (int) { algebraic t = *this; operator--(); return t; }
    
    template <typename TP2, typename TS2, typename Poly2, typename Seg2, typename Cfg2>
    algebraic& operator+= (const algebraic<TP2, TS2, Poly2, Seg2, Cfg2>& x)
    {
        ARAGELI_ASSERT_0(is_normal() && x.is_normal());
        algebraic_plus(poly_m, seg_m, x.poly_m, x.seg_m, poly_m, seg_m);
        ARAGELI_ASSERT_0(is_normal());
        return *this;
    }
    
    template <typename TP2, typename TS2, typename Poly2, typename Seg2, typename Cfg2>
    algebraic& operator-= (const algebraic<TP2, TS2, Poly2, Seg2, Cfg2>& x)
    {
        ARAGELI_ASSERT_0(is_normal() && x.is_normal());
        algebraic_minus(poly_m, seg_m, x.poly_m, x.seg_m, poly_m, seg_m);
        ARAGELI_ASSERT_0(is_normal());
        return *this;
    }
    
    template <typename TP2, typename TS2, typename Poly2, typename Seg2, typename Cfg2>
    algebraic& operator*= (const algebraic<TP2, TS2, Poly2, Seg2, Cfg2>& x)
    {
        ARAGELI_ASSERT_0(is_normal() && x.is_normal());
        algebraic_multiplies(poly_m, seg_m, x.poly_m, x.seg_m, poly_m, seg_m);
        ARAGELI_ASSERT_0(is_normal());
        return *this;
    }
    
    template <typename TP2, typename TS2, typename Poly2, typename Seg2, typename Cfg2>
    algebraic& operator/= (const algebraic<TP2, TS2, Poly2, Seg2, Cfg2>& x)
    {
        ARAGELI_ASSERT_0(is_normal() && x.is_normal());
        algebraic_divides(poly_m, seg_m, x.poly_m, x.seg_m, poly_m, seg_m);
        ARAGELI_ASSERT_0(is_normal());
        return *this;
    }

    bool is_null () const { return is_equal_const(null<TP>()); }
    
    bool is_unit () const { return is_equal_const(unit<TP>()); }
    bool is_opposite_unit () const { return is_equal_const(opposite_unit<TP>()); }
    

    bool operator! () const { return is_null(); }
    operator bool () const { return !is_null(); }

private:

    bool is_equal_const (const TP& x) const
    {
        ARAGELI_ASSERT_0(is_normal());
        return
            x >= seg_m.first() && x <= seg_m.second() &&
            Arageli::is_null(poly_m.subs(x));
    }


};


template <typename TP, typename TS, typename Poly, typename Seg, typename Cfg>
struct factory<algebraic<TP, TS, Poly, Seg, Cfg> >
{
private:

    typedef algebraic<TP, TS, Poly, Seg, Cfg> T;

public:

    static const bool is_specialized = true;

    static const T& unit ()
    {
        static const T unit_s = T(Arageli::unit<TS>());
        return unit_s;
    }
    
    static const T& unit (const T& x)
    { return unit(); }

    static const T& opposite_unit ()
    {
        static const T opposite_unit_s = T(Arageli::opposite_unit<TS>());
        return opposite_unit_s;
    }
    
    static const T& opposite_unit (const T& x)
    { return opposite_unit(); }

    static const T& null ()
    {
        static const T null_s = T(Arageli::null<TS>());
        return null_s;
    }
    
    static const T& null (const T& x)
    { return null(); }

};


#define ARAGELI_ALGEBRAIC_BINARY_OP(OP, OPASS)                                      \
    template                                                                        \
    <                                                                               \
        typename TP1, typename TS1, typename Poly1, typename Seg1, typename Cfg1,   \
        typename TP2, typename TS2, typename Poly2, typename Seg2, typename Cfg2    \
    >                                                                               \
    inline algebraic<TP1, TS1, Poly1, Seg1, Cfg1> operator OP                       \
    (                                                                               \
        algebraic<TP1, TS1, Poly1, Seg1, Cfg1> a,                                   \
        const algebraic<TP2, TS2, Poly2, Seg2, Cfg2>& b                             \
    )                                                                               \
    { return a OPASS b; }

ARAGELI_ALGEBRAIC_BINARY_OP(+, +=)
ARAGELI_ALGEBRAIC_BINARY_OP(-, -=)
ARAGELI_ALGEBRAIC_BINARY_OP(*, *=)
ARAGELI_ALGEBRAIC_BINARY_OP(/, /=)
ARAGELI_ALGEBRAIC_BINARY_OP(%, %=)


#define ARAGELI_ALGEBRAIC_LOGIC_OP(OP)                              \
    template                                                        \
    <                                                               \
        typename TP1, typename TS1, typename Poly1, typename Seg1, typename Cfg1,   \
        typename TP2, typename TS2, typename Poly2, typename Seg2, typename Cfg2    \
    >                                                               \
    inline bool operator OP                                         \
    (                                                                               \
        const algebraic<TP1, TS1, Poly1, Seg1, Cfg1>& a,                            \
        const algebraic<TP2, TS2, Poly2, Seg2, Cfg2>& b                             \
    )                                                                               \
    { return cmp(a, b) OP 0; }                      \
                                                                    \
    template                                                        \
    <                                                               \
        typename TP1, typename TS1, typename Poly1, typename Seg1, typename Cfg1,   \
        typename TS2    \
    >                                                               \
    inline bool operator OP                                         \
    (                                                                               \
        const algebraic<TP1, TS1, Poly1, Seg1, Cfg1>& a,                            \
        const TS2& b                                \
    )                                                                               \
    { return cmp(a, algebraic<TP1, TS2, Poly1, Seg1, Cfg1>(b)) OP 0; }                      \
                                                                    \
    template                                                        \
    <                                                               \
        typename TS1,   \
        typename TP2, typename TS2, typename Poly2, typename Seg2, typename Cfg2    \
    >                                                               \
    inline bool operator OP                                         \
    (                                                                               \
        const TS1& a,                           \
        const algebraic<TP2, TS2, Poly2, Seg2, Cfg2>& b                             \
    )                                                                               \
    { return cmp(algebraic<TP2, TS1, Poly2, Seg2, Cfg2>(a), b) OP 0; }



ARAGELI_ALGEBRAIC_LOGIC_OP(==)
ARAGELI_ALGEBRAIC_LOGIC_OP(!=)
ARAGELI_ALGEBRAIC_LOGIC_OP(<)
ARAGELI_ALGEBRAIC_LOGIC_OP(<=)
ARAGELI_ALGEBRAIC_LOGIC_OP(>)
ARAGELI_ALGEBRAIC_LOGIC_OP(>=)


template                                                        
<                                                               
    typename TP1, typename TS1, typename Poly1, typename Seg1, typename Cfg1,   
    typename TP2, typename TS2, typename Poly2, typename Seg2, typename Cfg2    
>                                                               
int cmp                     
(                                                                               
    const algebraic<TP1, TS1, Poly1, Seg1, Cfg1>& a,                            
    const algebraic<TP2, TS2, Poly2, Seg2, Cfg2>& b                             
)
{
    algebraic<TP1, TS1, Poly1, Seg1, Cfg1> c = a - b;
    if(c.is_null())return 0;
    int lsign = sign(c.polynom().subs(c.segment().first()));
    while(sign(c.segment().first())*sign(c.segment().second()) <= 0)
        interval_root_dichotomy(c.polynom(), lsign, c.segment());
    return sign(c.segment().first());
}


template                                                        
<                                                               
    typename TP1, typename TS1, typename Poly1, typename Seg1, typename Cfg1,   
    typename TP2
>                                                               
inline int cmp                      
(                                                                               
    const algebraic<TP1, TS1, Poly1, Seg1, Cfg1>& a,                            
    const TP2& b                                
)
{ return cmp(a, algebraic<TP2, TS1, Poly1, Seg1, Cfg1>(b)); }


template                                                        
<                                                               
    typename TP1,   
    typename TP2, typename TS2, typename Poly2, typename Seg2, typename Cfg2    
>                                                               
inline int cmp                      
(                                                                               
    const TP1& a,                           
    const algebraic<TP2, TS2, Poly2, Seg2, Cfg2>& b                             
)
{ return cmp(algebraic<TP1, TS2, Poly2, Seg2, Cfg2>(a), b); }


template
<
    typename Out,
    typename TP1, typename TS1, typename Poly1, typename Seg1, typename Cfg1
>
inline Out& operator<< (Out& out, const algebraic<TP1, TS1, Poly1, Seg1, Cfg1>& x)
{
    out << "(" << x.polynom() << ", " << x.segment() << ")";
    return out;
}


//template <typename TP1, typename TS1, typename Poly1, typename Seg1, typename Cfg1>
//algebraic<TP1, TS1, Poly1, Seg1, Cfg1>::algebraic (const char* s)
//{
//  std::istringstream buf(s);
//  buf >> *this;
//}




} // namesapce Arageli


//#ifdef ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE
//  #define ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_ALGEBRAIC
//  #include "algebraic.cpp"
//  #undef  ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_ALGEBRAIC
//#endif

#endif

---