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

2018-10-14 в 15:08, admin, рубрики: 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 (но это не точно)
В интернете мало доступной информации по этой теме

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

Ряды Тейлора (Википедия)
Полиномы (автор функции: «asvp»)
Линейная интерполяция предварительно рассчитанных результатов Math.Sin и Math.Cos (автор: я)

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

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

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

Ряды Тейлора

Плюсы:

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

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

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

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

Минусы:

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

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

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

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

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

Доработанный (скорость: 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.841471184637935 — Sin_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.255835595715179 — Sin_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

Источник

Информация

Обсуждаемое

Рекомендуем

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

Предисловие

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

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

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

Ряды Тейлора

Полиномы

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

Архив

Информация

Обсуждаемое

Рекомендуем

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

Предисловие

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

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

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

Ряды Тейлора

Полиномы

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

Рекомендованный контент

Новости

Актуальные темы

Архив