Вычисления с плавающей точкой на этапе компиляции

Как известно, в C++ нельзя производить сложные вычисления с плавающей точкой на стадии компиляции. Я решил попробовать избавиться от этого досадного недостатка. Цель, к которой мы будем идти, на примере вычисления корня:

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

Корень числа вычислится на этапе компиляции. Представление числа хранится как отношение двух целых чисел, поэтому чтобы получить значение, нужно обращаться через метод get();

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

1.41421356

Для хранения чисел с плавающей точкой будем использовать соотношение
двух целых чисел — дробь.

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; }
};

Например, чтобы объявить число 2, нужно объявить тип rational_t<2,1> (две первых),

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

Для удобства будем использовать хитрый макрос RATIONAL:

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

Рассмотрим его на примере: RATIONAL(12,34) — объявляет дробь 12.34 (12 целых 34 сотых)

1) A1##A2 склеивает два аргумента в 1234
2) sizeof(#A2)-1 = sizeof("34")-1 = 3 - 1 = 2
3) pow<10, 2>::value = 10 в степени 2 = 100

Таким образом, RATIONAL(12,34) объявляет тип rational_t<1234, 100> (то есть 1234/100 = 12.34)
Функция pow пишется так:

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;
};

Опишем арифметические операции для шаблона rational_t.
Сложение дробей: a1/b1 + a2/b2 = (a1*b2 + a2*b1)/(b1*b2), следовательно:

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

При постоянном сложении числитель и знаменатель будут постоянно по модулю расти. Чтобы не произошло переполнение при очередном сложении, необходимо дробь сокращать. Перепишем функцию сложения:

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;
};

Метафункция reduce принимает несокращенную дробь и возвращает сокращенную. Приведём остальные операции:

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;
};

Операция сравнения:

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

Определим метафункцию, определяющую необходимость сокращения. Максимальное допустимое число для хранения в 64-битном формате будет (2^64)^0.5/2 = 1ll<<31, следовательно:

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

value означает необходимость сокращения дроби, иначе может произойти переполнение при выполнении операции.
Дробь нужно сначала сокращать точно — по делимости на НОД (reduce_accurate),
если сокращение не удалось, то сокращать неточно, деля числитель и
знаменатель нацело на 2 (reduce_inaccurate).
Объявим метафункцию аккуратного сокращения:

template <bool, class R>
struct reduce_accurate;

Неаккуратное сокращение, после него делается попытка снова сократить
аккуратно, если, конечно, это требуется:

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;
};

Если не требуется, то неаккуратное сокращение возвращает то же значение:

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

Для аккуратного сокращения в комментариях предложили использовать НОД, с ним оказалось быстрее, хотя ожидалось, что нет.
Вычисление НОД(спасибо 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;
};

Функция сокращения делит числитель и знаменатель на НОД и, если необходимо, выполняет неточное сокращение.

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;
  };
};

Если точного сокращения не достаточно, то выполнить неточное сокращение:

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

Перейдём к самому интересному, напишем алгоритм вычисления корня по методу Ньютона.
Эта реализация не претендует на точность и приведена в качестве примера.

основной алгоритм:
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)
}

С использованием 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 — кол-во шагов в алгоритме Ньютона.

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;
};

Пример:

#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;
}

Пример цельным файлом: pastebin.com/ea7S2KTd
Проектом: github.com/korkov/rational

_{UPD: это обновлённая версия, спасибо всем за советы и исправления.}