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

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

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


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

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

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


    Т.е. нам нужно будет уметь считать факториал числа, возводить в степень, суммировать и работать с дробными числами… и всё это с помощью шаблонов 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 (Наиболее общий делитель) двух чисел — нам подходит рекурсивный алгоритм:
    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/%i\n", 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/%i\n", result::N, result::D);
    
    	// десятичное представление
    	printf("%f\n", result::N / static_cast<float>(result::D));
    }
    

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

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

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

    PS: отдельно извиняюсь за орфографические, синтаксические ошибки — был в сонном несостоянии, недоумении и восторге.

    Similar posts

    Ads
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More

    Comments 24

      +6
      1) «заколотало», «чисельные», «додаём» — о значении некоторых перлов вашего новояза остаётся только догадываться, и то безуспешно. Установите в браузер проверку орфографии, пожалуйста. Мне больно за читателей и других форумов, в которые вы пишете, причём чаще, чем на Хабре.

      2) способы проигнорировать функцию main() при выполнении программы существуют, да.

      3) Вводная по метапрограммированию — тоже любопытна.
        +1
        извините, исправил.
          +1
          1) «чисельные», «додаём» — это все кальки с украинской математической терминологии, автор украиноязычен и слушал лекции по матану только на украинском языке. Никакого новояза.
          +22
          Читать Майерса и Александреску это хорошо. Но я бы за такие развлечения в коде проекта без явно указанной необходимости в таком коде, отрывал руки без права пришития обратно.

          Шаблоны очень сложны в отладке и сообщения об ошибках выдаваемые gcc на шаблонах не страдают информативностью. Подобный код хорош, когда он написан раз и его трогать больше не надо. Или когда действительно встаёт задача какие-то вещи для оптимизации вынести в compile-time. В остальном — такой код — лютое нечитабельное и неподдерживаемое говно.

          Я не говорю про то, что шаблоны это плохо, но вот к метапрограммированию на них надо относиться очень аккуратно.
            +4
            Как хорошо, что благодаря constexpr, все эти танцы с шаблонами вскоре уйдут в прошлое.
              0
              Поддерживаю. Статья хороший пример мощи шаблонов и отличный стимул перейти на стандарт 11 года.
              +20
              Ну что вы, теперь так e вычисляить уже не модно!
              Вот он, новый способ:

              #include <iostream>
              
              constexpr int factorial(int n)
              {
                  return n > 1? factorial(n-1) * n : 1;
              }
              
              constexpr double cpow(double v, int n)
              {
                  return n ==0? 1: v * cpow(v, n-1);
              }
              
              constexpr double factor(double v, int n)
              {
                  return cpow(v, n) / factorial(n);
              }
              
              constexpr double sum(double v, int n, int max)
              {
                  return factor(v, n) + (max == n? 0: sum(v, n + 1, max));
              }
              
              constexpr double cexp(double v, int accuracy)
              {
                  return sum(v, 0, accuracy);
              }
              
              int main(int argc, const char** argv) 
              {
                  constexpr double e = cexp(1, 8);
                  std::cout.precision(15);
                  std::cout << e << "\n";
              }
              
              
                0
                И функцию Аккермана? ;-)
              • UFO just landed and posted this here
                  0
                  От непосредственного вычисления значения факториала можно избавиться путем выражения последующего члена ряда Тейлора через предыдущий. Насколько я помню, могу ошибаться, fpu так и вычисляет тригонометрию.
                    0
                    Насчет тригонометрии — вроде давно таблицы Брадиса во флэше записаны, вычислять только поправку нужно.
                    • UFO just landed and posted this here
                        +1
                        У вас неправильная сумма ряда. Вот результат на Maple:

                        2.7182818284590452353602874713526624977572470937000 — точное значение:
                        2.7182818284590452353602874713526624976385970411899 — Ряд Тейлора до 1/32!

                        Цепную дробь пока не посчитал.
                          +1
                          2.7182818284590452353602874713558030927124686274627 — цепная дробь (33 члена разложения)
                          • UFO just landed and posted this here
                              0
                              Зато если считать в рациональных числах (не сокращая их) и ограничивать числитель и знаменатель какой-нибудь константой, то цепные дроби дадут почти вдвое лучший результат по числу знаков. Так, если взять ограничение 2^32, то цепную дробь можно считать до [...,14,1,1], и получается 848456353/312129649 (17 верных знаков). Ряд Тейлора считается до 1/12!, и получается 1302061345/479001600 (9 верных знаков). Удивительно, что сумма обратных факториалов иногда сокращается.
                    +10
                    — сейчас я покажу вам особую, шаблонную магию!
                    — в рот мне ноги! нет-нет-нет, александреску, немедленно верни все как было!!!
                      +5
                      Не обращайте внимания на граммар-нази и лиц, лишенных чувства юмора. Статья здоровская! Только напишите в начале крупными буквами «Не пытайтесь повторить это в реальных проектах»:)
                        –1
                        Почему статья не в хабе «Ненормальное программирование». Требую перенести!
                        Шаблоны в С++ это как подъязык в языке. Т.е. сейчас вы решили задачу поиска экспоненты на специфическом языке, а не на С++. Статья хорошая в плане того, что бы взять на заметку и когда приспичит писать шаблоны, посмотреть и вспомнить как можно извратиться :)
                          0
                          Думаю, именно из-за такого безобразия в метапрограммировании на шаблонах Александреску и сбежал в D.
                          Compile-time аналог на D (используя цепную дробь):
                          module ecompute;
                          
                          import std.stdio;
                          
                          // Вариант 1, используя CTFE
                          static double computeE1(int n)
                          {
                          	double innerCompute(int i, int m)
                          	{
                          		if(m <= 0)
                          		{
                          			return 0.0;
                          		}
                          		return 1.0/(
                          			    1.0 + 1.0/(
                          			      i + 1.0/(
                          			    1.0 + innerCompute(i + 2, m - 1))));
                          	}
                          	return 2.0 + innerCompute(2, n);
                          }
                          
                          // Вариант 2, олдскул
                          template computeE2(int n)
                          {
                          	private template innerCompute(int i, int m)
                          	{
                          		static if(m <= 0)
                          		{
                          			enum innerCompute = 0.0;
                          		} else
                          		{
                          			enum innerCompute = 1.0/(
                          					1.0 + 1.0/(
                          					  i + 1.0/(
                          					1.0 + innerCompute!(i + 2, m - 1))));
                          		}
                          	}
                          	enum computeE2 = 2.0 + innerCompute!(2, n);
                          }
                          
                          void main()
                          {
                          	// Вариант 1
                          	static immutable e1 = computeE1(20);
                          	pragma(msg, e1);
                          	printf("Case 1, e = %1.20f \n",e1);
                          
                          	// Вариант 2
                          	static immutable e2 = computeE2!(20);
                          	pragma(msg, e2);
                          	printf("Case 2, e = %1.20f \n",e2);
                          }
                          

                          D упрощает compile-time программирование, а люди идут в своих извращениях дальше: compile-time raytracer.
                            +2
                            Я тоже около 15 лет думал, что пишется «substract». Так и писал. А оказалось — «subtract». Видимо, это какое-то популярное заблуждение.
                            0
                            Возможно стоит добавить следующие ссылки:
                            Template metaprogramming на википедии
                            Boost.MPL — фреймворк для создания программ, которые производят вычисления во время компиляции с помощью шаблонов.
                              0
                              Ох, я тоже что-то подобное тут писал =)
                              habrahabr.ru/post/124963/

                              Only users with full accounts can post comments. Log in, please.