Вычисляем значение числа e на этапе компиляции

Проглядывая книжку «Эффективное использование C++» ^[1], Скотта Мейерса, которая ( и я никого не удивлю ) достойна всяческих похвал, меня очень тронуло, то с какой возбуждённостю, вдохновлённостю, трепетом ( может мне показалось? ) автор говорит о шаблонах и их возможностях. Приведу маленький кусочек:

Метапрограммирование шаблонов ( template metaprogramming — TMP ) — это процесс написания основанных на шаблонах программ на C++, исполняемых во время компиляции. На минуту задумайтесь об этом: шаблонная метапрограмма — это программа, написанная на C++, которая исполняется внутри компилятора C++…
Было доказано, что технология TMP предоставляет собой полную машину Тьюринга ^[2], то есть обладает достаточной мощь для любых вычислений...

Да уж… сердце заколотало, в очередной раз удивился — только подумать — полная машина Тьюринга со всемы вытекающими последствиями… Как по мне, это просто невероятно и удивительно… хотя, кто его знает…

Предлагаю посмотреть на совсем уж маленький кусочек мира больших возможностей и невероятных приключений — попробуем вычислить на этапе компиляци значение, небезызвестного, числа e ^[3].

Как её вичислять (с некоторой погрешностью) — подскажет ряд Тейлора ^[4], а точнее ряд Маклорена:

Т.е. нам нужно будет уметь считать факториал числа, подносить в степень, суммировать и работать с дробными числами… и всё это с помощью шаблонов C++.
Для начало хотелось бы разобраться с дробными числами — нужно как-то сохранять числитель и знаменатель, и иметь также доступ к ним ( N — Numerator, D — Denominator ):

template<int n, int d>
struct Fractional
{
	enum { N = n, D = d };
};

Всё просто, но как насчёт нулевого знаменателя? Попробуем это:

template<int n, int d>
struct Fractional
{
private:
	enum { NonZeroDenominator = n / d };

public:
	enum { N = n, D = d };
};

Используем:

typedef Fractional<9, 0> number;
// ...
int temp = number::D;

В случае с msvc10 мы получим что-то вроде error C2057: expected constant expression — невнятно, но если пойти к месту ошибки — то как раз увидим переменную NonZeroDenominator — уже хоть что-то…

Итак, сохранять 2 числа умеем, а как же насчёт сокращения дробей? Тут надо уже научиться находить gcd (Наиболее общий делитель) ^[5] двух чисел — нам подходит рекурсивный алгоритм:

int gcd(int a, int b)
{
    if(b == 0) return a;
    return gcd(b, a % b);
}

который превращается с помощь шаблонов в:

template<int n1, int n2>
struct GCD
{
	enum { value = GCD<n2, n1 % n2>::value };
};

template<int n1>
struct GCD<n1, 0>
{
	enum { value = n1 };
};

Всё просто, не так ли? — пишем наиболее общую реализацию шаблона и делаем частную специализацию для частных случаев (если второе чило ноль — результат — это первое число).
С помощью всего выше написанного, делаем окончательную версию дробного числа:

template<int n, int d>
struct Fractional
{
private:
	enum { NonZeroDenominator = n / d };
	enum { gcd = GCD<n, d>::value };

public:
	enum { N = n / gcd, D = d / gcd };
};

С помощью известных формул — делим, множим, отнимаем, додаём наши числа:

//
// Divide
//
template<typename n, typename d>
struct Divide
{
};

template<int n1, int d1, int n2, int d2>
struct Divide<Fractional<n1, d1>, Fractional<n2, d2> >
{
private:
	typedef Fractional<n1, d1> n;
	typedef Fractional<n2, d2> d;

public:
	typedef Fractional<n::N * d::D, n::D * d::N> value;
};

//
// Multiple
//
template<typename n, typename d>
struct Multiple
{
};

template<int n1, int d1, int n2, int d2>
struct Multiple<Fractional<n1, d1>, Fractional<n2, d2> >
{
private:
	typedef Fractional<n1, d1> n;
	typedef Fractional<n2, d2> d;

public:
	typedef Fractional<n::N * d::N, n::D * d::D> value;
};

