diff --git a/bench/main-bench.cpp b/bench/main-bench.cpp index e137e22..53beda7 100644 --- a/bench/main-bench.cpp +++ b/bench/main-bench.cpp @@ -3,6 +3,10 @@ // ensure this is >= to MAX_NUM_DIGITS_XX for all benchmarks, else we will get // a precision error and the benchmarks will not be meaningful. +<<<<<<< HEAD +std::optional boost::real::const_precision_iterator::maximum_precision = 10000; +======= std::optional boost::real::const_precision_iterator::maximum_precision = 10000; +>>>>>>> upstream/master BENCHMARK_MAIN(); diff --git a/include/real/const_precision_iterator.hpp b/include/real/const_precision_iterator.hpp index 633ce9f..e654035 100644 --- a/include/real/const_precision_iterator.hpp +++ b/include/real/const_precision_iterator.hpp @@ -28,18 +28,17 @@ namespace boost { // fwd decl class real; - // same typedef is also found in real_data.hpp typedef std::variant real_number; - /// the default max precision to use if the user hasn't provided one. - const unsigned int DEFAULT_MAX_PRECISION = 10; + const size_t DEFAULT_MAX_PRECISION = 10; class const_precision_iterator { public: /** * @brief Optional user-provided maximum precision for all const_precision_iterators. */ - static std::optional maximum_precision; + + static std::optional maximum_precision; /// @TODO look into STL-style iterators // typedef std::forward_iterator_tag iterator_category; @@ -57,10 +56,10 @@ namespace boost { real_number * _real_ptr; /// current iterator precision - int _precision; + size_t _precision; /// local max precision, is used if set to > 0 by user - unsigned int _maximum_precision = 0; + size_t _maximum_precision = 0; interval _approximation_interval; @@ -97,7 +96,7 @@ namespace boost { * They may also set a general maximum precision (the static optional value). * Preference is given: _maximum_precision > maximum_precision > DEFAULT_MAX_PRECISION */ - unsigned int max_precision() const { + size_t max_precision() const { if((_maximum_precision == 0) && !(maximum_precision)) return DEFAULT_MAX_PRECISION; else if (_maximum_precision == 0) @@ -106,7 +105,7 @@ namespace boost { return _maximum_precision; } - void set_maximum_precision(unsigned int maximum_precision) { + void set_maximum_precision(size_t maximum_precision) { this->_maximum_precision = maximum_precision; } @@ -145,7 +144,7 @@ namespace boost { if (first_digit == 9) { this->_approximation_interval.upper_bound.digits.push_back(1); this->_approximation_interval.upper_bound.exponent++; - } else if (this->_precision < (int)real.digits().size()) { + } else if (this->_precision < real.digits().size()) { this->_approximation_interval.upper_bound.digits.push_back(first_digit + 1); } else { this->_approximation_interval.upper_bound.digits.push_back(first_digit); @@ -172,8 +171,9 @@ namespace boost { }, [this] (real_operation& real) { - // we don't need to init operands here - they *SHOULD* already be at cbegin or > + // we shouldn't need to init operands here - they *SHOULD* already be at cbegin or > update_operation_boundaries(real); + // _maximum_precision = std::max(real.get_lhs_itr().max_precision(), real.get_rhs_itr().max_precision()); }, [] (auto& real) { throw boost::real::bad_variant_access_exception(); @@ -198,8 +198,8 @@ namespace boost { *this = const_precision_iterator(a); this->iterate_n_times(this->max_precision() - 1); }, - [this] (real_operation& real) { - init_operation_itr(real, true); + [this, cend] (real_operation& real) { + init_operation_itr(real, cend); update_operation_boundaries(real); }, [] (auto & real) { @@ -262,7 +262,7 @@ namespace boost { void iterate_n_times(int n) { std::visit( overloaded { // perform operation on whatever is held in variant [this, &n] (real_explicit& real) { - if (this->_precision >= (int)real.digits().size()) { + if (this->_precision >= real.digits().size()) { return; } @@ -273,14 +273,13 @@ namespace boost { // If the explicit number just reaches the full precision (the end) // then set both boundaries are equals. - if (this->_precision + n >= (int)real.digits().size()) { + if (this->_precision + n >= real.digits().size()) { for(int i = this->_precision; i < (int)real.digits().size(); i++) { this->_approximation_interval.lower_bound.push_back(real.digits()[i]); } this->_approximation_interval.upper_bound = this->_approximation_interval.lower_bound; - } else { // If the explicit number didn't reaches the full precision (the end) @@ -316,7 +315,7 @@ namespace boost { this->_approximation_interval.upper_bound.normalize_left(); this->check_and_swap_boundaries(); - this->_precision = std::min(this->_precision + n, (int)real.digits().size()); + this->_precision = std::min(this->_precision + n, real.digits().size()); }, [this, &n] (real_algorithm& real) { // If the number is negative, bounds are interpreted as mirrored: diff --git a/include/real/exact_number.hpp b/include/real/exact_number.hpp index ffdd9cf..7486177 100644 --- a/include/real/exact_number.hpp +++ b/include/real/exact_number.hpp @@ -4,14 +4,20 @@ #include #include #include +#include +#include namespace boost { namespace real { struct exact_number { + using exponent_t = int; + std::vector digits = {}; - int exponent = 0; + exponent_t exponent = 0; bool positive = true; + static const int DIGIT_MAX = 9; + static const int BASE = 10; static bool aligned_vectors_is_lower(const std::vector &lhs, const std::vector &rhs) { // Check if lhs is lower than rhs @@ -171,6 +177,244 @@ namespace boost { this->normalize(); } + /// calculates *this / divisor + /// + /// @brief a binary-search type method for dividing exact_numbers. + void divide_vector(exact_number divisor, unsigned int max_precision) { + /// @TODO: replace this with something more efficient, + // it also completely recalculates on each precision increase + // instead, could use previous information to make better "guesses" + // for our iteration scheme. + + boost::real::exact_number numerator; + boost::real::exact_number left; + boost::real::exact_number right; + boost::real::exact_number residual; + boost::real::exact_number tmp; + boost::real::exact_number half ({5}, 0); + boost::real::exact_number distance; + ///@TODO ensure exponent doesn't overflow + boost::real::exact_number min_boundary_n({DIGIT_MAX}, -1 * max_precision, false); + boost::real::exact_number min_boundary_p({DIGIT_MAX}, -1 * max_precision); + + bool positive = ((*this).positive == divisor.positive); + numerator = (*this).abs(); + divisor = divisor.abs(); + + // we ignore exponents, then set them in the end. + // (a * 10^x) / (b*10^y) = (a*10 / b*10) * 10^((x-1)-(y-1)) + // 100/20 -> 1/2 * (10) + int exponent_dif = (numerator.exponent - 1) - (divisor.exponent - 1); + + numerator.exponent = 1; + divisor.exponent = 1; + + // std::cout << "num/den" << numerator.as_string() << ", " << divisor.as_string() << '\n'; + + tmp.digits = {1}; + tmp.exponent = 1; + + if (divisor == tmp) { + this->exponent = exponent_dif + 1; + (*this).positive = positive; + return; + } + + if (divisor == (*this)) { + (*this) = tmp; + (*this).positive = positive; + return; + } + + exact_number zero = exact_number(); + + if (divisor == zero) { + throw(boost::real::divide_by_zero()); + } + + if ((*this) == zero) { + return; + } + + // N < D --> 0 < abs(Q) < 1 + if (numerator < divisor) { + left = exact_number(); // 0 + right = tmp; // 1 + } else { // assuming D > 1. N > D ---> 1 < N / D < N + left = tmp; // 1 + right = numerator; + } + + // distance = (right - left) / 2 + distance = (right - left) * half; + (*this) = left + distance; + // N/D = Q -> QD -N = 0 + residual = (*this) * divisor - numerator; + if (residual == zero) { + this->exponent += exponent_dif; + this->positive = positive; + return; + } + + // calculate the result + // continue the loop while we are still inaccurate (up to max precision) + + // we would prefer to use || than &&, but we get problems because || + // statement may never evaluate to false. Again, we must look more closely at precision + // if we use && (to always terminate in n steps), then we must + // verify our answers are within +- epsilon of the solution. + while ((residual.abs() > min_boundary_p) && + (distance.exponent >= min_boundary_p.exponent)) { + // result too small, try halfway between left and (*this) + if (residual < min_boundary_n) { + left = (*this); + } + + // distance is halved + distance = half * distance; + distance.normalize(); + + // truncate insignificant digits of distance + while (distance.size() > max_precision + 1) { + distance.digits.pop_back(); + } + + // iterate (*this) + (*this) = left + distance; + + // truncate insignificant digits of (*this) + while ((*this).size() > max_precision + 1) { + (*this).digits.pop_back(); + } + + // recalculate residual N/D = Q ---> QD - N = residual + residual = ((*this) * divisor) - numerator; + residual.normalize(); + } // end while + // at this point, we may have exited early. we should be within +- min_boundary_p of the solution. + + // truncate (*this) + this->normalize(); + while ((*this).size() > max_precision) { + (*this).digits.pop_back(); + } + + // recalculate residual for the final (*this) value + exact_number residual_o = ((*this) * divisor) - numerator; + + // note we have to normalize before comparison, because -0.0 != zero .. + residual_o.normalize(); + + if (residual_o != zero) { // then, we are not fully accurate + // we try seeing if we can make the residual equal zero by adding/subtracting epsilon + exact_number tmp_lower = (*this); + tmp_lower.round_down(); + exact_number tmp_upper = (*this); + tmp_upper.round_up(); + + residual = tmp_lower * divisor - numerator; + residual.normalize(); + + if (residual == zero) { // rounded lower is better + (*this) = tmp_lower; + + if (positive) { + (*this).positive = true; + } + else { + (*this).positive = false; + } + + (*this).normalize(); + this->exponent += exponent_dif; + return; + } + + residual = tmp_upper * divisor - numerator; + residual.normalize(); + + if (residual == zero) { // rounded up is better + (*this) = tmp_upper; + + if (positive) { + (*this).positive = true; + } + else { + (*this).positive = false; + } + + (*this).normalize(); + this->exponent += exponent_dif; + return; + } + // at this point, it is impossible to make the residual 0 + } // else old residual == 0 + + if (positive) { + (*this).positive = true; + } + else { + (*this).positive = false; + } + + this->exponent += exponent_dif; + (*this).normalize(); + } + + + void round_up_abs() { + int index = digits.size() - 1; + bool keep_carrying = true; + + while((index > 0) && keep_carrying) { // bring the carry back + if(this->digits[index] != DIGIT_MAX) { + ++this->digits[index]; + keep_carrying = false; + } else { // digits[index] == 9, we keep carrying + this->digits[index] = 0; + } + --index; + } + + if ((index == 0) && keep_carrying) { // i.e., .999 should become 1.000 + if(this->digits[index] == DIGIT_MAX) { + this->digits[index] = 0; + this->push_front(1); + ++this->exponent; + } else { + ++this->digits[index]; + } + } + } + + void round_up() { + if (positive) + this->round_up_abs(); + } + + void round_down() { + if (!positive) + this->round_up_abs(); + } + + void round_down_abs() { + int index = digits.size() - 1; + bool keep_carrying = true; + + while((index > 0) && keep_carrying) { // bring the carry back + if(this->digits[index] != 0) { + --this->digits[index]; + keep_carrying = false; + } else { // digits[index] == 0, we keep carrying + this->digits[index] = DIGIT_MAX; + } + --index; + } + // we should be ok at this point because the first number in digits should != 0 + if (index == 0 && keep_carrying) + --this->digits[index]; + } + /** * @brief *default constructor*: It constructs a representation of the number zero. */ @@ -179,6 +423,74 @@ namespace boost { /// ctor from vector of digits, integer exponent, and optional bool positive exact_number(std::vector vec, int exp, bool pos = true) : digits(vec), exponent(exp), positive(pos) {}; + /// ctor from any integral type + template::value>> + exact_number(I x) { + if (x < 0) + positive = false; + else + positive = true; + + exponent = 0; + while (((x % BASE) != 0) || (x != 0)) { + exponent++; + push_front(std::abs(x%BASE)); + x /= BASE; + } + } + + explicit exact_number (const std::string& number) { + std::regex decimal("((\\+|-)?[[:digit:]]*)(\\.(([[:digit:]]+)?))?((e|E)(((\\+|-)?)[[:digit:]]+))?"); + if (!std::regex_match (number, decimal)) + throw boost::real::invalid_string_number_exception(); + //Know at this point that representation is valid + std::string decimal_part = regex_replace(number, decimal, "$5"); + std::string integer_part = regex_replace(number, decimal, "$1"); + std::string exp = regex_replace(number, decimal, "$8"); + int add_exponent = exp.length() == 0 ? 0 : std::stoi(exp); + if (integer_part[0] == '+') { + positive = true; + integer_part = integer_part.substr(1); + } + else if (integer_part[0] == '-') { + positive = false; + integer_part = integer_part.substr(1); + } + integer_part = regex_replace(integer_part, std::regex("(0?+)([[:digit:]]?+)"), "$2"); + size_t i = decimal_part.length() - 1; + while (decimal_part[i] == '0' && i > 0) { + --i; + } + decimal_part = decimal_part.substr(0, i + 1); + //decimal and integer parts are stripped of zeroes + exponent = integer_part.length() + add_exponent; + if (decimal_part.empty()) { + i = integer_part.length() - 1; + while (integer_part[i] == '0' && i > 0) + --i; + integer_part = integer_part.substr(0, i + 1); + } + if (integer_part.empty()) { + i = 0; + while (decimal_part[i] == '0' && i < decimal_part.length()) { + ++i; + --exponent; + } + decimal_part = decimal_part.substr(i); + } + if (integer_part.empty() && decimal_part.empty()) { + digits = {0}; + exponent = 0; + return; + } + for (const auto& c : integer_part ) { + digits.push_back(c - '0'); + } + for (const auto& c : decimal_part ) { + digits.push_back(c - '0'); + } + } + /** * @brief *Copy constructor:* It constructs a new boost::real::exact_number that is a copy of the * other boost::real::exact_number. @@ -205,10 +517,10 @@ namespace boost { */ bool operator<(const exact_number& other) const { std::vector zero = {0}; - if (this->digits == zero) { - if (other.digits == zero || !other.positive) - return false; - else + if (this->digits == zero || this->digits.empty()) { + return !(other.digits == zero || !other.positive || other.digits.empty()); + } else { + if ((other.digits == zero || other.digits.empty()) && !this->positive) return true; } @@ -240,10 +552,10 @@ namespace boost { */ bool operator>(const exact_number& other) const { std::vector zero = {0}; - if (this->digits == zero) { - if (other.digits == zero || other.positive) - return false; - else + if (this->digits == zero || this->digits.empty()) { + return !(other.digits == zero || other.positive || other.digits.empty()); + } else { + if ((other.digits == zero || other.digits.empty()) && this->positive) return true; } @@ -266,6 +578,14 @@ namespace boost { return other.exponent > this->exponent; } + bool operator>=(const exact_number& other) const { + return (!((*this) < other)); + } + + bool operator<=(const exact_number& other) const { + return (!((*this) > other)); + } + /** * @brief *Equality comparator operator:* It compares the *this boost::real::exact_number with the other * boost::real::exact_number to determine if *this and other are equals or not. @@ -277,8 +597,12 @@ namespace boost { return !(*this < other || other < *this); } - static exact_number abs(const exact_number& b) { - exact_number result = b; + bool operator!=(const exact_number& other) const { + return !(*this == other); + } + + exact_number abs() { + exact_number result = (*this); result.positive = true; return result; } @@ -290,7 +614,7 @@ namespace boost { if (positive == other.positive) { result = *this; result.add_vector(other); - } else if (abs(other) < abs(*this)) { + } else if (other.abs() < this->abs()) { result = *this; result.subtract_vector(other); } else { @@ -311,7 +635,7 @@ namespace boost { result = *this; result.add_vector(other); } else { - if (abs(other) < abs(*this)) { + if (other.abs() < this->abs()) { result = *this; result.subtract_vector(other); } else { @@ -334,7 +658,7 @@ namespace boost { return result; } - void operator*=(exact_number &other) { + void operator*=(exact_number other) { *this = *this * other; } @@ -373,8 +697,8 @@ namespace boost { } } } else { - int digit_amount = std::max(this->exponent, (int) this->digits.size()); + for (int i = 0; i < digit_amount; ++i) { if (i == this->exponent) { @@ -392,8 +716,6 @@ namespace boost { result.pop_back(); } } - - return result; } @@ -487,6 +809,14 @@ namespace boost { return this->digits.size(); } + bool is_integral() { + if (exponent < 0) { + return false; + } else { + /// @TODO: ensure it doesn't overflow + return (exponent_t) digits.size() <= exponent; + } + } }; } } diff --git a/include/real/irrational_helpers.hpp b/include/real/irrational_helpers.hpp index 38fff84..185ef59 100644 --- a/include/real/irrational_helpers.hpp +++ b/include/real/irrational_helpers.hpp @@ -3,10 +3,14 @@ #include #include +#include +#include namespace boost { namespace real { namespace irrational { + // used for extra precision. should be replaced with something more definitive in the future. + inline const int PLACEHOLDER = 10; /** * @brief The function returns the n-th digit of the champernowne number in the @@ -46,6 +50,56 @@ namespace boost { return binary[binary.size() - 1 - (index - n)]; } + + /// @TODO: figure out how to avoid unnecessary recalculation by saving + // some information in real_algorithm + + // e^x = sum_{k=0}^\infty x^k / k! = 1 + x + x^2/2! + x^3/3! + x^4/4! + // calculates e^x, where x = m/l + template + int exponential_get_nth_digit(unsigned int n) { + std::vector const one {1}; + std::vector const zero {0}; + + exact_number one_n {one, 1}; + exact_number last_term {one, 1}; // x^k_n / k_n! + exact_number current_value; // summation up to k_n + + exact_number k; + exact_number k_fac; + exact_number const mn (m); + exact_number const ln (l); + std::cout << mn.as_string() << ", " << ln.as_string() << '\n'; + exact_number x = mn; + + /// @TODO look into possible precision problems + x.divide_vector(ln, n+PLACEHOLDER); + + // prepare to calculate k=2 term + k = one_n + one_n; + last_term = x; + current_value = one_n + x; + + // keep getting terms from the taylor series until the terms go below our precision bound + /// @TODO: ensure cast doesn't overflow + while(last_term.exponent >= 1-(int)(n+PLACEHOLDER)) { + last_term *= x; + last_term.divide_vector(k, n+PLACEHOLDER); + current_value += last_term; + + k = k + one_n; + } + + return current_value.digits[n]; + } + + int pi_get_nth_digit(unsigned int n) { + if (n > 2000) + throw(pi_precision_exception()); + // 2000 digits of pi + exact_number pi("3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809525720106548586327886593615338182796823030195203530185296899577362259941389124972177528347913151557485724245415069595082953311686172785588907509838175463746493931925506040092770167113900984882401285836160356370766010471018194295559619894676783744944825537977472684710404753464620804668425906949129331367702898915210475216205696602405803815019351125338243003558764024749647326391419927260426992279678235478163600934172164121992458631503028618297455570674983850549458858692699569092721079750930295532116534498720275596023648066549911988183479775356636980742654252786255181841757467289097777279380008164706001614524919217321721477235014144197356854816136115735255213347574184946843852332390739414333454776241686251898356948556209921922218427255025425688767179049460165346680498862723279178608578438382796797668145410095388378636095068006422512520511739298489608412848862694560424196528502221066118630674427862203919494504712371378696095636437191728746776465757396241389086583264599581339047802759009"); + return pi[n]; + } } } } diff --git a/include/real/real.hpp b/include/real/real.hpp index e10d4fb..910f2e1 100644 --- a/include/real/real.hpp +++ b/include/real/real.hpp @@ -55,12 +55,13 @@ namespace boost { * operator "==" but for those cases where the class is not able to decide the value of the * result before reaching the maximum precision, a precision_exception is thrown. */ - - class real { private: std::shared_ptr _real_p; + + // ctor from shared_ptr to (already init) real_data. used in check_and_distribute. + real(std::shared_ptr x) : _real_p(x){}; public: @@ -91,7 +92,6 @@ namespace boost { */ real(const std::string& number) : _real_p(std::make_shared(real_explicit(number))) {}; - /** * @brief Initializer list constructor @@ -102,7 +102,6 @@ namespace boost { : _real_p(std::make_shared(real_explicit(digits, digits.size()))) {}; - /** * @brief *Signed initializer list constructor:* Creates a boost::real::real * instance that represents the number where the positive parameter is used to set the @@ -189,7 +188,7 @@ namespace boost { return _real_p->get_real_number(); } - const_precision_iterator get_real_itr() const { + const_precision_iterator& get_real_itr() const { return _real_p->get_precision_itr(); } @@ -239,6 +238,191 @@ namespace boost { return ret; } + /// @TODO: do we need different ctors to be more efficient? rvalue AND lvalue ref? + + // a constant used in the print tree helper function + static const int PRINT_SPACE = 5; + // a helper function for observing the trees within a real + // note the tree is displayed with a left 90 degree rotation + // (i.e., root on the left, nodes on the next leftmost level) + void print_tree(int space = PRINT_SPACE) { + std::visit(overloaded { + [space] (const real_explicit& real) { + for (int i = PRINT_SPACE; i < space; i++) + std::cout << ' '; + std::cout << real.get_exact_number().as_string() << '\n'; + }, + [space] (const real_algorithm& real) { + for (int i = PRINT_SPACE; i < space; i++) + std::cout << ' '; + std::cout << "alg\n"; + }, + [&space] (const real_operation& real) { + ((boost::real::real) real.rhs()).print_tree(space + PRINT_SPACE); + std::cout << '\n'; + + for (int i = PRINT_SPACE; i < space; i++) + std::cout << ' '; + + switch(real.get_operation()) { + case OPERATION::ADDITION: + std::cout << "+"; + break; + case OPERATION::SUBTRACTION: + std::cout << "-"; + break; + case OPERATION::MULTIPLICATION: + std::cout << "*"; + break; + case OPERATION::DIVISION: + std::cout << "/"; + break; + } + std::cout << '\n'; + + ((boost::real::real) real.lhs()).print_tree(space + PRINT_SPACE); + }, + [] (auto& real) { + throw boost::real::bad_variant_access_exception(); + } + }, _real_p->get_real_number()); + } + + // a helper function for distributing when performing addition/subtraction + // the returned bool tells us whether we distributed or not, which is mostly useful when assign_and_return_void is true + // since std::optional would = std::null_opt if assign_and_return_void is false. + std::pair> check_and_distribute(real & other, bool assign_and_return_void, OPERATION op) { + // The following simplifies using the distributive property, when numbers have the same pointers. + // We could do comparison by value, but this may force more computation than is necessary for the user, + // since it's difficult to determine whether values are the same + + std::shared_ptr a; + std::shared_ptr b; + std::shared_ptr x; + + if(auto op_ptr = std::get_if(this->_real_p->get_real_ptr())) { + if (auto op_ptr2 = std::get_if(other._real_p->get_real_ptr())) { // both of real_operation + if ((op_ptr->get_operation() == OPERATION::MULTIPLICATION) && op_ptr2->get_operation() == OPERATION::MULTIPLICATION) { + if (op_ptr->lhs() == op_ptr2->lhs()) { // x * a + x * b = (a + b) * x + a = op_ptr->rhs(); + b = op_ptr2->rhs(); + x = op_ptr->lhs(); + + } else if (op_ptr->lhs() == op_ptr2->rhs()) { // x * a + b * x = (a + b) * x + a = op_ptr->rhs(); + b = op_ptr2->lhs(); + x = op_ptr->lhs(); + + } else if (op_ptr->rhs() == op_ptr2->lhs()) { // a * x + x * b = (a + b) * x + a = op_ptr->lhs(); + b = op_ptr2->rhs(); + x = op_ptr->rhs(); + + } else if (op_ptr->rhs() == op_ptr2->rhs()) { // a * x + b * x = (a + b) * x + a = op_ptr->lhs(); + b = op_ptr2->lhs(); + x = op_ptr->rhs(); + } else { + return std::make_pair(false, std::nullopt); + } + + real a_op_b; + if(op == OPERATION::ADDITION) { + a_op_b = real(a) + real(b); + } + else if (op == OPERATION::SUBTRACTION) { + a_op_b = real(a) + real(b); + } else { + throw invalid_distribution_operation_exception(); + } + + if(assign_and_return_void) { + this->_real_p = std::make_shared(real_operation(a_op_b._real_p, x, OPERATION::MULTIPLICATION)); + return std::make_pair(true, std::nullopt); + } else { + return std::make_pair(true, real(real_operation(a_op_b._real_p, x, OPERATION::MULTIPLICATION))); + } + } + } else { // rhs is not an operation + if (op_ptr->get_operation() == OPERATION::MULTIPLICATION) { + if (other._real_p == op_ptr->lhs()) { // (a * x) + a -> (x + 1) * a + a = op_ptr->lhs(); + x = op_ptr->rhs(); + + } else if (other._real_p == op_ptr->rhs()) {// (x * a) + a -> (x + 1) * a + a = op_ptr->rhs(); + x = op_ptr->lhs(); + } else { + return std::make_pair(false, std::nullopt); + } + + real one ("1"); + real x_op_1; + if(op == OPERATION::ADDITION) { + x_op_1 = real(x) + one; + } + else if (op == OPERATION::SUBTRACTION) { + x_op_1 = real(x) + one; + } else { + throw invalid_distribution_operation_exception(); + } + + if(assign_and_return_void) { + this->_real_p = std::make_shared(real_operation(x_op_1._real_p, a, OPERATION::MULTIPLICATION)); + return std::make_pair(true, std::nullopt); + } else { + return std::make_pair(true, real(real_operation(x_op_1._real_p, a, OPERATION::MULTIPLICATION))); + } + } + } + } else if(auto op_ptr = std::get_if(&other._real_p->get_real_number())) { // lhs is not an operation, but rhs is + if (op_ptr->get_operation() == OPERATION::MULTIPLICATION) { + if (this->_real_p == op_ptr->lhs()) { // a + (a * x) -> (x + 1) * a + a = this->_real_p; + x = op_ptr->rhs(); + + } else if (this->_real_p == op_ptr->rhs()) {// a + (x * a) -> (x + 1) * a + a = this->_real_p; + x = op_ptr->lhs(); + } else { + return std::make_pair(false, std::nullopt); + } + + real x_op_1; + real one ("1"); + if(op == OPERATION::ADDITION) { + x_op_1 = real(x) + one; + } + else if (op == OPERATION::SUBTRACTION) { + x_op_1 = real(x) + one; + } else { + throw invalid_distribution_operation_exception(); + } + + if(assign_and_return_void) { + this->_real_p = std::make_shared(real_operation(x_op_1._real_p, a, OPERATION::MULTIPLICATION)); + return std::make_pair(true, std::nullopt); + } else { + return std::make_pair(true, real(real_operation(x_op_1._real_p, a, OPERATION::MULTIPLICATION))); + } + } + } else { // neither is an operation + if ((this->_real_p == other._real_p) && (op == OPERATION::ADDITION)) { // a + a = 2 * a + std::shared_ptr two = std::make_shared(real_explicit("2")); + + if(assign_and_return_void) { + this->_real_p = std::make_shared(real_operation(two, this->_real_p, OPERATION::MULTIPLICATION)); + return std::make_pair(true, std::nullopt); + } else { + return std::make_pair(true, real(real_operation(two, this->_real_p, OPERATION::MULTIPLICATION))); + } + } + } + + // at this point, we cannot distribute. + return std::make_pair(false, std::nullopt); + } + /** * @brief Sets this real_data to that of the operation between this previous * real_data and other real_data. @@ -247,14 +431,13 @@ namespace boost { * * @param other - the right side operand boost::real::real number. */ - void operator+=(real& other) { - // do not want to overwrite *this->_real_p if others are pointing to it - // if others are pointing to it, point to a copy in newly allocated memory, then create operation - if(this->_real_p.use_count() > 1) { - this->_real_p = std::make_shared(real_data(*this->_real_p)); + void operator+=(real other) { + auto [is_simplified, result] = check_and_distribute(other, true, OPERATION::ADDITION); + + if (!is_simplified) { + this->_real_p = + std::make_shared(real_operation(this->_real_p, other._real_p, OPERATION::ADDITION)); } - this->_real_p = - std::make_shared(real_operation(this->_real_p, other._real_p, OPERATION::ADDITION)); } /** @@ -265,7 +448,12 @@ namespace boost { * @return A copy of the new boost::real::real number representation. */ real operator+(real other) { - return real(real_operation(this->_real_p, other._real_p, OPERATION::ADDITION)); + auto [is_simplified, result] = check_and_distribute(other, false, OPERATION::ADDITION); + if (is_simplified) { + return result.value(); + } else { + return real(real_operation(this->_real_p, other._real_p, OPERATION::ADDITION)); + } } /** @@ -274,13 +462,13 @@ namespace boost { * * @param other - the right side operand boost::real::real number. */ - void operator-=(real& other) { - if(this->_real_p.use_count() > 1) { - // if others are pointing to it, point to a copy in newly allocated memory, then create operation - this->_real_p = std::make_shared(real_data(*this->_real_p)); + void operator-=(real other) { + auto [is_simplified, result] = check_and_distribute(other, true, OPERATION::SUBTRACTION); + + if(!is_simplified) { + this->_real_p = + std::make_shared(real_operation(this->_real_p, other._real_p, OPERATION::SUBTRACTION)); } - this->_real_p = - std::make_shared(real_operation(this->_real_p, other._real_p, OPERATION::SUBTRACTION)); } /** @@ -291,7 +479,12 @@ namespace boost { * @return A copy of the new boost::real::real number representation. */ real operator-(real other) { - return real(real_operation(this->_real_p, other._real_p, OPERATION::SUBTRACTION)); + auto [is_simplified, result] = check_and_distribute(other, false, OPERATION::SUBTRACTION); + if (is_simplified) { + return result.value(); + } else { + return real(real_operation(this->_real_p, other._real_p, OPERATION::SUBTRACTION)); + } } /** @@ -300,10 +493,7 @@ namespace boost { * * @param other - the right side operand boost::real::real number. */ - void operator*=(real& other) { - if(this->_real_p.use_count() > 1) { - this->_real_p = std::make_shared(real_data(*this->_real_p)); - } + void operator*=(real other) { this->_real_p = std::make_shared(real_operation(this->_real_p, other._real_p, OPERATION::MULTIPLICATION)); } @@ -319,19 +509,34 @@ namespace boost { return real(real_operation(this->_real_p, other._real_p, OPERATION::MULTIPLICATION)); } + /** + * @brief Sets this real_data to that of the operation between + * this previous real_data and other real_data. + * + * @param other - the right side operand boost::real::real number. + */ + void operator/=(real other) { + this->_real_p = + std::make_shared(real_operation(this->_real_p, other._real_p, OPERATION::DIVISION)); + } + + /** + * @brief Creates a new boost::real::real representing the product + * of *this and other + * + * @param other - the right side operand boost::real::real number. + * @return A copy of the new boost::real::real number representation. + */ + real operator/(real other) { + return real(real_operation(this->_real_p, other._real_p, OPERATION::DIVISION)); + } + /** * @brief Assigns *this to other * @param other - the boost::real::real number to copy. */ - void operator=(real& other) { - if(this->_real_p.use_count() > 1) { - // if this is being referenced to, point to new memory before assigning - // i.e., if A = B + B, and we do B = D, we first make B point elsewhere. - // so that A != D + D - this->_real_p = std::make_shared(); - } - this->_real_p = - std::make_shared(*other._real_p); + void operator=(real other) { + this->_real_p = other._real_p; } /** @@ -339,10 +544,6 @@ namespace boost { * @param number - a valid string representing a number. */ void operator=(const std::string& number) { - if(this->_real_p.use_count() > 1) { - // if this is being referenced to, point to new memory before assigning - this->_real_p = std::make_shared(); - } this->_real_p = std::make_shared(real_explicit(number)); } @@ -505,6 +706,4 @@ inline boost::real::real operator "" _r(const char* x, size_t len) { return boost::real::real(x); } - - #endif //BOOST_REAL_HPP diff --git a/include/real/real_algorithm.hpp b/include/real/real_algorithm.hpp index 4b5a324..48e305e 100644 --- a/include/real/real_algorithm.hpp +++ b/include/real/real_algorithm.hpp @@ -106,4 +106,4 @@ namespace boost { } -#endif //BOOST_REAL_REAL_ALGORITHM_HPP +#endif //BOOST_REAL_REAL_ALGORITHM_HPP \ No newline at end of file diff --git a/include/real/real_data.hpp b/include/real/real_data.hpp index 16d73b1..b5c1366 100644 --- a/include/real/real_data.hpp +++ b/include/real/real_data.hpp @@ -36,6 +36,10 @@ namespace boost { return _real; } + real_number const * get_real_ptr() const { + return &_real; + } + const_precision_iterator& get_precision_itr() { return _precision_itr; } @@ -58,16 +62,15 @@ namespace boost { case OPERATION::SUBTRACTION: - this->_approximation_interval.lower_bound = - ro.get_lhs_itr().get_interval().lower_bound - - ro.get_rhs_itr().get_interval().upper_bound; + this->_approximation_interval.lower_bound = + ro.get_lhs_itr().get_interval().lower_bound - + ro.get_rhs_itr().get_interval().upper_bound; - this->_approximation_interval.upper_bound = - ro.get_lhs_itr().get_interval().upper_bound - - ro.get_rhs_itr().get_interval().lower_bound; + this->_approximation_interval.upper_bound = + ro.get_lhs_itr().get_interval().upper_bound - + ro.get_rhs_itr().get_interval().lower_bound; break; - case OPERATION::MULTIPLICATION: { bool lhs_positive = ro.get_lhs_itr().get_interval().positive(); bool rhs_positive = ro.get_rhs_itr().get_interval().positive(); @@ -162,7 +165,124 @@ namespace boost { } break; } + case OPERATION::DIVISION: { + exact_number zero = exact_number(); + exact_number residual; + exact_number quotient; + exact_number numerator; + exact_number denominator; + + // if the interval contains zero, iterate until it doesn't, or until max_precision. + while (!ro.get_rhs_itr().get_interval().positive() && + !ro.get_rhs_itr().get_interval().positive() && + _precision <= this->max_precision()) + ++(*this); + + // if the interval contains zero after iterating until max precision, throw, + // because this causes one side of the result interval to tend towards +/-infinity + if (!ro.get_rhs_itr().get_interval().positive() && + !ro.get_rhs_itr().get_interval().negative()) + throw boost::real::divergent_division_result_exception(); + + // Q = N/D + // first, the upper boundary + if (ro.get_lhs_itr().get_interval().positive()) { + if (ro.get_rhs_itr().get_interval().positive()) { + numerator = ro.get_lhs_itr().get_interval().upper_bound; + denominator = ro.get_rhs_itr().get_interval().lower_bound; + } else { + numerator = ro.get_lhs_itr().get_interval().lower_bound; + denominator = ro.get_rhs_itr().get_interval().upper_bound; + } + } else if (ro.get_lhs_itr().get_interval().negative()) { + if (ro.get_rhs_itr().get_interval().positive()) { + numerator = ro.get_lhs_itr().get_interval().upper_bound; + denominator = ro.get_rhs_itr().get_interval().lower_bound; + } else if (ro.get_rhs_itr().get_interval().negative()) { + numerator = ro.get_lhs_itr().get_interval().lower_bound; + denominator = ro.get_rhs_itr().get_interval().upper_bound; + } + } else { + if (ro.get_rhs_itr().get_interval().positive()) { + numerator = ro.get_lhs_itr().get_interval().upper_bound; + denominator = ro.get_rhs_itr().get_interval().upper_bound; + } else if (ro.get_rhs_itr().get_interval().negative()) { + numerator = ro.get_lhs_itr().get_interval().lower_bound; + denominator = ro.get_rhs_itr().get_interval().lower_bound; + } + } + + // calculate the upper bound + quotient = numerator; + quotient.divide_vector(denominator, this->max_precision()); + + residual = quotient * denominator - numerator; + residual.normalize(); + quotient.normalize(); + + this->_approximation_interval.upper_bound = quotient; + + if (residual.abs() > zero) { + this->_approximation_interval.upper_bound.round_up(); + } + + // if both operands are numbers (not intervals), then we can skip doing the lower bound separately + if (ro.get_rhs_itr().get_interval().is_a_number() && ro.get_lhs_itr().get_interval().is_a_number()) { + _approximation_interval.lower_bound = quotient; + + if (residual == zero) { + _approximation_interval.upper_bound = _approximation_interval.lower_bound; + } else { + _approximation_interval.lower_bound.round_down(); + } + return; + } + // lower boundary + if (ro.get_lhs_itr().get_interval().positive()) { + if (ro.get_rhs_itr().get_interval().positive()) { + numerator = ro.get_lhs_itr().get_interval().lower_bound; + denominator = ro.get_rhs_itr().get_interval().upper_bound; + } else { + numerator = ro.get_lhs_itr().get_interval().upper_bound; + denominator = ro.get_rhs_itr().get_interval().lower_bound; + } + } else if (ro.get_lhs_itr().get_interval().negative()) { + if (ro.get_rhs_itr().get_interval().positive()) { + numerator = ro.get_lhs_itr().get_interval().lower_bound; + denominator = ro.get_rhs_itr().get_interval().upper_bound; + } else if (ro.get_rhs_itr().get_interval().negative()) { + numerator = ro.get_lhs_itr().get_interval().upper_bound; + denominator = ro.get_rhs_itr().get_interval().lower_bound; + } + } else { + if (ro.get_rhs_itr().get_interval().positive()) { + numerator = ro.get_lhs_itr().get_interval().lower_bound; + denominator = ro.get_rhs_itr().get_interval().lower_bound; + } else if (ro.get_rhs_itr().get_interval().negative()) { + numerator = ro.get_lhs_itr().get_interval().upper_bound; + denominator = ro.get_rhs_itr().get_interval().upper_bound; + } + } + + quotient = _approximation_interval.lower_bound; + + quotient = numerator; + quotient.divide_vector(denominator, this->max_precision()); + + residual = quotient * denominator - numerator; + residual.normalize(); + quotient.normalize(); + + this->_approximation_interval.lower_bound = quotient; + + if (residual.abs() > zero) { + this->_approximation_interval.lower_bound.round_down(); + } + + break; + } + default: throw boost::real::none_operation_exception(); } @@ -184,11 +304,13 @@ namespace boost { inline void const_precision_iterator::operation_iterate_n_times(real_operation &ro, int n) { /// @warning there could be issues if operands have different precisions/max precisions - if (ro.get_lhs_itr()._precision < this->_precision + n) + if (ro.get_lhs_itr()._precision < this->_precision + n) { ro.get_lhs_itr().iterate_n_times(n); + } - if (ro.get_rhs_itr()._precision < this->_precision + n) + if (ro.get_rhs_itr()._precision < this->_precision + n) { ro.get_rhs_itr().iterate_n_times(n); + } this->_precision += n; diff --git a/include/real/real_exception.hpp b/include/real/real_exception.hpp index 527c9f7..2656924 100644 --- a/include/real/real_exception.hpp +++ b/include/real/real_exception.hpp @@ -39,6 +39,38 @@ namespace boost { return "Cannot perform this method on this real variant type"; } }; + + struct divide_by_zero : public std::exception { + + const char * what () const throw () override { + return "Divison by zero is undefined"; + } + }; + + struct invalid_denominator : public std::exception { + + const char * what () const throw () override { + return "Divison with denominators 0 < d < 1 is currently undefined."; + } + }; + + struct divergent_division_result_exception : public std::exception { + const char * what () const throw () override { + return "The divisor approximation interval contains 0, so the quotient is unbounded"; + } + }; + + struct pi_precision_exception : public std::exception { + const char * what () const throw () override { + return "pi is currently undefined for precision > 2000 digits"; + } + }; + + struct invalid_distribution_operation_exception : public std::exception { + const char * what () const throw () override { + return "Distribution does not work for this operation."; + } + }; } } diff --git a/include/real/real_explicit.hpp b/include/real/real_explicit.hpp index 798e137..b745d58 100644 --- a/include/real/real_explicit.hpp +++ b/include/real/real_explicit.hpp @@ -26,7 +26,7 @@ namespace boost { // Number representation as a vector of digits with an integer part and a sign (+/-) // TODO: Add normalizations to the constructors exact_number explicit_number; - + public: /** @@ -51,58 +51,8 @@ namespace boost { * * @throws boost::real::invalid_string_number exception */ - explicit real_explicit(const std::string& number) { - std::regex decimal("((\\+|-)?[[:digit:]]*)(\\.(([[:digit:]]+)?))?((e|E)(((\\+|-)?)[[:digit:]]+))?"); - if (!std::regex_match (number, decimal)) - throw boost::real::invalid_string_number_exception(); - //Know at this point that representation is valid - std::string decimal_part = regex_replace(number, decimal, "$5"); - std::string integer_part = regex_replace(number, decimal, "$1"); - std::string exp = regex_replace(number, decimal, "$8"); - int add_exponent = exp.length() == 0 ? 0 : std::stoi(exp); - if (integer_part[0] == '+') { - explicit_number.positive = true; - integer_part = integer_part.substr(1); - } - else if (integer_part[0] == '-') { - explicit_number.positive = false; - integer_part = integer_part.substr(1); - } - integer_part = regex_replace(integer_part, std::regex("(0?+)([[:digit:]]?+)"), "$2"); - size_t i = decimal_part.length() - 1; - while (decimal_part[i] == '0' && i > 0) { - --i; - } - decimal_part = decimal_part.substr(0, i + 1); - //decimal and integer parts are stripped of zeroes - int exponent = integer_part.length() + add_exponent; - if (decimal_part.empty()) { - i = integer_part.length() - 1; - while (integer_part[i] == '0' && i > 0) - --i; - integer_part = integer_part.substr(0, i + 1); - } - if (integer_part.empty()) { - i = 0; - while (decimal_part[i] == '0' && i < decimal_part.length()) { - ++i; - --exponent; - } - decimal_part = decimal_part.substr(i); - } - if (integer_part.empty() && decimal_part.empty()) { - explicit_number.digits = {0}; - explicit_number.exponent = 0; - return; - } - explicit_number.exponent = exponent; - for (const auto& c : integer_part ) { - explicit_number.digits.push_back(c - '0'); - } - for (const auto& c : decimal_part ) { - explicit_number.digits.push_back(c - '0'); - } - } + explicit real_explicit(const std::string& number) : explicit_number(number) + {}; /** * @brief *Initializer list constructor with exponent:* Creates a boost::real::real_explicit @@ -153,6 +103,10 @@ namespace boost { return explicit_number.digits; } + exact_number get_exact_number() const { + return explicit_number; + } + /** * @brief Constructs a new boost::real::real_explicit::const_precision_iterator that iterates the number * approximation intervals in increasing order according to the approximation precision. @@ -172,7 +126,6 @@ namespace boost { if (n < explicit_number.digits.size()) { return explicit_number.digits.at(n); } - return 0; } @@ -183,6 +136,8 @@ namespace boost { * @return a reference of *this with the new represented number. */ real_explicit& operator=(const real_explicit& other) = default; + + }; } } diff --git a/include/real/real_operation.hpp b/include/real/real_operation.hpp index 7efffcd..048e89f 100644 --- a/include/real/real_operation.hpp +++ b/include/real/real_operation.hpp @@ -19,7 +19,7 @@ namespace boost{ * * @warning due to the recursive nature of real_operation, destruction may cause stack overflow */ - enum class OPERATION{ADDITION, SUBTRACTION, MULTIPLICATION}; + enum class OPERATION{ADDITION, SUBTRACTION, MULTIPLICATION, DIVISION}; class real_operation{ private: @@ -37,7 +37,7 @@ namespace boost{ */ real_operation(std::shared_ptr &lhs, std::shared_ptr &rhs, OPERATION op) : _lhs(lhs), _rhs(rhs), _operation(op) {}; - OPERATION get_operation() { + OPERATION get_operation() const { return _operation; } @@ -46,6 +46,15 @@ namespace boost{ /// fwd decl'd, defined in real_data const_precision_iterator& get_rhs_itr(); + + + std::shared_ptr rhs() const { + return _rhs; + } + + std::shared_ptr lhs() const { + return _lhs; + } }; } } diff --git a/test/include/test_helpers.hpp b/test/include/test_helpers.hpp index 4a5cbd0..e609feb 100644 --- a/test/include/test_helpers.hpp +++ b/test/include/test_helpers.hpp @@ -4,7 +4,7 @@ #include #include -std::optional boost::real::const_precision_iterator::maximum_precision = 10; +std::optional boost::real::const_precision_iterator::maximum_precision = 10; namespace Catch { template<> diff --git a/test/real_iterator_division_test.cpp b/test/real_iterator_division_test.cpp new file mode 100644 index 0000000..5c2db3b --- /dev/null +++ b/test/real_iterator_division_test.cpp @@ -0,0 +1,140 @@ +#include + +#include +#include + +TEST_CASE("Operator / boost::real::const_precision_iterator") { // assumes max precision is 10. + + SECTION("144/12") { + boost::real::real a("144"); // [100, 200] + boost::real::real b("12"); // [10 , 20] + boost::real::real result = a/b; // [100, 200] / [10, 20] = [100/20, 200/10] = [5, 20] + + auto result_it = result.get_real_itr().cbegin(); + + boost::real::interval expected_interval({}); + expected_interval.lower_bound.positive = true; + expected_interval.upper_bound.positive = true; + expected_interval.lower_bound.exponent = 1; + expected_interval.upper_bound.exponent = 2; + expected_interval.lower_bound.digits = {5}; + expected_interval.upper_bound.digits = {2}; + CHECK(expected_interval == result_it.get_interval()); + + ++result_it; // [140, 150] / [12, 12] = [140/12, 150/12] = [11.666...6, 12.5] + expected_interval.lower_bound.exponent = 2; + expected_interval.upper_bound.exponent = 2; + expected_interval.lower_bound.digits = {1,1,6,6,6,6,6,6,6,6}; + expected_interval.upper_bound.digits = {1,2,5}; + CHECK(expected_interval == result_it.get_interval()); + + ++result_it; // [144, 144] / [12, 12] = 12 + expected_interval.lower_bound.exponent = 2; + expected_interval.upper_bound.exponent = 2; + expected_interval.lower_bound.digits = {1,2}; + expected_interval.upper_bound.digits = {1,2}; + CHECK(expected_interval == result_it.get_interval()); + } + + SECTION("12/144") { // currently doesn't work for intermediate steps. We can fix it later. it gets close + // to expected, but not exact + boost::real::real a("144"); // [100, 200] + boost::real::real b("12"); // [10 , 20] + boost::real::real result = b/a; + // 12 / 144 + // [10, 20] / [100, 200] = [10/200, 20/100] = [.05, .2] + auto result_it = result.get_real_itr().cbegin(); + + boost::real::interval expected_interval({}); + expected_interval.lower_bound.positive = true; + expected_interval.upper_bound.positive = true; + expected_interval.lower_bound.exponent = -1; + expected_interval.upper_bound.exponent = 0; + expected_interval.lower_bound.digits = {5}; + expected_interval.upper_bound.digits = {2}; + CHECK(expected_interval == result_it.get_interval()); + + result_it = result.get_real_itr().cend(); + // [12, 12] / [144, 144] = [.08333...3, .0833...4] + expected_interval.lower_bound.positive = true; + expected_interval.upper_bound.positive = true; + expected_interval.lower_bound.exponent = -1; + expected_interval.upper_bound.exponent = -1; + expected_interval.lower_bound.digits = {8,3,3,3,3,3,3,3,3,3}; + expected_interval.upper_bound.digits = {8,3,3,3,3,3,3,3,3,4}; + CHECK(expected_interval == result_it.get_interval()); + } + + SECTION("1/3") { + boost::real::real a("1"); // [1, 1] + boost::real::real b("3"); // [3, 3] + boost::real::real result = a/b; // [1,1] / [3,3] + + auto result_it = result.get_real_itr().cbegin(); // same as .cend(), in this case + + boost::real::interval expected_interval({}); + expected_interval.lower_bound.positive = true; + expected_interval.upper_bound.positive = true; + expected_interval.lower_bound.exponent = 0; + expected_interval.upper_bound.exponent = 0; + expected_interval.lower_bound.digits = {3,3,3,3,3,3,3,3,3,3}; + expected_interval.upper_bound.digits = {3,3,3,3,3,3,3,3,3,4}; + CHECK(expected_interval == result_it.get_interval()); + } + + SECTION("-1/3") { + boost::real::real a("-1"); // [1, 1] + boost::real::real b("3"); // [3, 3] + boost::real::real result = a/b; // [1,1] / [3,3] + + auto result_it = result.get_real_itr().cbegin(); // same as .cend(), in this case + + boost::real::interval expected_interval({}); + expected_interval.lower_bound.positive = false; + expected_interval.upper_bound.positive = false; + expected_interval.lower_bound.exponent = 0; + expected_interval.upper_bound.exponent = 0; + expected_interval.lower_bound.digits = {3,3,3,3,3,3,3,3,3,4}; + expected_interval.upper_bound.digits = {3,3,3,3,3,3,3,3,3,3}; + CHECK(expected_interval == result_it.get_interval()); + } + + SECTION("-15/10") { + boost::real::real a("-15"); // [-1, -2] + boost::real::real b("10"); // [1, 2] + boost::real::real result = a/b; // [-1,-2] / [1,2] -> [-2, -1] + + auto result_it = result.get_real_itr().cbegin(); + + boost::real::interval expected_interval({}); + expected_interval.lower_bound.positive = false; + expected_interval.upper_bound.positive = false; + expected_interval.lower_bound.exponent = 1; + expected_interval.upper_bound.exponent = 1; + expected_interval.lower_bound.digits = {2}; + expected_interval.upper_bound.digits = {1}; + CHECK(expected_interval == result_it.get_interval()); + + ++result_it; + expected_interval.lower_bound.digits = {1,5}; + expected_interval.upper_bound.digits = {1,5}; + CHECK(expected_interval == result_it.get_interval()); + } + + SECTION("1/2") { + boost::real::real a("1"); // [1, 1] + boost::real::real b("2"); // [3, 3] + boost::real::real result = a/b; // [1,1] / [3,3] + + auto result_it = result.get_real_itr().cbegin(); // same as .cend(), in this case + + boost::real::interval expected_interval({}); + expected_interval.lower_bound.positive = true; + expected_interval.upper_bound.positive = true; + expected_interval.lower_bound.exponent = 0; + expected_interval.upper_bound.exponent = 0; + expected_interval.lower_bound.digits = {5}; + expected_interval.upper_bound.digits = {5}; + CHECK(expected_interval == result_it.get_interval()); + } +} \ No newline at end of file