include/realroot/extended.hpp

Go to the documentation of this file.

/********************************************************************
 *   This file is part of the source code of the realroot library.
 *   Author(s): P. Dobrowolski, Warsaw University of Technology, Poland
 *   $Id: extended.hpp,v 1.1 2012/02/22 20:00:00 pdobrowolski Exp $
 ********************************************************************/
#ifndef REALROOT_EXTENDED_H_
#define REALROOT_EXTENDED_H_

#include <cassert>
#include <realroot/scalar.hpp>

namespace mmx {

// this must hold for two extended types to be operable
#define extended_check(A, B)                        \
    assert((A.c == 0 || B.c == 0 || A.c == B.c) &&  \
           (A.c != 0 || A.b == 0) &&                \
           (B.c != 0 || B.b == 0))

// return a common root extension of two potentially extended types
#define common_root(A, B) \
    (A.c != 0 ? A.c : B.c)

template<class NT_>
struct extended
{
    // a + b sqrt(c)
    // note: c must not be a square!
    typedef NT_             NT;
    typedef scalar<NT>      T;
    typedef extended<NT>    Self;

    T a, b, c;

    extended()
    {
    }

    extended(const T &a_)
        : a(a_),
          b(T(0)),
          c(T(0))
    {
    }

    extended(const T &a_, const T &b_, const T &c_)
        : a(a_),
          b(b_),
          c(c_)
    {
    }

    static int signof(const T &x)
    {
        if (x > 0)
            return 1;
        else if (x < 0)
            return -1;
        else
            return 0;
    }

    static int compare(const Self &lhs, const Self &rhs)
    {
        extended_check(lhs, rhs);
        Self delta = lhs - rhs;
        // equal
        if (delta.a == 0 && delta.b == 0) return 0;
        // simple
        if (delta.a == 0) return signof(delta.b);
        if (delta.b == 0) return signof(delta.a);
        // less
        if (delta.a < 0 && delta.b < 0) return -1;
        // greater
        if (delta.a > 0 && delta.b > 0) return 1;
        // less or greater, but mixed signs
        T aa = delta.a * delta.a;
        T bbc = delta.b * delta.b * delta.c;
        if (aa > bbc)
            return signof(delta.a);
        else
            return signof(delta.b);
        return 0;
    }
};

template<class NT> inline
bool operator == (const extended<NT> &lhs, const extended<NT> &rhs)
{
    extended_check(lhs, rhs);
    return lhs.a == rhs.a && lhs.b == rhs.b;
}

template<class NT> inline
bool operator != (const extended<NT> &lhs, const extended<NT> &rhs)
{
    return !(lhs == rhs);
}

template<class NT> inline
extended<NT> operator *(const extended<NT> &lhs, const extended<NT> &rhs)
{
    extended_check(lhs, rhs);
    typename extended<NT>::T root = common_root(lhs, rhs);
    return extended<NT>(lhs.a * rhs.a + lhs.b * rhs.b * root,
                        lhs.a * rhs.b + lhs.b * rhs.a,
                        root);
}

template<class NT> inline
extended<NT> operator /(const extended<NT> &lhs, const extended<NT> &rhs)
{
    extended_check(lhs, rhs);
    typename extended<NT>::T root = common_root(lhs, rhs);
    typename extended<NT>::T denominator = rhs.a * rhs.a - rhs.b * rhs.b * root;
    return extended<NT>((lhs.a * rhs.a - lhs.b * rhs.b * lhs.c) / denominator,
                        (lhs.b * rhs.a - lhs.a * rhs.b) / denominator,
                        root);
}

template<class NT> inline
extended<NT> operator +(const extended<NT> &lhs, const extended<NT> &rhs)
{
    extended_check(lhs, rhs);
    typename extended<NT>::T root = common_root(lhs, rhs);
    return extended<NT>(lhs.a + rhs.a,
                        lhs.b + rhs.b,
                        root);
}

template<class NT> inline
extended<NT> operator -(const extended<NT> &lhs, const extended<NT> &rhs)
{
    extended_check(lhs, rhs);
    typename extended<NT>::T root = common_root(lhs, rhs);
    return extended<NT>(lhs.a - rhs.a,
                        lhs.b - rhs.b,
                        root);
}

template<class NT> inline
extended<NT> operator -(const extended<NT> &lhs)
{
    return extended<NT>(-lhs.a,
                        -lhs.b,
                        lhs.c);
}

template<class NT> inline
extended<NT> &operator <<= (extended<NT> &x, long int s)
{
    x.a <<= s;
    x.b <<= s;
    return x;
}

template<class NT> inline
extended<NT> &operator += (extended<NT> &x, const extended<NT> &other)
{
    x = x + other;
    return x;
}

template<class NT> inline
extended<NT> &operator -= (extended<NT> &x, const extended<NT> &other)
{
    x = x - other;
    return x;
}

template<class NT> inline
extended<NT> &operator *= (extended<NT> &x, const extended<NT> &other)
{
    x = x * other;
    return x;
}

template<class NT> inline
extended<NT> &operator /= (extended<NT> &x, const extended<NT> &other)
{
    x = x / other;
    return x;
}

template<class NT> inline
bool operator <(const extended<NT> &lhs, const extended<NT> &rhs)
{
    return extended<NT>::compare(lhs, rhs) < 0;
}

template<class NT> inline
bool operator <=(const extended<NT> &lhs, const extended<NT> &rhs)
{
    return extended<NT>::compare(lhs, rhs) <= 0;
}

template<class NT> inline
bool operator >(const extended<NT> &lhs, const extended<NT> &rhs)
{
    return extended<NT>::compare(lhs, rhs) > 0;
}

template<class NT> inline
bool operator >=(const extended<NT> &lhs, const extended<NT> &rhs)
{
    return extended<NT>::compare(lhs, rhs) >= 0;
}

#undef extended_check
#undef common_root
} //namespace mmx

#endif // REALROOT_EXTENDED_H_

---