real.hpp

#ifndef BOOST_REAL_REAL_HPP
#define BOOST_REAL_REAL_HPP

#include <iostream>
#include <vector>
#include <algorithm>
#include <initializer_list>
#include <utility>

#include <real/real_exception.hpp>
#include <real/real_helpers.hpp>
#include <real/real_explicit.hpp>
#include <real/real_algorithm.hpp>


namespace boost {
    namespace real {

        class real {

            // Available operations
            enum class OPERATION {ADDITION, SUBTRACT, MULTIPLICATION};
            enum class KIND {EXPLICIT, OPERATION, ALGORITHM};

            KIND _kind;

            // Explicit number
            real_explicit _explicit_number;

            // Algorithmic number
            real_algorithm _algorithmic_number;

            // Composed number
            OPERATION _operation;
            real* _lhs_ptr = nullptr;
            real* _rhs_ptr = nullptr;

            // Precision
            unsigned int _maximum_precision = 0;

            void copy_operands(const real& other) {
                if (other._lhs_ptr != nullptr) {
                    this->_lhs_ptr = new real(*other._lhs_ptr);
                }

                if (other._rhs_ptr != nullptr) {
                    this->_rhs_ptr = new real(*other._rhs_ptr);
                }
            }

        public:

            static unsigned int maximum_precision;

            class const_precision_iterator {
            private:

                // Internal number to iterate
                real const* _real_ptr = nullptr;

                // Explicit number iterator
                boost::real::real_explicit::const_precision_iterator _explicit_it;

                // Algorithmic number iterator
                boost::real::real_algorithm::const_precision_iterator _algorithmic_it;

                // If the number is a composition, the const_precision_iterator uses the operand iterators
                const_precision_iterator* _lhs_it_ptr = nullptr;
                const_precision_iterator* _rhs_it_ptr = nullptr;

