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

    Как известно, в 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: это обновлённая версия, спасибо всем за советы и исправления.
    Ads
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More

    Comments 28

      +2
      Ох жесть… Я дочитал)
      Это конечно очень круто, что можно на этапе компиляции что-то сложное вычислять, но мне все-же проще посчитать на калькуляторе за 10 секунд, чем неделю писать и как-то отлаживать такую вот штуковину.
        +3
        Это конечно очень круто, что можно на этапе компиляции что-то сложное вычислять, но мне все-же проще посчитать на калькуляторе за 10 секунд
        Вы проглядели последний тэг :)
        как-то отлаживать такую вот штуковину.
        Боюсь это невозможно. Делается же все в compile-time.
          0
          Вот и я про то)
          +4
          вообще сложно что-то отлаживать на шаблонах в C++
          +2
          Интересное продолжение!

          > Для хранения чисел с плавающей точкой будем использовать соотношение двух целых чисел — дробь.
          Не-а. Плавающая точка это представление чисел в виде мантиссы и порядка. У вас получилась арифметика рациональных чисел, пусть и неточная.

          Единственное что могу добавить: для сокращения лучше бы искать НОД. Или слишком долго компилируется?
            0
            >Интересное продолжение!
            Спасибо :)

            >для сокращения лучше бы искать НОД. Или слишком долго компилируется?
            После неточного сокращения (reduce_inaccurate) делается попытка снова сократить точно, то есть снова бы пришлось искать НОД, поэтому я сразу отказался от этого. Простые числа до 100 считаются же всего 1 раз.
              0
              Простые числа — да, один раз. Но перебор и попытки разделить нацело выполняются каждый раз. Так как максимальное количество шагов в алгоритме Эвклида будет меньше сотни для 64-разрядных чисел и pi(100)=25, разница в производительности должна быть не очень большой, а выигрыш в точности получить можно.
                +1
                Я написал сокращение через НОД, если интересно вот, то что изменилось:
                pastebin.com/kRy21zPM

                Результаты времени компиляции без НОД для sqrt(1000):

                real 0m5.859s
                user 0m2.628s
                sys 0m0.148s
                31.622776825553
                -1.41587300855894e-05

                С НОД:
                real 0m13.373s
                user 0m6.540s
                sys 0m0.296s
                31.622776825553
                -1.41587300855894e-05

                То есть с НОД дольше (но для меньших чисел алгоритм с НОД выигрывает)
                точность -1.41587300855894e-05 не изменилась.
                  +1
                  А если реализовать через остаток от деления (алгоритм Эвклида, а не двоичный алгоритм Эвклида): pastebin.com/H9TvWdeW

                  На g++ 4.6.1 без НОД sqrt(10000): real 0m1.359s
                  С НОД: real 0m0.358s
                    0
                    здорово, обновил библиотеку
                      +4
                      While you are at it, NOD -> GCD (greatest common divisor).
                +1
                А что, НОД так долго ищется? Там же оценка сверху на количество шагов Log_2(n+m), если я правильно помню. Что для 64-битных чисел дает 65 шагов.
            +2
            Месье знает толк, ага. Идея Александреску живут и побеждает!
              +10
              Стоит заметить что пользователи GCC >=4.5 могут вместо всего вышеприведённого написать constexpr-функцию (часть стандарта C++0x) и в ней спокойно работать с желаемым типом данных, будь то структура или число с плавающей точкой.
                0
                Есть библиотека в бусте boost::rational. Зарелизили в 1.47 (последней). Делает то же.
                  +5
                  названия совпадают, но boost::rational это совершенно другое,
                  и появилась она не в 1.47, в 1.31 уже точно была
                    0
                    Спасибо за информацию, запомню.
                  0
                  just for fun вещь конечно. но как упоминалось выше мисье знает толк в извращениях
                    0
                    Увидел пару интересных приёмов кунг-фу (например определение числа цифр в числе на этапе компиляции = sizeof(#n)-1 ). В общем, спасибо!
                      +1
                      Кстати, очень аккуратный и изящный код. А извращенцем зовут ниасиляторы! :P
                        +1
                        основной алгоритм:
                        sqrt_eval(p, res, x) {
                        t1 = x/res
                        t2 = res+t1
                        tmp = t2/2
                        less_val = tmp < res? tmp: res
                        if (p-1 == 0)
                        return less_val;
                        return sqrt_eval(p-1, less_val, x)
                        }
                        С какой целью выбирается наименьшее значение?
                          0
                          согласен, он не нужен, спасибо
                          переписал с алгоритма для нахождения корня в целых числах,
                          а там, чтобы обеспечить sqrt(x)*sqrt(x) =< x
                            0
                            Прошу прощения, смотрите комментарий ниже.
                          0
                          Я поясню свой вопрос: выбирая наименьшее значение из текущего и следующего приближений, в случае если вы выбираете начальное приближение < sqrt(x), алгоритм вернёт его же после произвольного числа итераций.
                            0
                            Хорошая разминка для мозга, в принципе все элементарно. Вот когда я решил написать свою библиотеку compile-time для работы со списками, векторами данных, для работы со строками во время компиляции, чуть мозг себе не сломал.
                              0
                              boost::ratio же, не?
                              А на практике удобнее boost::rational — он не на этапе компиляции, но зато имеет меньшее число ограничений.
                                0
                                boost::ratio это то, что нужно, спасибо.
                                К своему оправданию могу сказать, что boost::ratio вышла в последней версии 1.47, я использовал 1.46 :(
                                Комментатор выше, который говорил про boost::rational, скорее всего имел именно boost::ratio.

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