Calculations with a floating point at a compilation stage

As you know, in C ++ you cannot perform complex floating-point calculations at the compilation stage. I decided to try to get rid of this annoying flaw. The goal to which we will go, for example, the calculation of the root:

typedef RATIONAL(2,0) x; typedef sqrt<x>::type result; 

The root of the number is calculated at the compilation stage. The representation of the number is stored as a ratio of two integers, so in order to get the value, you need to call through the get () method;

 std::cout << result::get() << std::endl; 

1.41421356

For storing floating point numbers, we will use the ratio
two integers - fraction.

 template <int64_t A, int64_t B> struct rational_t { const static int64_t a = A, b = B; static double get() { return (double)a/b; } }; 

For example, to declare the number 2, you need to declare the type rational_t ​​<2.1> (the first two)

 rational_t<3,2> -> 3/2 rational_t<56,10> -> 56/10 = 5.6 rational_t<3,100> -> 3/100 = 0.03 

For convenience, we will use the clever RATIONAL macro:

 #define RATIONAL(A1, A2) rational_t<(int)(A1##A2), pow<10, sizeof(#A2)-1>::value> 

Consider it on an example: RATIONAL (12,34) - announces the fraction 12.34 (12, 34 hundredths)

 1) A1 ## A2 sticks together two arguments in 1234
 2) sizeof (# A2) -1 = sizeof ("34") - 1 = 3 - 1 = 2
 3) pow <10, 2> :: value = 10 to the power of 2 = 100

Thus, RATIONAL (12,34) declares the type rational_t ​​<1234, 100> (that is, 1234/100 = 12.34)
The pow function is written like this:

 template <int V, unsigned D> struct pow { const static int value = V * pow<V, D - 1>::value; }; template <int V> struct pow<V, 0> { const static int value = 1; }; 

We describe the arithmetic operations for the rational_t ​​template.
The addition of fractions: a1 / b1 + a2 / b2 = (a1 * b2 + a2 * b1) / (b1 * b2), therefore:

 template <class R1, class R2> struct plus { typedef rational_t<R1::a * R2::b + R2::a * R1::b, R1::b * R2::b> type; }; 

With constant addition, the numerator and denominator will constantly grow in absolute value. To avoid overflow at the next addition, it is necessary to reduce the fraction. Rewrite the addition function:

 template <class R1, class R2> struct plus { typedef rational_t<R1::a * R2::b + R2::a * R1::b, R1::b * R2::b> type1; typedef typename reduce<type1>::type type; }; 

The reduce meta-function takes an unabridged fraction and returns an abbreviated one. We present the remaining operations:

 template <class R1, class R2> struct minus { typedef rational_t<R1::a * R2::b - R2::a * R1::b, R1::b * R2::b> type1; typedef typename reduce<type1>::type type; }; template <class R1, class R2> struct mult { typedef rational_t<R1::a * R2::a, R1::b * R2::b> type1; typedef typename reduce<type1>::type type; }; template <class R1, class R2> struct divide { typedef rational_t<R1::a * R2::b, R1::b * R2::a> type1; typedef typename reduce<type1>::type type; }; 

Comparison operation:

 template <class R1, class R2> struct less { static const bool value = (R1::a * R2::b - R2::a * R1::b) < 0; }; 

We define a metafunction defining the need for reduction. The maximum allowable number for storage in 64-bit format will be (2 ^ 64) ^ 0.5 / 2 = 1ll << 31, therefore:

 template <class R> struct require_reduce { const static int64_t max = (1ll<<31); const static bool value = (R::a >= max) || (R::b >= max); }; 

value means the need to reduce the fraction, otherwise there may be an overflow during the operation.
Fraction must first be reduced exactly - by divisibility by GCD (reduce_accurate),
if the reduction was not possible, then the reduction is inaccurate, dividing the numerator and
denominator is a complete 2 (reduce_inaccurate).
Let's declare the neatly reduced metafunction:

 template <bool, class R> struct reduce_accurate; 

Inaccurate reduction, after it an attempt is made to reduce again
neatly if, of course, this is required:

 template <bool, class R> struct reduce_inaccurate { typedef rational_t<(R::a >> 1), (R::b >> 1)> type_; typedef typename reduce_accurate<require_reduce<type_>::value, type_>::type type; }; 

If not required, an inaccurate abbreviation returns the same value:

 template <class R> struct reduce_inaccurate<false, R> { typedef R type; }; 

For a careful abbreviation in the comments, they suggested using NOD, it turned out to be faster with it, although it was expected that it would not.
Calculation NOD (thanks gribozavr ):

 template <int64_t m, int64_t n> struct gcd; template <int64_t a> struct gcd<a, 0> { static const int64_t value = a; }; template <int64_t a, int64_t b> struct gcd { static const int64_t value = gcd<b, a % b>::value; }; 

The reduction function divides the numerator and denominator by the GCD and, if necessary, performs an inaccurate reduction.

 template <bool, class R> struct reduce_accurate { template <bool, class R> struct reduce_accurate { const static int64_t new_a = R::a / gcd<R::a, R::b>::value; const static int64_t new_b = R::b / gcd<R::a, R::b>::value; typedef rational_t<new_a, new_b> new_type; typedef typename reduce_inaccurate<require_reduce<new_type>::value, new_type>::type type; }; }; 

If an exact reduction is not enough, then perform an inaccurate reduction:

 template <class R> struct reduce_accurate<false, R> { typedef typename reduce_inaccurate<require_reduce<R>::value, R>::type type; }; 

We turn to the most interesting, write the algorithm for calculating the root according to Newton's method.
This implementation does not claim accuracy and is given as an example.

  : sqrt_eval(p, res, x) { t1 = x/res t2 = res+t1 tmp = t2/2 if (p-1 == 0) return tmp; return sqrt_eval(p-1, tmp, x) } 

Using rational_t:

 template <int64_t p, class res, class x> struct sqrt_eval { typedef typename divide<x, res>::type t1; typedef typename plus<res, t1>::type t2; typedef typename divide<t2, rational_t<2,1> >::type tmp; typedef typename sqrt_eval<p-1, tmp, x>::type type; }; template <class res, class x> struct sqrt_eval<0, res, x> { typedef res type; }; 

')

 sqrt(x) { res = (x + 1)/2 return sqrt_eval(15, res, x) } 

15 - number of steps in the Newton algorithm.

 template <class x> struct sqrt { typedef typename divide< typename plus<x, rational_t<1,1> >::type, rational_t<2,1> >::type res; typedef typename sqrt_eval<15, res, x>::type type; }; template <int64_t a> struct sqrt< rational_t<0, a> > { static const int64_t value = 0; }; 

Example:

 #include <iostream> #include "rational_algo.hpp" int main() { std::cout.precision(15); const double s = rational::sqrt<RATIONAL(2,0)>::type::get(); std::cout << s << std::endl; std::cout << 2-s*s << std::endl; return 0; } 

An example of a solid file: pastebin.com/ea7S2KTd
Project: github.com/korkov/rational

_{UPD: this is an updated version, thanks everyone for the tips and fixes.}