                void calculate_operation_boundaries() {

                    switch (this->_real_ptr->_operation) {

                        case OPERATION::ADDITION:
                            boost::real::helper::add_boundaries(
                                    this->_lhs_it_ptr->approximation_interval.lower_bound,
                                    this->_rhs_it_ptr->approximation_interval.lower_bound,
                                    this->approximation_interval.lower_bound
                            );

                            boost::real::helper::add_boundaries(
                                    this->_lhs_it_ptr->approximation_interval.upper_bound,
                                    this->_rhs_it_ptr->approximation_interval.upper_bound,
                                    this->approximation_interval.upper_bound
                            );
                            break;


                        case OPERATION::SUBTRACT:
                            boost::real::helper::subtract_boundaries(
                                    this->_lhs_it_ptr->approximation_interval.lower_bound,
                                    this->_rhs_it_ptr->approximation_interval.upper_bound,
                                    this->approximation_interval.lower_bound
                            );

                            boost::real::helper::subtract_boundaries(
                                    this->_lhs_it_ptr->approximation_interval.upper_bound,
                                    this->_rhs_it_ptr->approximation_interval.lower_bound,
                                    this->approximation_interval.upper_bound
                            );
                            break;


                        case OPERATION::MULTIPLICATION: {
                            bool lhs_positive = this->_lhs_it_ptr->approximation_interval.positive();
                            bool rhs_positive = this->_rhs_it_ptr->approximation_interval.positive();
                            bool lhs_negative = this->_lhs_it_ptr->approximation_interval.negative();
                            bool rhs_negative = this->_rhs_it_ptr->approximation_interval.negative();

                            if (lhs_positive && rhs_positive) { // Positive - Positive
                                boost::real::helper::multiply_boundaries(
                                        this->_lhs_it_ptr->approximation_interval.lower_bound,
                                        this->_rhs_it_ptr->approximation_interval.lower_bound,
                                        this->approximation_interval.lower_bound
                                );

                                boost::real::helper::multiply_boundaries(
                                        this->_lhs_it_ptr->approximation_interval.upper_bound,
                                        this->_rhs_it_ptr->approximation_interval.upper_bound,
                                        this->approximation_interval.upper_bound
                                );

                            } else if (lhs_negative && rhs_negative) { // Negative - Negative
                                boost::real::helper::multiply_boundaries(
                                        this->_lhs_it_ptr->approximation_interval.upper_bound,
                                        this->_rhs_it_ptr->approximation_interval.upper_bound,
                                        this->approximation_interval.lower_bound
                                );

                                boost::real::helper::multiply_boundaries(
                                        this->_lhs_it_ptr->approximation_interval.lower_bound,
                                        this->_rhs_it_ptr->approximation_interval.lower_bound,
                                        this->approximation_interval.upper_bound
                                );
                            } else if (lhs_negative && rhs_positive) { // Negative - Positive
                                boost::real::helper::multiply_boundaries(
                                        this->_lhs_it_ptr->approximation_interval.lower_bound,
                                        this->_rhs_it_ptr->approximation_interval.upper_bound,
                                        this->approximation_interval.lower_bound
                                );

                                boost::real::helper::multiply_boundaries(
                                        this->_lhs_it_ptr->approximation_interval.upper_bound,
                                        this->_rhs_it_ptr->approximation_interval.lower_bound,
                                        this->approximation_interval.upper_bound
                                );

                            } else if (lhs_positive && rhs_negative) { // Positive - Negative
                                boost::real::helper::multiply_boundaries(
                                        this->_lhs_it_ptr->approximation_interval.upper_bound,
                                        this->_rhs_it_ptr->approximation_interval.lower_bound,
                                        this->approximation_interval.lower_bound
                                );

                                boost::real::helper::multiply_boundaries(
                                        this->_lhs_it_ptr->approximation_interval.lower_bound,
                                        this->_rhs_it_ptr->approximation_interval.upper_bound,
                                        this->approximation_interval.upper_bound
                                );

                            } else { // One is around zero all possible combinations are be tested

                                boundary current_boundary;

                                // Lower * Lower
                                boost::real::helper::multiply_boundaries(
                                        this->_lhs_it_ptr->approximation_interval.lower_bound,
                                        this->_rhs_it_ptr->approximation_interval.lower_bound,
                                        current_boundary
                                );

                                this->approximation_interval.lower_bound = current_boundary;
                                this->approximation_interval.upper_bound = current_boundary;

                                // Upper * upper
                                boost::real::helper::multiply_boundaries(
                                        this->_lhs_it_ptr->approximation_interval.upper_bound,
                                        this->_rhs_it_ptr->approximation_interval.upper_bound,
                                        current_boundary
                                );

                                if (current_boundary < this->approximation_interval.lower_bound) {
                                    this->approximation_interval.lower_bound = current_boundary;
                                }

                                if (this->approximation_interval.upper_bound < current_boundary) {
                                    this->approximation_interval.upper_bound = current_boundary;
                                }

                                // Lower * upper
                                boost::real::helper::multiply_boundaries(
                                        this->_lhs_it_ptr->approximation_interval.lower_bound,
                                        this->_rhs_it_ptr->approximation_interval.upper_bound,
                                        current_boundary
                                );

                                if (current_boundary < this->approximation_interval.lower_bound) {
                                    this->approximation_interval.lower_bound = current_boundary;
                                }

                                if (this->approximation_interval.upper_bound < current_boundary) {
                                    this->approximation_interval.upper_bound = current_boundary;
                                }

                                // Upper * lower
                                boost::real::helper::multiply_boundaries(
                                        this->_lhs_it_ptr->approximation_interval.upper_bound,
                                        this->_rhs_it_ptr->approximation_interval.lower_bound,
                                        current_boundary
                                );

                                if (current_boundary < this->approximation_interval.lower_bound) {
                                    this->approximation_interval.lower_bound = current_boundary;
                                }

                                if (this->approximation_interval.upper_bound < current_boundary) {
                                    this->approximation_interval.upper_bound = current_boundary;
                                }
                            }
                            break;
                        }

                        default:
                            throw boost::real::none_operation_exception();
                    }
                }

            public:

                // Number approximation_interval boundaries
                boost::real::interval approximation_interval;

                const_precision_iterator() = default;

                const_precision_iterator(const const_precision_iterator& other) = default;

                explicit const_precision_iterator(real const* real_number) : _real_ptr(real_number) {

                    switch (this->_real_ptr->_kind) {

                        case KIND::EXPLICIT:
                                this->_explicit_it = this->_real_ptr->_explicit_number.cbegin();
                            this->approximation_interval = this->_explicit_it.approximation_interval;
                            break;

                        case KIND::ALGORITHM:
                            this->_algorithmic_it = this->_real_ptr->_algorithmic_number.cbegin();
                            this->approximation_interval = this->_algorithmic_it.approximation_interval;
                            break;

                        case KIND::OPERATION:
                            this->_lhs_it_ptr = new const_precision_iterator(this->_real_ptr->_lhs_ptr->cbegin());
                            this->_rhs_it_ptr = new const_precision_iterator(this->_real_ptr->_rhs_ptr->cbegin());
                            this->calculate_operation_boundaries();
                            break;
                    }
                }

