Получение приближенных значений тригонометрических функций быстро. (комментарии) (3 стр) / Общее / Форум / Программирование игр / GameDev.ru

}:+()___ [Smile]
А ещё таблицы обычно хуже SIMD'ифицируются.
Я вон предлагал желающим обогнать полином таблицей:
https://gamedev.ru/code/forum/?id=204567&page=9&m=5202338#m121

> Важно, что нахождение коэффициентов аппроксимации — задача однократная, поэтому можно хоть неделю молотить тупым перебором, лишь бы в рантайме бесплатно получить увеличение точности.
Э, так можно далеко зайти.

FordPerfect

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

Ты прав, я тоже об этом думал. Дело в том, что добавление одного слагаемого к разложению Чебышёва обычно сразу увеличивает точность на 1-2 десятичных порядка, и на этом фоне уменьшение максимальной погрешности от перехода к равномерному приближению выглядит не очень заметно. Кроме того, есть ещё понятие среднеквадратической погрешности, так вот она изменится ещё менее заметно (если вообще уменьшится). Но если делать суперважную библиотеку, то может и стоит потратить на это время.

Тут ещё что важно. Как я уже сказал, критерий точности на отрезке может быть определён по-разному, т.е. нельзя однозначно равномерный критерий считать лучше других (но как ориентир - пожалуй). И, наконец, эти нюансы сводятся почти на нет при округлениях промежуточных значений и конечного результата - ведь точность аппроксимации берётся всегда с запасом.

На практике имеет смысл совмещать и использовать сплайны невысокого порядка

Я это и имел в виду, говоря об интерполяции. Применительно к таблицам сплайны принято называть интерполяцией (она бывает не только линейной). Возможно ошибаюсь, но, по-моему, именно так реализованы стандартные библиотечные функции в C++.

Mirrel

насколько это будет производительнее решения через треугольник?

Не совсем понял вопрос. Проще померить производительность в конкретном коде.

Aenigma
Это точно мне был ответ?

> Возможно ошибаюсь, но, по-моему, именно так реализованы стандартные библиотечные функции в C++.
Зависит, какая реализация библиотеки, и какая функция.
exp2f от GCC, таблица+полином:
https://gamedev.ru/flame/forum/?id=234120&page=3&m=4704472#m33
Для синуса - вроде видел (в GCC? fdlibm?) чисто полином.
В Intel Math Kernel Library, вроде, для синуса таблица.

Да, в http://www.netlib.org/fdlibm/k_sin.c - полином.

Aenigma
> Но если делать суперважную библиотеку, то может и стоит потратить на это время.
Я себе для \(L_\infty\) оптимизации наваял программку, которая итеративно делает взвешенную \(L_2\) оптимизацию с автоматическим выбором веса точек: если точка вылезла за допустимый диапазон погрешности — то увеличиваем ей вес, повторяем пока все точки не влезли в коридор. Еще сверху автоматический выбор ширины коридора прикрутил. Очевидно, что алгоритм оставляет желать лучшего, однако достаточно простой и позволяет раскладывать на любой набор базисных функций, а не только одномерные полиномы.

> Как я уже сказал, критерий точности на отрезке может быть определён по-разному
Да, я в своей программе, например, могу задавать вес ошибки и искать минимумы как абсолютной, так и относительной погрешности.

> И, наконец, эти нюансы сводятся почти на нет при округлениях промежуточных значений и конечного результата - ведь точность аппроксимации берётся всегда с запасом.
Если цель оптимизации — это скорость, то точность обычно много меньше максимально возможной и нюансы округления не важны. Когда же стоит цель написать библиотечную функцию с максимальной точностью, то как раз занимаются перебором битов в окрестности оптимума с целью найти идеальный набор коэффициентов.

> Я это и имел в виду, говоря об интерполяции.
Обычно под интерполяцией понимают аппроксимацию, которая проходит через контрольные точки, в отличие от сплайнов.

FordPerfect, да, я в последнем сообщении перепутал автора, ну да ладно...

}:+()___ [Smile]

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

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

Выкладываю ещё порцию кода.

