Аппроксимация обратных тригонометрических функций - PullRequest
10 голосов
/ 11 сентября 2011

Я должен реализовать asin, acos и atan в среде, где у меня есть только следующие математические инструменты:

  • синус
  • косинус
  • элементарная арифметика с фиксированной запятой (числа с плавающей запятой недоступны)

У меня также уже есть достаточно хорошая функция квадратного корня.

Можно ли использовать их для реализации достаточно эффективных обратных тригонометрических функций?

Мне не нужна слишком большая точность (числа с плавающей запятой в любом случае имеют очень ограниченную точность), подойдет базовое приближение.

Я уже наполовину решил заняться поиском в таблице, но я хотел бы знать, есть ли какой-нибудь более аккуратный вариант (для которого не требуется несколько сотен строк кода только для реализации базовой математики).

EDIT:

Чтобы прояснить ситуацию : мне нужно запускать функцию сотни раз за кадр со скоростью 35 кадров в секунду.

Ответы [ 9 ]

8 голосов
/ 11 сентября 2011

В среде с фиксированной запятой (S15.16) я успешно использовал алгоритм CORDIC (общее описание см. В Википедии) для вычисления atan2 (y, x), затем вывел asin () и acos () из этого, используя wellИзвестные функциональные тождества, которые включают квадратный корень:

asin(x) = atan2 (x, sqrt ((1.0 + x) * (1.0 - x)))
acos(x) = atan2 (sqrt ((1.0 + x) * (1.0 - x)), x)

Оказывается, что найти полезное описание ИСПРАВЛЕНИЯ CORDIC для atan2 () на double сложнее, чем я думал.Следующий веб-сайт содержит достаточно подробное описание, а также обсуждает два альтернативных подхода, полиномиальную аппроксимацию и справочные таблицы:

http://ch.mathworks.com/examples/matlab-fixed-point-designer/615-calculate-fixed-point-arctangent

2 голосов
/ 31 декабря 2012

Должно быть легко добавить следующий код в фиксированную точку.В нем используется рациональное приближение для вычисления арктангенса, нормализованного к интервалу [0 1) (вы можете умножить его на Pi / 2, чтобы получить реальный арктангенс).Затем вы можете использовать хорошо известные тождества , чтобы получить арксинус / арккос из арктангенса.

normalized_atan(x) ~ (b x + x^2) / (1 + 2 b x + x^2)

where b = 0.596227

Максимальная ошибка составляет 0,1620º

#include <stdint.h>
#include <math.h>

// Approximates atan(x) normalized to the [-1,1] range
// with a maximum error of 0.1620 degrees.

float norm_atan( float x )
{
    static const uint32_t sign_mask = 0x80000000;
    static const float b = 0.596227f;

    // Extract the sign bit
    uint32_t ux_s  = sign_mask & (uint32_t &)x;

    // Calculate the arctangent in the first quadrant
    float bx_a = ::fabs( b * x );
    float num = bx_a + x * x;
    float atan_1q = num / ( 1.f + bx_a + num );

    // Restore the sign bit
    uint32_t atan_2q = ux_s | (uint32_t &)atan_1q;
    return (float &)atan_2q;
}

// Approximates atan2(y, x) normalized to the [0,4) range
// with a maximum error of 0.1620 degrees

float norm_atan2( float y, float x )
{
    static const uint32_t sign_mask = 0x80000000;
    static const float b = 0.596227f;

    // Extract the sign bits
    uint32_t ux_s  = sign_mask & (uint32_t &)x;
    uint32_t uy_s  = sign_mask & (uint32_t &)y;

    // Determine the quadrant offset
    float q = (float)( ( ~ux_s & uy_s ) >> 29 | ux_s >> 30 ); 

    // Calculate the arctangent in the first quadrant
    float bxy_a = ::fabs( b * x * y );
    float num = bxy_a + y * y;
    float atan_1q =  num / ( x * x + bxy_a + num );

    // Translate it to the proper quadrant
    uint32_t uatan_2q = (ux_s ^ uy_s) | (uint32_t &)atan_1q;
    return q + (float &)uatan_2q;
} 

В случаевам нужна большая точность, есть рациональная функция 3-го порядка:

normalized_atan(x) ~ ( c x + x^2 + x^3) / ( 1 + (c + 1) x + (c + 1) x^2 + x^3)

where c = (1 + sqrt(17)) / 8

, которая имеет максимальную погрешность аппроксимации 0,00811º

2 голосов
/ 12 сентября 2011

Вам нужна большая точность для функции arcsin(x)? Если нет, вы можете вычислить arcsin в N узлах и сохранить значения в памяти. Я предлагаю использовать линейное приближение. если x = A*x_(N) + (1-A)*x_(N+1), то x = A*arcsin(x_(N)) + (1-A)*arcsin(x_(N+1)), где известно arcsin(x_(N)).

2 голосов
/ 11 сентября 2011

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

например:
arcsin(z) = Sigma((2n!)/((2^2n)*(n!)^2)*((z^(2n+1))/(2n+1))) где n в [0,бесконечность)

1 голос
/ 25 сентября 2014

Отправляя сюда мой ответ от этого другого похожего вопроса.

У nVidia есть несколько замечательных ресурсов, которые я использовал для своих собственных нужд, несколько примеров: acos asin atan2 и т. Д. И т. Д.

Эти алгоритмы дают достаточно точные результаты. Вот прямой пример Python с вставленной в него копией кода:

import math
def nVidia_acos(x):
    negate = float(x<0)
    x=abs(x)
    ret = -0.0187293
    ret = ret * x
    ret = ret + 0.0742610
    ret = ret * x
    ret = ret - 0.2121144
    ret = ret * x
    ret = ret + 1.5707288
    ret = ret * math.sqrt(1.0-x)
    ret = ret - 2 * negate * ret
    return negate * 3.14159265358979 + ret

А вот результаты для сравнения:

nVidia_acos(0.5)  result: 1.0471513828611643
math.acos(0.5)    result: 1.0471975511965976

Это довольно близко! Умножьте на 57.29577951, чтобы получить результаты в градусах, что также соответствует их формуле "градусы".

1 голос
/ 11 сентября 2011

Expression as definite integrals

http://en.wikipedia.org/wiki/Inverse_trigonometric_functions#Expression_as_definite_integrals

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

Infinite series

1 голос
/ 11 сентября 2011

Может быть, какая-то интеллектуальная грубая сила, такая как Ньютон Рапсон.

Так что для решения asin () вы идете с самым крутым спадом на грех ()

0 голосов
/ 19 ноября 2015

Только непрерывные функции аппроксимируются полиномами. И arcsin (x) разрывается в точке x = 1. то же самое arccos (x). Но сокращение диапазона до интервала 1, sqrt (1/2) в этом случае позволяет избежать этой ситуации. У нас есть arcsin (x) = pi / 2-arccos (x), arccos (x) = pi / 2-arcsin (x). Вы можете использовать matlab для минимаксного приближения. Приблизительно только в диапазоне [0, sqrt (1/2) )] (если угол для этого arcsin равен запросу больше, чем sqrt (1/2), найдите функцию cos (x) .arctangent только для x <1.arctan (x) = pi / 2-arctan (1 / x). </p>

0 голосов
/ 12 сентября 2011

Используйте полиномиальное приближение.Подбор наименьших квадратов проще всего (в Microsoft Excel есть), а приближение Чебышева является более точным.

Этот вопрос уже рассматривался ранее: Как работают тригонометрические функции?

...