Boost accurate rationals demo

boost_accurate_rationals_demo.cpp

#include <iostream>
#include <boost/multiprecision/cpp_int.hpp>

typedef boost::multiprecision::cpp_int bigint;
typedef boost::rational<bigint> rational;

/*
  Online running: https://ideone.com/zTH8pK

  This is example of absolutely accurate rational arithmetics with boost.
  Including double/float to rational absolutely accurate convertsion,
  because this conversion don't implemented in boost:
  http://www.boost.org/doc/libs/1_60_0/libs/rational/rational.html#Conversion%20from%20floating%20point
*/

template<typename T>
rational to_rational(T f) {
	int exp;
	f = std::frexp(f, &exp);  // now f in range: (-1;-0.5] or [0.5; 1)
	
	bigint n = 0;
	bigint d = 1;

	if (f < 0) {
		f *= -1;
		d = -1;
	}
	// now f in range: [0.5; 1)

	while (f != 0.0) {
		f *= 2.0;
		n *= 2;
		d *= 2;
		if (f >= 1.0) {
			n += 1;
			f -= 1.0;
		}
	}

	while (exp > 0) {
		n *= 2;
		exp--;
	}
	while (exp < 0) {
		d *= 2;
		exp++;
	}
	return rational(n, d);
}

int main() {
	double big = 1 << 30;
	big = big * big * big * big;  // 2^120
	rational xs[] = {to_rational(123.456), to_rational(-0.000239), to_rational(big)};
	for (int j = 0; j < 3; j++) {
		std::cout << "xs[" << j << "] = " << boost::rational_cast<double>(xs[j]) << " = " << xs[j] << std::endl;
	}
	
	rational mul = to_rational(239.0);
	rational add = to_rational(2012.0);
	
	for (int i = 0; i < 100; i++) {
		for (int j = 0; j < 3; j++) {
			xs[j] = mul * xs[j] + add;
		}
	}
	for (int i = 0; i < 100; i++) {
		for (int j = 0; j < 3; j++) {
			xs[j] = (xs[j] - add) / mul;
		}
	}
	
	std::cout << "No precision losses:" << std::endl;
	for (int j = 0; j < 3; j++) {
		std::cout << "xs[" << j << "] = " << boost::rational_cast<double>(xs[j]) << " = " << xs[j] << std::endl;
	}
}