// Вычисляет азимут точки (x, y). Аналог команды fpatan или функции atan2, но в случае неопределённости вида 0/0 возвращает 0
_declspec(naked) float _vectorcall Arctg2(float y, float x) // xmm0=y, xmm1=x
{
  static const float ct[8] = // Таблица констант
  {
     -1.0f,                  // -1
     1.59940588f,            // a0
     0.472401142f,           // a0*a2
     0.939805329f,           // a0*(a2+b1)
     1.82871687f,            // a0*a1
     2.36185217f,            // a0*(a1+b0)
     0.0151078226f,          // a0*a3
     0.0844737813f           // a0*(a3+b2)
  };
  _asm
  {
    unpcklps xmm1, xmm0      // xmm1 = y # x
    mov edx, offset ct       // В edx - адрес таблицы констант
    movaps xmm2, xmm1        // xmm2 = y # x
    movmskps ecx, xmm1       // Биты 1 и 0 ecx содержат знаковые биты y и x соответственно
    mulps xmm2, xmm2         // xmm1 = y^2 # x^2
    movshdup xmm3, xmm2      // xmm3 = y^2
    movups xmm4, [edx]       // xmm4 = a0*(a2+b1) # a0*a2 : a0 # -1
    comiss xmm2, xmm3        // cf=1, если |y|>|x|
    movlhps xmm1, xmm2       // xmm1 = y^2 # x^2 : y # x
    movups xmm3, [edx+16]    // xmm3 = a0*(a3+b2) # a0*a3 : a0*(a1+b0) # a0*a1
    pshufd xmm2, xmm1, 177   // xmm2 = x^2 # y^2 : x # y
    sbb eax, eax             // eax = -1, если |y|>|x|; и eax = 0, если |y|<=|x|
    divps xmm2, xmm1         // xmm2 = x^2/y^2 # y^2/x^2 : x/y # y/x
    jnb arctg_cont           // Перейти, если |y|<=|x|
    movshdup xmm2, xmm2      // При |y|>|x| рассмотреть обратные отношения x^2/y^2 и x/y
    mulss xmm2, xmm4         // и поменять знак величины x/y
      arctg_cont:            // xmm2 = y^2/x^2 : y/x, если |y|<=|x|; или  x^2/y^2 : x/y, если |y|>=|x|
    mov edx, 0x7F921FB6      // В edx - двоичное представление числа pi/2, сдвинутое влево
    ucomiss xmm2, xmm2       // pf=1, если неопределённость
    jp arctg_exit            // Вернуть y в случае неопределённости
    shl ecx, 31              // Знак x - в бите 31 ecx, знак y - в cf
    pshufd xmm0, xmm2, 170   // Заполнить xmm0 значением z=y^2/x^2 или z=x^2/y^2 - аргумент многочлена
    shufps xmm4, xmm4, 229   // xmm4/a0 = a2+b1 # a2 : 1 # 1
    rcr edx, 1               // edx = pi/2*sgn(y)
    mulps xmm3, xmm0         // xmm3/a0 = (a3+b2)*z # a3*z : (a1+b0)*z # a1*z
    addps xmm3, xmm4         // xmm3/a0 = (a3+b2)*z+(a2+b1) # a3*z+a2 : (a1+b0)*z+1 # a1*z+1
    and eax, edx             // eax = pi/2*sgn(y), если |y|>|x|, иначе eax = 0
    mulps xmm0, xmm0         // xmm0 = z^2 # z^2 : z^2 # z^2
    movhlps xmm1, xmm3       // xmm1/a0 = (a3+b2)*z+(a2+b1) # a3*z+a2
    movd xmm4, edx           // xmm4 = pi/2*sgn(y)
    mulps xmm0, xmm1         // xmm0/a0 = (a3+b2)*z^3+(a2+b1)*z^2 # a3*z^3+a2*z^2
    addps xmm0, xmm3         // xmm0/a0 = (a3+b2)*z^3+(a2+b1)*z^2+(a1+b0)*z+1 # a3*z^3+a2*z^2+a1*z+1
    xor eax, ecx             // eax = pi/2*sgn(x/y), если |y|>|x|, иначе eax = 0
    movshdup xmm1, xmm0      // xmm0 = P; xmm1 = Q
    movd xmm3, eax           // xmm3 = pi/2*sgn(x/y), если |y|>|x|, иначе xmm3 = 0
    divss xmm0, xmm1         // xmm0 = P/Q
    jecxz arctg_low          // Если x>0, добавлять +-pi не надо
    addss xmm4, xmm4         // xmm4 = pi*sgn(y)
    addss xmm3, xmm4         // Учесть добавку +-pi для II и III четверти
      arctg_low:             // Завершающие действия
    mulss xmm0, xmm2         // xmm0 = y/x*P/Q, если |y|<=|x|; иначе -x/y*P/Q 
    addss xmm0, xmm3         // xmm0 = arctg(y/x) с учётом квадранта
      arctg_exit:            // Метка для выхода из функции
    ret                      // Вернуться
  }
}

