Синус без стандартной библиотеки

В этой небольшой статье, которая будет представлена скорее как вычислительный эксперимент, я покажу каким образом можно вычислить синус некоего аргумента с плавающей точкой одинарной точности. Не ждите от этой статьи большего, просто это один из опытов, который полезно сделать и оставить в качестве описания, как полезного трюка, нежели реального кода.

Но…

Для начала опытов, позаимствуем шаблон для вычисления факториала из статьи «Compile time функциональное программирование в D» c habr.ru за авторством пользователя NCrashed:

// Параметризиуем только целочисленными значениями без знака
template factorial(uint n)
{
    // Шаблон-помошник 
    private template inner(ulong acc, uint n)
    {
        // В шаблонах мы можем использовать только static версию if
        static if(n == 0)
            enum inner = acc; // хвост
        else
            enum inner = inner!(acc * n, n-1 ); // уходим в глубину
    }
    // Вызываем внутренний шаблон с аккумулятором равным 1
    enum factorial = inner!(1, n);
}

Как видите, это стандартное вычисление факториала с помощью хвостовой рекурсии, реализованное в шаблонах, поскольку циклы напрямую нам недоступны. В примере используется одноименный шаблон-помощник, который помогает реализовать рекурсивный вызов шаблона посредством замены результата его уже «размотанным» телом, вычисляемым на каждом этапе подстановки значения аккумулятора и на единицу уменьшенной версии начального аргумента.

Шаблон факториала упомянут не просто так и будет использоваться нами для того, чтобы реализовать вычисление синуса с использованием замечательного разложения тригонометрических функций (в нашем случае, эта функция синус) в ряд Макклорена.

При этом базовый паттерн построения шаблона практически совпадает с уже приведенным кодом за авторством NCrashed:

/// синус на шаблонах
template taylorSinus(float x, uint n)
{
    private template inner(float acc, uint n)
    {
        static if (n == 0)
            enum inner = acc;
        else
            enum inner = inner!(
                acc + (((-1.0f) ^^ n) * ((x ^^ (2 * n + 1)) / factorial!(2 * n + 1))),
                n - 1
            );
    }
    enum taylorSinus = inner!(x, n);
}

/// упрощенный вариант
template sinus(float x)
{
    alias sinus = taylorSinus!(x, 10);
}

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

Для испытания кода мы вычислим синус для двух значений (значения брались из головы, так что не обессудьте) стандартным способом и с помощью шаблонов, выведя результат вычисления шаблонного кода с помощью директивы pragma(msg):

void main()
{
    import std.math : sin;
    import std.stdio : writeln;

   pragma(msg, "Шаблонный синус от 0.5: ", sinus!0.5f);
   pragma(msg, "Шаблонный синус от 0.205: ", sinus!0.205f);
   writeln("Синус от 0.5: ", sin(0.5f));
   writeln("Синус от 0.205: ", sin(0.205f));
}

Результат:

Шаблонный синус от 0.5: 0.479426F
Шаблонный синус от 0.205: 0.203567F
Шаблонный синус от 0.5: 0.479426F
Шаблонный синус от 0.205: 0.203567F
Синус от 0.5: 0.479426
Синус от 0.205: 0.203567

Как видите, реализовать вычисления на стадии компиляции не так сложно, а иногда и полезно, а я рекомендую вам познакомится с упоминавшейся статьей для погружения в мир CTFE (compile time function evaluation, вычисления на стадии компиляции) и шаблонов.

Добавить комментарий