Функция Math.Sin (double) для GPU

в 15:08, , рубрики: C#, double, gpgpu, Math.Sin, unity3d, математика, ряды Тейлора

Предисловие

Мне понадобилось вычислять дугу с повышенной точностью на процессоре видеокарты в режиме реального времени.

Автор не ставил перед собой цель превзойти стандартную функцию System.Math.Sin() (C#) и ее не достиг.

Результат работы и мой выбор (для тех, кому не хочется читать):

Sin_3(rad)

using System;
class Math_d
{
    const double PI025      = Math.PI / 4;
    /// <summary> 2^17 = 131072 (1 мб), ошибка меньше 10000 (с конца), при 2^21 = 22097152 (16 мб) минимальная ошибка +-1 (с конца) (почти нет ошибки) </summary>
    const int length_mem    = 22097152;
    const double step_rad   = PI025 / length_mem;
    /// <summary> Массив предварительно рассчитанного sin, используемый для линейной интерполяции. </summary>
    static double[] mem_sin;
    /// <summary> Массив предварительно рассчитанного cos, используемый для линейной интерполяции. </summary>
    static double[] mem_cos;

    /// <summary> Инициирует данные, вроде массивов sin, необходимые в ходе расчетов. </summary>
    public static void Initialise()
    {
        Ini_Mem_Sin();
        Ini_Mem_Cos();
    }
    public static double Sin_3(double rad)
    {
        double rad_025;
        int i;

        i = (int)(rad / PI025);     //Находим индекс части
        rad_025 = rad - PI025 * i;  //Находим отступ от начала части (в радианах)
        i = i & 7;                  //Находим остаток от деления индекса на 8

        //Ищем индекс кусочка круга
        switch (i)
        {
            case 0:
                return Sin_Lerp(rad_025);
            case 1:
                return Cos_Lerp(PI025 - rad_025);
            case 2:
                return Cos_Lerp(rad_025);
            case 3:
                return Sin_Lerp(PI025 - rad_025);
            case 4:
                return -Sin_Lerp(rad_025);
            case 5:
                return -Cos_Lerp(PI025 - rad_025);
            case 6:
                return -Cos_Lerp(rad_025);
            case 7:
                return -Sin_Lerp(PI025 - rad_025);
        }

        return 0;
    }

    /// <summary> Подготавливает массив sin для линейной интерполяции </summary>
    static void Ini_Mem_Sin()
    {
        double rad;
        mem_sin = new double[length_mem];
        for (int i = 0; i < length_mem; i++)
        {
            rad = (i * PI025) / length_mem;
            mem_sin[i] = Math.Sin(rad);
        }
    }
    /// <summary> Подготавливает массив cos для линейной интерполяции </summary>
    static void Ini_Mem_Cos()
    {
        double rad;
        mem_cos = new double[length_mem];
        for (int i = 0; i < length_mem; i++)
        {
            rad = (i * PI025) / length_mem;
            mem_cos[i] = Math.Cos(rad);
        }
    }
    /// <summary> Использует линейную интерполяцию для поиска sin от 0 до pi/4. </summary>
    /// <param name="rad"> Угол в радианах от 0 до pi/4. </param>
    static double Sin_Lerp(double rad)
    {
        int i_0;
        int i_1;
        double i_0d;
        double percent;
        double a;
        double b;
        double s;

        percent = rad / PI025;
        i_0d    = percent * length_mem;
        i_0     = (int)i_0d;
        i_1     = i_0 + 1;
        a       = mem_sin[i_0];
        b       = mem_sin[i_1];
        s       = i_0d - i_0;

        return Math_d.Lerp(a, b, s);
    }
    /// <summary> Использует линейную интерполяцию для поиска cos от 0 до pi/4. </summary>
    /// <param name="rad"> Угол в радианах от 0 до pi/4. </param>
    static double Cos_Lerp(double rad)
    {
        int i_0;
        int i_1;
        double i_0d;
        double percent;
        double a;
        double b;
        double s;

        percent = rad / PI025;
        i_0d = percent * length_mem;
        i_0 = (int)i_0d;
        i_1 = i_0 + 1;
        a = mem_cos[i_0];
        b = mem_cos[i_1];
        s = i_0d - i_0;

        return Math_d.Lerp(a, b, s);
    }
    /// <summary> Производит линейную интерполяцию между двумя значениями. (return a + s * (b - a)) </summary>
    /// <param name="a"> Начальное значение. </param>
    /// <param name="b"> Конечное значение. </param>
    /// <param name="s"> Процент интерполяции. 0 = a, 1 = b, 0.5 = середина между a и b. </param>
    public static double Lerp(double a, double b, double s)
    {
        return a + s * (b - a);
    }
}

Причины публикации

  • В HLSL языке нет стандартной функции Sin для double (но это не точно)
  • В интернете мало доступной информации по этой теме

Рассмотренные подходы

Анализируемые параметры

  • Точность по отношению к Math.Sin
  • Скорость по отношению к Math.Sin

Кроме анализа, мы улучшим их быстродействие.

Ряды Тейлора

Плюсы:

  • Высочайшая точность

    Эта функция, используемая для расчетов значения Sin, может использоваться для расчета бесконечно точного значения Sin. Чем больше итераций она претерпевает, тем точнее на выходе получается значение (в гипотезе). В практике программирования стоит учитывать ошибки округления вычислений в зависимости от используемых типов параметров (double, float, decimal и другие).
  • Рассчитывает любой угол

    В качестве аргумента в функцию можно внести любое значение, поэтому нет необходимости мониторить входные параметры.
  • Независимость

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

Минусы:

  • Очень низкая скорость (4-10%)

    Требует много итераций для того, чтобы точность приблизилась к точности Math.Sin, в результате чего работает в 25 раз медленнее стандартной функции.
  • Чем больше угол, тем ниже точность

    Чем больше угол введен в функцию, тем больше итераций нужно чтобы достичь той же точности, что и Math.Sin.

Оригинальный вид (скорость: 4%):

Стандартная функция подразумевает вычисления факториалов, а также возведение в степень каждую итерацию.

image

Доработанный (скорость: 10%):

Вычисление некоторых степеней можно сократить в цикле (a *= aa;), а другие факториалы предварительно вычислить и вынести в массив, при этом изменение знаков (+, -, +, ...) можно не возводить в степень и также сократить их вычисление, используя предыдущие значения.

Результатом получится следующий код:

Sin(rad, steps)

    // <summary> Массив факториалов, используемый в функции Fact </summary>
    static double[] fact;

    /// <summary> Рассчитывает с высокой точностью синус угла в радианах с помощью рядов Тейлора. 
    /// Чем больше rad, тем ниже точность. 
    /// Скорость (к Math): 4% (fps) при steps = 17 </summary>
    /// <param name="rad"> Угол в радианах. Лучше всего рассчитывать угол до pi/4. </param>
    /// <param name="steps"> Количество шагов: чем больше, тем точнее результат. Для pi/4 и точности E-15 достаточно 8. </param>
    public static double Sin(double rad, int steps)
    {
        double ret;
        double a;   //Угол, возведенный в нужную степень 
        double aa;  //Угол * угол
        int i_f;    //Индекс факториала
        int sign;   //Знак (колеблется от - до +, при этом первый раз = +)

        ret     = 0;
        sign    = -1;
        aa      = rad * rad;
        a       = rad;
        i_f     = 1;

        //Определение тригонометрических функций через ряды Тейлора
        for (int n = 0; n < steps; n++)
        {
            sign    *= -1;
            ret     += sign * a / Fact(i_f);
            a       *= aa;
            i_f     += 2;
        }

        return ret;
    }    
    /// <summary> Возвращает факториал (n!). Если n > fact.Length, возвращает -1. </summary>
    /// <param name="n"> Число, которое нужно возвести в факториал. </param>
    public static double Fact(int n)
    {
        if (n >= 0 && n < fact.Length)
        {
            return fact[n];
        }
        else
        {
            Debug.Log("Выход за пределы массива факториала. n = " + n + ", длина массива = " + fact.Length);
            return -1;
        }
    }
    /// <summary> Инициирует факториал. </summary>
    static void Init_Fact()
    {
        int steps;

        steps = 46;

        fact = new double[steps];

        fact[0] = 1;

        for (int n = 1; n < steps; n++)
        {
            fact[n] = fact[n - 1] * n;
        }
    }

Улучшенный вид (скорость: 19%):

Мы знаем, что чем меньше угол, тем меньше нужно итераций. Самый малый нужный нам угол = 0.25*PI, т.е. 45 градусов. Считая Sin и Cos в области 45 градусов мы можем получить все значения от -1 до 1 для Sin (в области 2 * PI). Для этого мы разделим круг (2*PI) на 8 частей и для каждой части укажем свой способ расчета синуса. Более того, чтобы ускорить расчет, откажемся от использования функции получения остатка (%) (для получения положения угла внутри зоны 45 градусов):

Sin_2(rad)

    // <summary> Массив факториалов, используемый в функции Fact </summary>
    static double[] fact;
    /// <summary> Надстройка над Sin </summary>
    /// <param name="rad"></param>
    public static double Sin_2(double rad)
    {
        double rad_025;
        int i;

        //rad = rad % PI2;    //% - довольно резурсозатратная функция. Если заменить, fps поднимется с 90 до 150 (при вызове 100 000 раз в кадр)
        //rad_025 = rad % PI025;
        i = (int)(rad / PI025);
        rad_025 = rad - PI025 * i;

        i = i & 7;  //Находим остаток от деления на 8 частей

        //Ищем индекс кусочка круга
        switch (i)
        {
            case 0:
                return Sin(rad_025, 8);
            case 1:
                return Cos(PI025 - rad_025, 8);
            case 2:
                return Cos(rad_025, 8);
            case 3:
                return Sin(PI025 - rad_025, 8);
            case 4:
                return -Sin(rad_025, 8);
            case 5:
                return -Cos(PI025 - rad_025, 8);
            case 6:
                return -Cos(rad_025, 8);
            case 7:
                return -Sin(PI025 - rad_025, 8);
        }

        return 0;
    }
    /// <summary> Рассчитывает с высокой точностью синус угла в радианах с помощью рядов Тейлора. 
    /// Чем больше rad, тем ниже точность. 
    /// Скорость (к Math): 10% (fps) при steps = 17 </summary>
    /// <param name="rad"> Угол в радианах. Лучше всего рассчитывать угол до pi/4. </param>
    /// <param name="steps"> Количество шагов: чем больше, тем точнее результат. Для pi/4 и точности E-15 достаточно 8. </param>
    public static double Sin(double rad, int steps)
    {
        double ret;
        double a;   //Угол, возведенный в нужную степень 
        double aa;  //Угол * угол
        int i_f;    //Индекс факториала
        int sign;   //Знак (колеблется от - до +, при этом первый раз = +)

        ret     = 0;
        sign    = -1;
        aa      = rad * rad;
        a       = rad;
        i_f     = 1;

        //Определение тригонометрических функций через ряды Тейлора
        for (int n = 0; n < steps; n++)
        {
            sign    *= -1;
            ret     += sign * a / Fact(i_f);
            a       *= aa;
            i_f     += 2;
        }

        return ret;
    }

    /// <summary> Рассчитывает с высокой точностью косинус угла в радианах с помощью рядов Тейлора. 
    /// Чем больше rad, тем ниже точность. 
    /// Скорость (к Math): 10% (fps), 26% (test) при steps = 17 </summary>
    /// <param name="rad"> Угол в радианах. Лучше всего рассчитывать угол до pi/4. </param>
    /// <param name="steps"> Количество шагов: чем больше, тем точнее результат. Для pi/4 и точности E-15 достаточно 8. </param>
    public static double Cos(double rad, int steps)
    {
        double ret;
        double a;
        double aa;  //Угол * угол
        int i_f;    //Индекс факториала
        int sign;   //Знак (колеблется от - до +, при этом первый раз = +)

        ret     = 0;
        sign    = -1;
        aa      = rad * rad;
        a       = 1;
        i_f     = 0;

        //Определение тригонометрических функций через ряды Тейлора
        for (int n = 0; n < steps; n++)
        {
            sign    *= -1;
            ret     += sign * a / Fact(i_f);
            a       *= aa;
            i_f     += 2;
        }

        return ret;
    }

    /// <summary> Возвращает факториал (n!). Если n > fact.Length, возвращает -1. </summary>
    /// <param name="n"> Число, которое нужно возвести в факториал. </param>
    public static double Fact(int n)
    {
        if (n >= 0 && n < fact.Length)
        {
            return fact[n];
        }
        else
        {
            Debug.Log("Выход за пределы массива факториала. n = " + n + ", длина массива = " + fact.Length);
            return -1;
        }
    }

    /// <summary> Инициирует факториал. </summary>
    static void Init_Fact()
    {
        int steps;

        steps = 46;

        fact = new double[steps];

        fact[0] = 1;

        for (int n = 1; n < steps; n++)
        {
            fact[n] = fact[n - 1] * n;
        }
    }

Полиномы

С данным способом я столкнулся на просторах интернета, автору нужна была быстрая функция поиска Sin для double пониженной точности (ошибка < 0.000 001) без использования библиотек предварительно вычисленных значений.

Плюсы:

  • Высокая скорость (9-84%)

    Изначально скинутый полином без изменений выдавал скорость в 9% от оригинального Math.Sin, что в 10 раз медленнее. Благодаря небольшим изменениям скорость резко возрастает до 84%, что неплохо, если закрыть глаза на точность.
  • Не требуется дополнительных предварительных вычислений и памяти

    Если вверху и далее внизу нам нужно составлять массивы переменных чтобы ускорить вычисления, то тут все ключевые коэффициенты были любезно рассчитаны и помещены в формулу самим автором в виде констант.
  • Точность выше чем у Mathf.Sin (float)

    Для сравнения:

    0.841471 — Mathf.Sin(1)(движок Unity);
    0.841470984807897 — Math.Sin(1)(стандартная функция C#);
    0.841470956802368 — sin(1)(GPU, язык hlsl);
    0.841471184637935Sin_0(1).

Минусы:

  • Не универсальна

    Нельзя настроить точность вручную потому что неизвестно каким инструментарием пользовался автор для вычисления данного полинома.
  • Зачем?

    Зачем автору понадобилась функция, которая не требует никаких массивов и у которой такая низкая (по сравнению с double) точность?

Оригинальный вид:

Sin_1(x)

/// <summary> Скорость (к Math): 9% (fps)</summary>
/// <param name="x"> Угол в радианах от -2*Pi до 2*Pi </param>
public static double Sin_1(double x)
{
    return
        0.9999997192673006 * x - 0.1666657564532464 * Math.Pow(x, 3) +
        0.008332803647181511 * Math.Pow(x, 5) - 0.00019830197237204295 * Math.Pow(x, 7) +
        2.7444305061093514e-6 * Math.Pow(x, 9) - 2.442176561869478e-8 * Math.Pow(x, 11) +
        1.407555708887347e-10 * Math.Pow(x, 13) - 4.240664814288337e-13 * Math.Pow(x, 15);
}

Улучшенный вид:

Sin_0(rad)

/// <summary> Скорость (к Math): 83% (fps)</summary>
/// <param name="rad"> Угол в радианах от -2*Pi до 2*Pi </param>
public static double Sin_0(double rad)
{
    double x;
    double xx;
    double ret;

    xx = rad * rad;
    x = rad;                      //1
    ret = 0.9999997192673006 * x;
    x *= xx;                        //3
    ret -= 0.1666657564532464 * x;
    x *= xx;                        //5
    ret += 0.008332803647181511 * x;
    x *= xx;                        //7
    ret -= 0.00019830197237204295 * x;
    x *= xx;                        //9
    ret += 2.7444305061093514e-6 * x;
    x *= xx;                        //11
    ret -= 2.442176561869478e-8 * x;
    x *= xx;                        //13
    ret += 1.407555708887347e-10 * x;
    x *= xx;                        //15
    ret -= 4.240664814288337e-13 * x;

    return ret;
}

Линейная интерполяция

Данный метод основан на линейной интерполяции между результатами двух записей в массиве.
Записи разделены на mem_sin и mem_cos, в них содержатся предварительно рассчитанные результаты стандартной функции Math.Sin и Math.Cos на отрезке входных параметров от 0 до 0.25*PI.

Принцип манипуляций с углом от 0 до 45 градусов не отличается от улучшенной версии рядов Тейлора, но при этом вызывается функция, которая находит — между какими двумя записями находится угол — и находит значение между ними.

Плюсы:

  • Высокая скорость (65%)

    Благодаря простоте алгоритма интерполяции, скорость достигает 65% от скорости Math.Sin. Считаю скорость >33% удовлетворительной.
  • Высочайшая точность

    Пример редкого случая отклонения:
    0.255835595715180 — Math.Sin;
    0.255835595715179Sin_3.
  • Быстрая нога

    Люблю эту функцию потому что она родилась в муках, написана мной и превзошла требования: скорость > 33%, точность выше 1е-14. Дам ей гордое имя — Vēlōx Pes.

Минусы:

  • Требует место в памяти

    Для работы необходимо предварительно вычислить два массива: для sin и для cos; каждый массив весит ~16мб (16*2=32мб)

Оригинальный вид:

Sin_3(rad)

    const double PI025      = Math.PI / 4;
    /// <summary> 2^17 = 131072 (1 мб), ошибка меньше 10000 (с конца), при 2^21 = 22097152 (16 мб) минимальная ошибка +-1 (с конца) (почти нет ошибки) </summary>
    const int length_mem    = 22097152;
    const double step_rad   = PI025 / length_mem;
    /// <summary> Массив предварительно рассчитанного sin, используемый для линейной интерполяции. </summary>
    static double[] mem_sin;
    /// <summary> Массив предварительно рассчитанного cos, используемый для линейной интерполяции. </summary>
    static double[] mem_cos;

    /// <summary> Инициирует данные, вроде массивов sin, необходимые в ходе расчетов. </summary>
    public static void Initialise()
    {
        Ini_Mem_Sin();
        Ini_Mem_Cos();
    }
    public static double Sin_3(double rad)
    {
        double rad_025;
        int i;

        i = (int)(rad / PI025);     //Находим индекс части
        rad_025 = rad - PI025 * i;  //Находим отступ от начала части (в радианах)
        i = i & 7;                  //Находим остаток от деления индекса на 8

        //Ищем индекс кусочка круга
        switch (i)
        {
            case 0:
                return Sin_Lerp(rad_025);
            case 1:
                return Cos_Lerp(PI025 - rad_025);
            case 2:
                return Cos_Lerp(rad_025);
            case 3:
                return Sin_Lerp(PI025 - rad_025);
            case 4:
                return -Sin_Lerp(rad_025);
            case 5:
                return -Cos_Lerp(PI025 - rad_025);
            case 6:
                return -Cos_Lerp(rad_025);
            case 7:
                return -Sin_Lerp(PI025 - rad_025);
        }

        return 0;
    }

    /// <summary> Подготавливает массив sin для линейной интерполяции </summary>
    static void Ini_Mem_Sin()
    {
        double rad;
        mem_sin = new double[length_mem];
        for (int i = 0; i < length_mem; i++)
        {
            rad = (i * PI025) / length_mem;
            mem_sin[i] = Math.Sin(rad);
        }
    }
    /// <summary> Подготавливает массив cos для линейной интерполяции </summary>
    static void Ini_Mem_Cos()
    {
        double rad;
        mem_cos = new double[length_mem];
        for (int i = 0; i < length_mem; i++)
        {
            rad = (i * PI025) / length_mem;
            mem_cos[i] = Math.Cos(rad);
        }
    }
    /// <summary> Использует линейную интерполяцию для поиска sin от 0 до pi/4. </summary>
    /// <param name="rad"> Угол в радианах от 0 до pi/4. </param>
    static double Sin_Lerp(double rad)
    {
        int i_0;
        int i_1;
        double i_0d;
        double percent;
        double a;
        double b;
        double s;

        percent = rad / PI025;
        i_0d    = percent * length_mem;
        i_0     = (int)i_0d;
        i_1     = i_0 + 1;
        a       = mem_sin[i_0];
        b       = mem_sin[i_1];
        s       = i_0d - i_0;

        return Math_d.Lerp(a, b, s);
    }
    /// <summary> Использует линейную интерполяцию для поиска cos от 0 до pi/4. </summary>
    /// <param name="rad"> Угол в радианах от 0 до pi/4. </param>
    static double Cos_Lerp(double rad)
    {
        int i_0;
        int i_1;
        double i_0d;
        double percent;
        double a;
        double b;
        double s;

        percent = rad / PI025;
        i_0d = percent * length_mem;
        i_0 = (int)i_0d;
        i_1 = i_0 + 1;
        a = mem_cos[i_0];
        b = mem_cos[i_1];
        s = i_0d - i_0;

        return Math_d.Lerp(a, b, s);
    }
    /// <summary> Производит линейную интерполяцию между двумя значениями. (return a + s * (b - a)) </summary>
    /// <param name="a"> Начальное значение. </param>
    /// <param name="b"> Конечное значение. </param>
    /// <param name="s"> Процент интерполяции. 0 = a, 1 = b, 0.5 = середина между a и b. </param>
    public static double Lerp(double a, double b, double s)
    {
        return a + s * (b - a);
    }

Автор: Normal_Mur

Источник

Поделиться

* - обязательные к заполнению поля