// Арксинус
_declspec(naked) float _vectorcall Arcsin(float x)
{
  static const float One = 1;
  _asm
  {
    movss xmm2, One        // xmm2 = 1
    movaps xmm1, xmm0      // xmm1 = x
    addss xmm1, xmm2       // xmm1 = 1+x
    subss xmm2, xmm0       // xmm2 = 1-x
    mulss xmm1, xmm2       // xmm1 = (1-x)(1+x) = 1-x^2
    sqrtss xmm1, xmm1      // xmm1 = (1-x^2)^0.5
    jmp Arctg2             // xmm0 = arctg(x/(1-x^2)^0.5) = arcsin x
  }
}

// Арккосинус
_declspec(naked) float _vectorcall Arccos(float x)
{
  static const float One = 1;
  _asm
  {
    movss xmm2, One        // xmm2 = 1
    movaps xmm1, xmm0      // xmm1 = x
    addss xmm0, xmm2       // xmm0 = 1+x
    subss xmm2, xmm1       // xmm2 = 1-x
    mulss xmm0, xmm2       // xmm0 = (1-x)(1+x) = 1-x^2
    sqrtss xmm0, xmm0      // xmm0 = (1-x^2)^0.5
    jmp Arctg2             // xmm0 = arcctg(x/(1-x^2)^0.5) = arccos x
  }
}

inline float XMScalarSinEst(float Value) noexcept
{
    // Map Value to y in [-pi,pi], x = 2*pi*quotient + remainder.
    float quotient = XM_1DIV2PI * Value;
    if (Value >= 0.0f)
    {
        quotient = static_cast<float>(static_cast<int>(quotient + 0.5f));
    }
    else
    {
        quotient = static_cast<float>(static_cast<int>(quotient - 0.5f));
    }
    float y = Value - XM_2PI * quotient;

    // Map y to [-pi/2,pi/2] with sin(y) = sin(Value).
    if (y > XM_PIDIV2)
    {
        y = XM_PI - y;
    }
    else if (y < -XM_PIDIV2)
    {
        y = -XM_PI - y;
    }

    // 7-degree minimax approximation
    float y2 = y * y;
    return (((-0.00018524670f * y2 + 0.0083139502f) * y2 - 0.16665852f) * y2 + 1.0f) * y;
}

Aenigma
> Не совсем понял вопрос. Проще померить производительность в конкретном коде.
Любой угол - это обычные координаты. Нам зачастую ни синус ни косинус не даны. Координаты мы знаем практически всегда.
И по принципу треугольника можно просто вычислять отношение противолежащего катета к гипотенузе, прилежащего катета к гипотенузе, противолежащего катета к прилежащему и противолежащего катета к прилежащему.

Точнее сразу значение, без вычисления тригонометрических функций.

В коде разложение тригонометрических ф-ций в ряд Тейлора?

Mirrel
Область применения тригонометрических функций выходит далеко за пределы задач по геометрии. Кроме того, если исходные данные представляют собой координаты, это совсем другая задача, сравнение тут неуместно.

Skunk, конечно, нет. Ряд Тейлора годится только для аппроксимации в небольшой окрестности точки. Прочитайте мои сообщения на предыдущей странице.