Math.Sin (double) function for GPU
Foreword
I needed to calculate the arc with increased accuracy on the video card processor in real time.
The author did not set himself the goal of exceeding the standard function System.Math.Sin () (C #) and did not achieve it.
The result of the work and my choice (for those who do not want to read):
Sin_3 (rad)
using System;
classMath_d
{
constdouble PI025 = Math.PI / 4;
///<summary> 2^17 = 131072 (1 мб), ошибка меньше 10000 (с конца), при 2^21 = 22097152 (16 мб) минимальная ошибка +-1 (с конца) (почти нет ошибки) </summary>constint length_mem = 22097152;
constint length_mem_M1 = length_mem - 1;
///<summary> Массив предварительно рассчитанного sin, используемый для линейной интерполяции. </summary>staticdouble[] mem_sin;
///<summary> Массив предварительно рассчитанного cos, используемый для линейной интерполяции. </summary>staticdouble[] mem_cos;
///<summary> Инициирует данные, вроде массивов sin, необходимые в ходе расчетов. </summary>publicstaticvoidInitialise()
{
Ini_Mem_Sin();
Ini_Mem_Cos();
}
///<summary> Использует линейную интерполяцию для поиска значения Cos, записанные в массив. </summary>///<param name="rad">Угол</param>publicstaticdoubleSin_3(double rad)
{
double rad_025;
int i;
//Обрабатываем отрицательные значения угла//if (rad < 0) { rad = -rad + Math.PI; }
i = (int)(rad / PI025); //Находим индекс части
rad_025 = rad - PI025 * i; //Находим отступ от начала части (в радианах)
i = i & 7; //Находим остаток от деления индекса на 8//Ищем индекс кусочка кругаswitch (i)
{
case0:
return Sin_Lerp(rad_025);
case1:
return Cos_Lerp(PI025 - rad_025);
case2:
return Cos_Lerp(rad_025);
case3:
return Sin_Lerp(PI025 - rad_025);
case4:
return -Sin_Lerp(rad_025);
case5:
return -Cos_Lerp(PI025 - rad_025);
case6:
return -Cos_Lerp(rad_025);
case7:
return -Sin_Lerp(PI025 - rad_025);
}
return0;
}
///<summary> Подготавливает массив sin для линейной интерполяции </summary>staticvoidIni_Mem_Sin()
{
double rad;
mem_sin = newdouble[length_mem];
for (int i = 0; i < length_mem; i++)
{
rad = (i * PI025) / length_mem_M1;
mem_sin[i] = Math.Sin(rad);
}
}
///<summary> Подготавливает массив cos для линейной интерполяции </summary>staticvoidIni_Mem_Cos()
{
double rad;
mem_cos = newdouble[length_mem];
for (int i = 0; i < length_mem; i++)
{
rad = (i * PI025) / length_mem_M1;
mem_cos[i] = Math.Cos(rad);
}
}
///<summary> Использует линейную интерполяцию для поиска sin от 0 до pi/4. </summary>///<param name="rad"> Угол в радианах от 0 до pi/4. </param>staticdoubleSin_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_M1;
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 Lerp(a, b, s);
}
///<summary> Использует линейную интерполяцию для поиска cos от 0 до pi/4. </summary>///<param name="rad"> Угол в радианах от 0 до pi/4. </param>staticdoubleCos_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_M1;
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 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>publicstaticdoubleLerp(double a, double b, double s)
{
return a + s * (b - a);
}
}
Reasons for posting
- There is no standard Sin for double function in the HLSL language (but this is not accurate)
- On the Internet there is little information available on this topic.
Considered approaches
- Taylor Rows (Wikipedia)
- Polynomials ( author of the function: "asvp" )
- Linear interpolation of pre-calculated Math.Sin and Math.Cos results (author: I)
Analyzed parameters
- Accuracy with respect to Math.Sin
- Speed relative to Math.Sin
In addition to the analysis, we will improve their speed.
Taylor Rows
Pros:
- Highest accuracyЭта функция, используемая для расчетов значения Sin, может использоваться для расчета бесконечно точного значения Sin. Чем больше итераций она претерпевает, тем точнее на выходе получается значение (в гипотезе). В практике программирования стоит учитывать ошибки округления вычислений в зависимости от используемых типов параметров (double, float, decimal и другие).
- Calculates any angleВ качестве аргумента в функцию можно внести любое значение, поэтому нет необходимости мониторить входные параметры.
- IndependenceНе требует предварительных вычислений (как функции, рассмотренные ниже), и часто является базой, на которой собираются более быстрые функции.
Minuses:
- Very low speed (4-10%)Требует много итераций для того, чтобы точность приблизилась к точности Math.Sin, в результате чего работает в 25 раз медленнее стандартной функции.
- The greater the angle, the lower the accuracy.Чем больше угол введен в функцию, тем больше итераций нужно чтобы достичь той же точности, что и Math.Sin.
Original form (speed: 4%):
The standard function involves calculating factorials, as well as raising to a power each iteration.
Modified (speed: 10%):
Calculation of some degrees can be reduced in a cycle (a * = aa;), and other factorials can be pre-calculated and put into an array, while changing the signs (+, -, +, ...) can not to raise to a power and also reduce their calculation using the previous values.
The result is the following code:
Sin (rad, steps)
// <summary> Массив факториалов, используемый в функции Fact </summary>staticdouble[] fact;
///<summary> Рассчитывает с высокой точностью синус угла в радианах с помощью рядов Тейлора. /// Чем больше rad, тем ниже точность. /// Скорость (к Math): 4% (fps) при steps = 17 </summary>///<param name="rad"> Угол в радианах. Лучше всего рассчитывать угол до pi/4. </param>///<param name="steps"> Количество шагов: чем больше, тем точнее результат. Для pi/4 и точности E-15 достаточно 8. </param>publicstaticdoubleSin(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>publicstaticdoubleFact(int n)
{
if (n >= 0 && n < fact.Length)
{
return fact[n];
}
else
{
Debug.Log("Выход за пределы массива факториала. n = " + n + ", длина массива = " + fact.Length);
return-1;
}
}
///<summary> Инициирует факториал. </summary>staticvoidInit_Fact()
{
int steps;
steps = 46;
fact = newdouble[steps];
fact[0] = 1;
for (int n = 1; n < steps; n++)
{
fact[n] = fact[n - 1] * n;
}
}
Improved view (speed: 19%):
We know that the smaller the angle, the less iteration is needed. The smallest angle we need is 0.25 * PI, i.e. 45 degrees. Considering Sin and Cos in the area of 45 degrees, we can get all the values from -1 to 1 for Sin (in the area of 2 * PI). To do this, we divide the circle (2 * PI) into 8 parts and for each part we indicate our own method of calculating the sine. Moreover, in order to speed up the calculation, we will abandon the use of the function for obtaining the residue (%) (to obtain the position of the angle inside the 45 degree zone):
Sin_2 (rad)
// <summary> Массив факториалов, используемый в функции Fact </summary>staticdouble[] fact;
///<summary> Надстройка над Sin </summary>///<param name="rad"></param>publicstaticdoubleSin_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)
{
case0:
return Sin(rad_025, 8);
case1:
return Cos(PI025 - rad_025, 8);
case2:
return Cos(rad_025, 8);
case3:
return Sin(PI025 - rad_025, 8);
case4:
return -Sin(rad_025, 8);
case5:
return -Cos(PI025 - rad_025, 8);
case6:
return -Cos(rad_025, 8);
case7:
return -Sin(PI025 - rad_025, 8);
}
return0;
}
///<summary> Рассчитывает с высокой точностью синус угла в радианах с помощью рядов Тейлора. /// Чем больше rad, тем ниже точность. /// Скорость (к Math): 10% (fps) при steps = 17 </summary>///<param name="rad"> Угол в радианах. Лучше всего рассчитывать угол до pi/4. </param>///<param name="steps"> Количество шагов: чем больше, тем точнее результат. Для pi/4 и точности E-15 достаточно 8. </param>publicstaticdoubleSin(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>publicstaticdoubleCos(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>publicstaticdoubleFact(int n)
{
if (n >= 0 && n < fact.Length)
{
return fact[n];
}
else
{
Debug.Log("Выход за пределы массива факториала. n = " + n + ", длина массива = " + fact.Length);
return-1;
}
}
///<summary> Инициирует факториал. </summary>staticvoidInit_Fact()
{
int steps;
steps = 46;
fact = newdouble[steps];
fact[0] = 1;
for (int n = 1; n < steps; n++)
{
fact[n] = fact[n - 1] * n;
}
}Polynomials
I encountered this method on the Internet, the author needed a quick search function Sin for doubles of lower accuracy (error <0.000 001) without using libraries of previously calculated values.
Pros:
- High speed (9-84%)Изначально скинутый полином без изменений выдавал скорость в 9% от оригинального Math.Sin, что в 10 раз медленнее. Благодаря небольшим изменениям скорость резко возрастает до 84%, что неплохо, если закрыть глаза на точность.
- No additional preliminary calculations and memory required.Если вверху и далее внизу нам нужно составлять массивы переменных чтобы ускорить вычисления, то тут все ключевые коэффициенты были любезно рассчитаны и помещены в формулу самим автором в виде констант.
- Accuracy is higher than 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).
Minuses:
- Not universalНельзя настроить точность вручную потому что неизвестно каким инструментарием пользовался автор для вычисления данного полинома.
- What for?Зачем автору понадобилась функция, которая не требует никаких массивов и у которой такая низкая (по сравнению с double) точность?
Original appearance:
Sin_1 (x)
///<summary> Скорость (к Math): 9% (fps)</summary>///<param name="x"> Угол в радианах от -2*Pi до 2*Pi </param>publicstaticdoubleSin_1(double x)
{
return0.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);
}
Superior view:
Sin_0 (rad)
///<summary> Скорость (к Math): 83% (fps)</summary>///<param name="rad"> Угол в радианах от -2*Pi до 2*Pi </param>publicstaticdoubleSin_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;
}
Linear interpolation
This method is based on linear interpolation between the results of two records in an array.
The records are divided into mem_sin and mem_cos, they contain the pre-calculated results of the standard function Math.Sin and Math.Cos on a segment of input parameters from 0 to 0.25 * PI.
The principle of manipulation with an angle from 0 to 45 degrees does not differ from the improved version of the Taylor series, but it calls a function that finds — between which two records there is an angle — and finds a value between them.
Pros:
- High speed (65%)Благодаря простоте алгоритма интерполяции, скорость достигает 65% от скорости Math.Sin. Считаю скорость >33% удовлетворительной.
- Highest accuracyПример редкого случая отклонения:
0.255835595715180 — Math.Sin;
0.255835595715179 — Sin_3. - Fast legЛюблю эту функцию потому что она родилась в муках, написана мной и превзошла требования: скорость > 33%, точность выше 1е-14. Дам ей гордое имя — Vēlōx Pes.
Minuses:
- Requires memory spaceДля работы необходимо предварительно вычислить два массива: для sin и для cos; каждый массив весит ~16мб (16*2=32мб)
Original appearance:
Sin_3 (rad)
classMath_d
{
constdouble PI025 = Math.PI / 4;
///<summary> 2^17 = 131072 (1 мб), ошибка меньше 10000 (с конца), при 2^21 = 22097152 (16 мб) минимальная ошибка +-1 (с конца) (почти нет ошибки) </summary>constint length_mem = 22097152;
constint length_mem_M1 = length_mem - 1;
///<summary> Массив предварительно рассчитанного sin, используемый для линейной интерполяции. </summary>staticdouble[] mem_sin;
///<summary> Массив предварительно рассчитанного cos, используемый для линейной интерполяции. </summary>staticdouble[] mem_cos;
///<summary> Инициирует данные, вроде массивов sin, необходимые в ходе расчетов. </summary>publicstaticvoidInitialise()
{
Ini_Mem_Sin();
Ini_Mem_Cos();
}
///<summary> Использует линейную интерполяцию для поиска значения Cos, записанные в массив. </summary>///<param name="rad">Угол</param>publicstaticdoubleSin_3(double rad)
{
double rad_025;
int i;
//Обрабатываем отрицательные значения угла//if (rad < 0) { rad = -rad + Math.PI; }
i = (int)(rad / PI025); //Находим индекс части
rad_025 = rad - PI025 * i; //Находим отступ от начала части (в радианах)
i = i & 7; //Находим остаток от деления индекса на 8//Ищем индекс кусочка кругаswitch (i)
{
case0:
return Sin_Lerp(rad_025);
case1:
return Cos_Lerp(PI025 - rad_025);
case2:
return Cos_Lerp(rad_025);
case3:
return Sin_Lerp(PI025 - rad_025);
case4:
return -Sin_Lerp(rad_025);
case5:
return -Cos_Lerp(PI025 - rad_025);
case6:
return -Cos_Lerp(rad_025);
case7:
return -Sin_Lerp(PI025 - rad_025);
}
return0;
}
///<summary> Подготавливает массив sin для линейной интерполяции </summary>staticvoidIni_Mem_Sin()
{
double rad;
mem_sin = newdouble[length_mem];
for (int i = 0; i < length_mem; i++)
{
rad = (i * PI025) / length_mem_M1;
mem_sin[i] = Math.Sin(rad);
}
}
///<summary> Подготавливает массив cos для линейной интерполяции </summary>staticvoidIni_Mem_Cos()
{
double rad;
mem_cos = newdouble[length_mem];
for (int i = 0; i < length_mem; i++)
{
rad = (i * PI025) / length_mem_M1;
mem_cos[i] = Math.Cos(rad);
}
}
///<summary> Использует линейную интерполяцию для поиска sin от 0 до pi/4. </summary>///<param name="rad"> Угол в радианах от 0 до pi/4. </param>staticdoubleSin_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_M1;
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 Lerp(a, b, s);
}
///<summary> Использует линейную интерполяцию для поиска cos от 0 до pi/4. </summary>///<param name="rad"> Угол в радианах от 0 до pi/4. </param>staticdoubleCos_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_M1;
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 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>publicstaticdoubleLerp(double a, double b, double s)
{
return a + s * (b - a);
}
}
UPD: Fixed the error in determining the index in Sin_Lerp (), Cos_Lerp (), Ini_Mem_Sin () and Ini_Mem_Cos ().