RDE_lib2/dyadic_interval.h

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

Copyright 2010 Terry Lyons, Stephen Buckley, Djalil Chafai, 
Greg Gyurkó and Arend Janssen. 

Distributed under the terms of the GNU General Public License, 
Version 3. (See accompanying file License.txt)

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



#pragma once

#ifndef __DYADIC_INTERVAL__
#define __DYADIC_INTERVAL__


#include <deque>
#include <iostream>
#include "dyadic.h"


class dyadic_interval : dyadic // use private derivation to avoid spurious conversions etc.
{
public:
        using dyadic::k;
        using dyadic::n;
        using dyadic::operator ++;

        template<typename S>
                dyadic_interval(S s): dyadic(s)
        {}

        template<typename S, typename T>
                dyadic_interval(S s, T t): dyadic(s,t)
        {}

        dyadic_interval(void)
                : dyadic()
        {
        }

        ~dyadic_interval(void)
        {
        }

        const dyadic & inf(void) const
        {
                return (const dyadic &) *this;
        }

        dyadic sup(void) const
        {
                return ++ dyadic( *this );
        }

        dyadic_interval & flip_interval(void)
        {
                k = k^1;
                return *this;
        }

        int aligned (void) const
        {
                return (k+1)%2;
        }

        dyadic_interval & shrink_interval_left(const unsigned int Arg = 1)
        {
                k <<= Arg;
                n += Arg;
                return *this;
        }

        dyadic_interval & shrink_interval_right(void)
        {
                ++(k <<= 1);
                n += 1;
                return *this;
        }

        dyadic_interval & expand_interval(const unsigned int Arg = 1)
        {
                (dyadic::operator< (dyadic())) ? ++k : k;
                k >>= Arg; // rounding always towards zero so must make sure that one is rounding the dyadic on the correct side of zero
                n -= Arg;
                (dyadic::operator< (dyadic())) ? --k : k;
                return *this;
        }


        bool operator < (const dyadic_interval & Arg) const 
        {
                // tree leaf order
                return (0 < (Arg.n - n))        
                        ? k <= (Arg.k >> (Arg.n - n)) 
                        : (k >> (n - Arg.n)) < Arg.k ; 
        }

        bool operator != (const dyadic_interval & Arg) const 
        {
                return (k != Arg.k || n != Arg.n);
        }
        
        bool operator == (const dyadic_interval & Arg) const
        {
                return (k == Arg.k && n == Arg.n);
        }

        bool contains (const dyadic_interval & Arg) const 
        {
                return 
                        // Arg must be shorter to be contained
                        (0 <= (Arg.n - n))      

                        ? k == (Arg.k >> (Arg.n - n)) 
                        : false ; 
        }



        static std::deque<dyadic_interval> to_dyadic_intervals(double inf, double sup, int tolerance)
        {
                dyadic_interval m_inf = dyadic((int)floor(ldexp(inf,tolerance)),tolerance);
                dyadic_interval m_sup = dyadic((int)floor(ldexp(sup,tolerance)),tolerance);
                std::deque<dyadic_interval> Ans;
                dyadic_interval temp(m_inf);
                while (!(temp.sup() > m_sup))
                {
                        temp.expand_interval();
                        if ( temp.inf() < m_inf )
                        {
                                Ans.push_back(temp.shrink_interval_right());
                                m_inf=temp.inf();
                        }
                }
                while ((temp.inf() < sup)&&(temp.sup() > m_sup))
                {
                        temp.shrink_interval_left();
                        if(!(temp.sup() > m_sup))
                        {
                                Ans.push_back(temp);
                                ++ (dyadic &) temp;
                        }
                }
                //all done
                return Ans;
        }

    static dyadic_interval dyadic_bracket(const double Arg, const int AbsolutePrecision=0)
        {
                dyadic_interval temp;
        temp.n=AbsolutePrecision;
                temp.k= (int) floor(ldexp(Arg,AbsolutePrecision));
                return temp;
        }
        
        friend inline std::ostream& operator << (std::ostream & os , const dyadic_interval & rhs)
        {
                os << "[" << (double) rhs.inf() << ", " << (double) rhs.sup() << ") ";
                return os;
        }
};

#endif // __DYADIC_INTERVAL__