//
// Substract
//
template<typename n, typename d>
struct Substract
{
};

template<int n1, int d1, int n2, int d2>
struct Substract<Fractional<n1, d1>, Fractional<n2, d2> >
{
private:
	typedef Fractional<n1, d1> n;
	typedef Fractional<n2, d2> d;

public:
	typedef Fractional<n::N * d::D - d::N * n::D, n::D * d::D> value;
};

//
// Add
//
template<typename n, typename d>
struct Add
{
};

template<int n1, int d1, int n2, int d2>
struct Add<Fractional<n1, d1>, Fractional<n2, d2> >
{
private:
	typedef Fractional<n1, d1> n;
	typedef Fractional<n2, d2> d;

public:
	typedef Fractional<n::N * d::D + d::N * n::D, n::D * d::D> value;
};

Снова же — пишем пустой набросок нашей "фунции", например Divide — отмечая, что она (функция) принимает 2 аргумента. И дальше с помощью частичной специализации шаблона уточняем, что хотим видеть не что нибудь, а именно идентификатор шаблона нужного нам, т.е. Divide<n1, n2>, к примеру. Использование:

	typedef Fractional<4, 20> n1;
	typedef Fractional<5, 32> n2;

	typedef Add<n1, n2>::value summ;
	printf("%i/%in", summ::N, summ::D);
    // 57/160

Также нам нужно возведение в степень и факториал, определение которых говорит само о себе:

//
// Factorial
//
template<int N>
struct Factorial
{
	enum { value = N * Factorial<N - 1>::value };
};

template<>
struct Factorial<0>
{
	enum { value = 1 };
};

//
// Power
//
template<int x, int n>
struct Pow
{
	enum { value = x * Pow<x, n - 1>::value };
};

template<int x>
struct Pow<x, 0>
{
	enum { value = 1 };
};

Итак, теперь у нас есть весь набор всего необходимого, чтобы реализовать формулу выше — понятно суммировать мы будем не до бесконечности, а сколько сможем, т.е. например, выражение Exp<4, 8>::value будет давать дробное число, которое численно равно экспоненте в 4й степени и суммирование произведено всего лишь до 8 (бесконечность рядом) члена включительно.

Проблема возникает в том, а как нам суммировать дробные числа, которые даже не являются чисельными значениями — это всего лишь типы! Да, они содержат в себе чисельные данные, но к ним еще нужно добраться в ходе подсчёта суммы ряда… Но решение есть и оно состоит в том, что мы можем достать из производного класса данные (и typedef-ы — самое важное) базового класса. Именно! — чтобы подсчитать сумму ряда, нам нужно будет наследваться и наследоваться, и наследоваться… в идеале до бесконечности.
Кусок кода:

//
// Exponent
//
template<int x, int n>
struct Exp :
	public Exp<x, n - 1>
{
private:
	typedef typename Exp<x, n - 1>::value previous;
protected:
	typedef Fractional<Pow<x, n>::value, Factorial<n>::value> current;
public:
	typedef typename Add<current, previous>::value value;
};

template<int x>
struct Exp<x, 0>
{
public:
	typedef Fractional<Pow<x, 0>::value, Factorial<0>::value> current;
public:
	typedef current value;
};

current содержит в себе значение одного члена ряда — т.е. один класс в этой целой иерархии классов, грубо говоря, предназначен для хранение значения одного члена. А с помощью того, что он может взять данные базового класса — т.е. значение предыдущего члена ряда — то всё это даёт нам то, что кроме значения одного отдельного элемента ряда, мы в текущем классе имеем сумму ряда до этого элемента (класса) включительно.

И что в итоге, а в итоге мы с гордостью можем написать следущее:

int main()
{
	// дробь
	typedef Exp<1, 8>::value result;
	printf("%i/%in", result::N, result::D);

	// десятичное представление
	printf("%fn", result::N / static_cast<float>(result::D));
}

на 32 битной машине больше 8ми членов ряда взять не получится — переполнение int.

Результат: 2.718279 (109601/40320).

Волшебство :)
Надеюсь, Вам было приятно. Спасибо за внимание.

Автор: Door

Источник ^[6]