                explicit const_precision_iterator(real const* real_number, bool cend) : _real_ptr(real_number) {

                    switch (this->_real_ptr->_kind) {

                        case KIND::EXPLICIT:
                            if (cend) {
                                this->_explicit_it = this->_real_ptr->_explicit_number.cend();
                            } else {
                                this->_explicit_it = this->_real_ptr->_explicit_number.cbegin();
                            }
                            this->approximation_interval = this->_explicit_it.approximation_interval;
                            break;

                        case KIND::ALGORITHM:
                            if (cend) {
                                this->_algorithmic_it = this->_real_ptr->_algorithmic_number.cend();
                            } else {
                                this->_algorithmic_it = this->_real_ptr->_algorithmic_number.cbegin();
                            }
                            this->approximation_interval = this->_algorithmic_it.approximation_interval;
                            break;

                        case KIND::OPERATION:
                            this->_lhs_it_ptr = new const_precision_iterator(this->_real_ptr->_lhs_ptr, cend);
                            this->_rhs_it_ptr = new const_precision_iterator(this->_real_ptr->_rhs_ptr, cend);
                            this->calculate_operation_boundaries();
                            break;
                    }
                }

                void operator++() {

                    switch (this->_real_ptr->_kind) {

                        case KIND::EXPLICIT:
                            ++this->_explicit_it;
                            this->approximation_interval = this->_explicit_it.approximation_interval;
                            break;

                        case KIND::ALGORITHM:
                            ++this->_algorithmic_it;
                            this->approximation_interval = this->_algorithmic_it.approximation_interval;
                            break;

                        case KIND::OPERATION:
                            // Composed number iteration
                            this->approximation_interval.lower_bound.clear();
                            this->approximation_interval.upper_bound.clear();

                            // Recursive iteration of the operands
                            ++(*this->_lhs_it_ptr);
                            ++(*this->_rhs_it_ptr);

                            // Final bound calculation
                            this->calculate_operation_boundaries();
                            break;
                    }
                };

                bool operator==(const const_precision_iterator& other) const {
                    // uninitialized iterators are never equals
                    if (this->_real_ptr == nullptr || other._real_ptr == nullptr) {
                        return false;
                    }

                    return (other._real_ptr == this->_real_ptr) && (other.approximation_interval == this->approximation_interval);
                }

                bool operator!=(const const_precision_iterator& other) const {
                    return !(*this == other);
                }
            };

            real() = default;

            real(const real& other)  :
                    _kind(other._kind),
                    _explicit_number(other._explicit_number),
                    _algorithmic_number(other._algorithmic_number),
                    _operation(other._operation) { this->copy_operands(other); };

            real(const std::string& number)
                    : _kind(KIND::EXPLICIT), _explicit_number(number) {}

            real(std::initializer_list<int> digits)
                    : _kind(KIND::EXPLICIT), _explicit_number(digits, digits.size()) {}


            real(std::initializer_list<int> digits, bool positive)
                    : _kind(KIND::EXPLICIT), _explicit_number(digits, digits.size(), positive) {}

            real(std::initializer_list<int> digits, int exponent)
                    : _kind(KIND::EXPLICIT), _explicit_number(digits, exponent) {};

            real(std::initializer_list<int> digits, int exponent, bool positive)
                    : _kind(KIND::EXPLICIT), _explicit_number(digits, exponent, positive) {};

            real(int (*get_nth_digit)(unsigned int), int exponent)
                    : _kind(KIND::ALGORITHM), _algorithmic_number(get_nth_digit, exponent) {}

            real(int (*get_nth_digit)(unsigned int),
                 int exponent,
                 bool positive)
                    : _kind(KIND::ALGORITHM),
                      _algorithmic_number(get_nth_digit, exponent, positive) {}

            ~real() {
                delete this->_lhs_ptr;
                this->_lhs_ptr = nullptr;

                delete this->_rhs_ptr;
                this->_rhs_ptr = nullptr;
            }

            unsigned int max_precision() const {
                if (this->_maximum_precision == 0) {
                    return boost::real::real::maximum_precision;
                }

                return this->_maximum_precision;
            }

            void set_maximum_precision(unsigned int maximum_precision) {
                this->_maximum_precision = maximum_precision;
            }

            const_precision_iterator cbegin() const {
                return const_precision_iterator(this);
            }

            const_precision_iterator cend() const {
                return const_precision_iterator(this, true);
            }

            /************** Operators ******************/

            int operator[](unsigned int n) const {
                int result;

                switch (this->_kind) {

                    case KIND::EXPLICIT:
                        result = this->_explicit_number[n];
                        break;

                    case KIND::ALGORITHM:
                        result = this->_algorithmic_number[n];
                        break;

                    case KIND::OPERATION:
                        throw boost::real::invalid_representation_exception();
                        break;
                }

                return result;
            };

