Предисловие
Мне понадобилось вычислять дугу с повышенной точностью на процессоре видеокарты в режиме реального времени.
Автор не ставил перед собой цель превзойти стандартную функцию 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); }
ссылка на оригинал статьи https://habr.com/post/426355/
Добавить комментарий