            real& operator+=(const real& other) {
                this->_lhs_ptr = new real(*this);
                this->_rhs_ptr = new real(other);
                this->_kind = KIND::OPERATION;
                this->_operation = OPERATION::ADDITION;
                return *this;
            }

            real operator+(const real& other) const {
                real result = *this;
                result += other;
                return result;
            }

            real& operator-=(const real& other) {
                this->_lhs_ptr = new real(*this);
                this->_rhs_ptr = new real(other);
                this->_kind = KIND::OPERATION;
                this->_operation = OPERATION::SUBTRACT;
                return *this;
            }

            real operator-(const real& other) const {
                real result = *this;
                result -= other;
                return result;
            }

            real& operator*=(const real& other) {
                this->_lhs_ptr = new real(*this);
                this->_rhs_ptr = new real(other);
                this->_kind = KIND::OPERATION;
                this->_operation = OPERATION::MULTIPLICATION;
                return *this;
            }

            real operator*(const real& other) const {
                real result = *this;
                result *= other;
                return result;
            }

            real& operator=(const real& other) {
                this->_kind = other._kind;
                this->_explicit_number = other._explicit_number;
                this->_operation = other._operation;
                this->copy_operands(other);
                return *this;
            }

            real& operator=(const std::string& number) {
                *this = real(number);
                return *this;
            }

            bool operator<(const real& other) const {
                auto this_it = this->cbegin();
                auto other_it = other.cbegin();

                unsigned int current_precision = std::max(this->max_precision(), other.max_precision());
                for (unsigned int p = 0; p < current_precision; ++p) {
                    // Get more precision
                    ++this_it;
                    ++other_it;

                    bool this_full_precision = this_it.approximation_interval.is_a_number();
                    bool other_full_precision = other_it.approximation_interval.is_a_number();
                    if (this_full_precision && other_full_precision) {
                        return this_it.approximation_interval < other_it.approximation_interval;
                    }

                    if (this_it.approximation_interval < other_it.approximation_interval) {
                        return true;
                    }

                    if (other_it.approximation_interval < this_it.approximation_interval) {
                        return false;
                    }
                }

                // If the precision is reached and the number ranges still overlap, then we cannot
                // know if they are equals or other is less than this and we throw an error.
                throw boost::real::precision_exception();
            }

            bool operator>(const real& other) const {
                auto this_it = this->cbegin();
                auto other_it = other.cbegin();

                unsigned int current_precision = std::max(this->max_precision(), other.max_precision());
                for (unsigned int p = 0; p < current_precision; ++p) {
                    // Get more precision
                    ++this_it;
                    ++other_it;

                    bool this_full_precision = this_it.approximation_interval.is_a_number();
                    bool other_full_precision = other_it.approximation_interval.is_a_number();
                    if (this_full_precision && other_full_precision) {
                        return this_it.approximation_interval > other_it.approximation_interval;
                    }

                    if (this_it.approximation_interval > other_it.approximation_interval) {
                        return true;
                    }

                    if (other_it.approximation_interval > this_it.approximation_interval) {
                        return false;
                    }
                }

                // If the precision is reached and the number ranges still overlap, then we cannot
                // know if they are equals or other is less than this and we throw an error.
                throw boost::real::precision_exception();
            }

            bool operator==(const real& other) const {
                auto this_it = this->cbegin();
                auto other_it = other.cbegin();

                unsigned int current_precision = std::max(this->max_precision(), other.max_precision());
                for (unsigned int p = 0; p < current_precision; ++p) {
                    // Get more precision
                    ++this_it;
                    ++other_it;

                    bool this_full_precision = this_it.approximation_interval.is_a_number();
                    bool other_full_precision = other_it.approximation_interval.is_a_number();
                    if (this_full_precision && other_full_precision) {
                        return this_it.approximation_interval == other_it.approximation_interval;
                    }

                    bool this_is_lower = this_it.approximation_interval < other_it.approximation_interval;
                    bool other_is_lower = other_it.approximation_interval < this_it.approximation_interval;
                    if (this_is_lower || other_is_lower) {
                        return false;
                    }
                }

                // If the precision is reached and the numbers full precision is not reached, then
                // we cannot know if they are equals or not.
                throw boost::real::precision_exception();
            }
        };
    }
}

std::ostream& operator<<(std::ostream& os, const boost::real::real& r) {
    os << r.cend().approximation_interval;
    return os;
}

#endif //BOOST_REAL_REAL_HPP