быстрая тригонометрия , ускорить вычисление arcsin/cos/tan

Господа! Помогите придумать алгоритм быстрого вычисления обратных тригонометрических функций (именно - arcsin, arccos, arctan). Результат нужен с точностью до градуса.
С синусами и косинусами я разобрался - вогнал их значения в массив длиной 360, потом при обращении к моей функции просто беру значение из массива (прирост производительности сильно ощутим).
Например fastsin(30) возвращает значение sintable[30]. C обратными функциями все гораздо сложнее. Поиск в массиве простым перебором работает раз в 1000 (а то и более) медленнее, чем стандартные более точные аналоги.

А значения для обратных функций у тебя только твои (из sintable[]) или любые?

Можно сделать таблицу обратных значений. Сложности возникнут только для аргументов близких к +/-1. Для аргументов больших sqrt(2)/2 пользуйся формулами приведения y=sqrt(1-x^2)

Можно линеаризовать функцию. Можно представить косинус как y=1-x^2/2, откуда x^2=2(1-y) и x=sqrt(2(1-y)) и составить таблицу для функции arccos(1-x^2/2), тогда зная косинус=y, вычисляешь x=sqrt(2(1-y)) и используя его как индекс, ищешь нужное значение.

Но первый способ проще.

Да и в твоем алгоритме поиск лучше осуществлять бинарный а не линейный, у тебя же монотоннае функции.

Проще не значит лучше. Будем пробовать. amk, спасибо за инфу.
Будут еще у кого какие мысли - кидайте, буду только рад.

Guest, значения можно брать из sintable, но не обязательно, в зависимости от алгоритма, который предложишь.

Если Вас не смущают большие массивы, то для arccos, например, тоже сделайте таблицу.
Так как h=cos(0)-cos(1 град)=0.0001523..., то понадобится массив длины 2/h=13132

Цитата Го_ @ 05.01.06, 09:43

Так как h=cos(0)-cos(1 град)=0.0001523..., то понадобится массив длины 2/h=13132

А разве шаг будет постоянный?
Может стоит найти граничные значения аргумента для каждого градуса, т.е. когда аркфункция будет равна 0.5, 1.5, 2.5 и т.д. и хранить именно их в таблице? Дальше при задании аргумента искать бинарным поиском промежуток, в который попало значение и доставать результат для этого промежутка. Должно быть быстро.

FlyNN, скажите, откуда у вас такая задача возникла?

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

MeG, для игры. Движок - говно, вот и приходится изуверством заниматься. Я вообще раньше не парился по поводу производительности. А сейчас вот занялся и много нового для себя открыл.

Посмотри "формула Тейлора для многочлена" - по сути представление фуекции в виде многочлена.

Pn(x)=f(x0) + f1(x0)/1!*(x-x0) + f2(x0)/2!*((x-x0)^2) + f3(x0)/3!*((x-x0)^3) + ... + fn(x0)/n!*((x-x0)^n)+Rx
 
Rx=fk(x0)/n!*((x-x0)^n), k=n+1

Rx - остаточный член
fn - n-ая производная от функции
n - натуральное число, определяет точность функции. Такой способ заложен вроде в калькуляторах, но если уменьшить точность, то...

FlyNN, возьми интеловскую математическую библиотеку заточенную под математическое аппаратное обеспечение и получишь наивысшее сочетание точности и производительности

Где ее взять?

CORDIC is an acronym for COrdinate Rotation DIgital Computer. It is a class of shift-add algorithms for rotating vectors in a plane. In a nutshell, the CORDIC rotator performs a rotation using a series of specific incremental rotation angles selected so that each is performed by a shift and add operation.

Rotation of unit vectors provides us with a way to accurately compute trig functions, as well as a mechanism for computing the magnitude and phase angle of an input vector. Vector rotation is also useful in a host of DSP applications including modulation and Fourier Transforms.

http://www.fpga-guru.com/cordic.htm
Там есть ссылка на pdf файл, в котором описан алгоритм.

Цитата FlyNN @ 12.01.06, 12:57

Где ее взять?

Самый простой способ, на www.intel.com
называется она Intel® Math Kernel Library (Intel® MKL)

Но, как оказалось, эта библиотека перестала быть бесплатной
Можно подписаться на не коммерческую версию только для линукса


Другой вариант, посмотреть в сторону ассемблера, всякие ряды Тейлоры зашиты в FPU процессоров семейства x86 лет двадцать уже и сейчас выполняются за считанные такты. Спроси в ассемблерном разделе как реализовать нужные тебе функции

Я с ассемблером на ВЫ. Даже если ПОДСКАЖУТ, как реализовать, вряд ли это у меня получится . А есть ли готовые модули (dll), которые реализуют эти ASM алгоритмы?

Цитата

возьми интеловскую математическую библиотеку заточенную под математическое аппаратное обеспечение и получишь наивысшее сочетание точности и производительности

Бегло посморел MKL 6, тригонометрии там, кажется, нет.

Попробуйте дробно-рациональную аппроксимацию arcsin (интервал x, -1..1), устроит - приведу остальные

function arcsin(X:double):double;
Type ArrND=array[1..30] of double;
var I:integer;
P,Q:double;
A, B:ArrND;
begin
A[1]:= 1.95387311727746E-0005;
A[2]:=-1.52852221394629E+0000;
A[3]:=-2.01776736316873E-0003;
A[4]:= 7.76296426337929E-0001;
B[1]:=-1.52652894854093E+0000;
B[2]:=-1.74609738502257E-0003;
P:=A[4];
for i:=3 downto 1 do P:=x*P+A[I];
Q:=X+B[2];
for I:=1 downto 1 do Q:=x*Q+B[I];
Result:=P/Q;
end;

Есть такие формулы, возможно они тебе помогут.
они приводят arcsin и arccos к arctan

Прикреплённая картинка

и

Прикреплённая картинка

MeG Именно по этим формулам и считает ArcSin и ArcCos дельфи. А ArcTan эта одна команда сопроцессора - FPATAN.

Добавлено 12.01.06, 14:50
P.S Тема перекачивала отсюда: падение роизводительности

Цитата

P.S Тема перекачивала отсюда: падение роизводительности

Я бы не сказал, что тема перекочевала. Связь действительно есть, а вот проблемы обсудаем совершенно разные.

Кстати, P.O.D, Ты не представляешь себе, как только что помог мне фразой:

Цитата

MeG Именно по этим формулам и считает ArcSin и ArcCos дельфи. А ArcTan эта одна команда сопроцессора - FPATAN.

Я реализовал мою функцию "CalcAngle" через ArcTan, тем самым, судя по приведенным выше формулам, обойдя вычисление сразу 2-х корней и добившись заметного увеличения скорости.

Вот только если все можно было бы считать череза arctan...

это то ?

float fcos( float x ) {
  float4 values;
  float fv;
  values.x = x * 0.15915494;
  fv = frac( values.x );
  values.y = 0.25 - fv;
  values.z = fv - 0.75;
  values.x = max( values.y, values.z );
  values.x = -37.207532016 * values.x * (-0.168868639 + values.x * values.x);
  return values.x;
}
float fsin( float x ) {
  float4 values;
  float fv;
  values.x = x * 0.15915494 + 0.2499999951;
  fv = frac( values.x );
  values.y = 0.25 - fv;
  values.z = fv - 0.75;
  values.x = max( values.y, values.z );
  values.x = -37.207532016 * values.x * (-0.168868639 + values.x * values.x);
  return values.x;
}
float ftan( float x ) {
  // NOTE : complete formula tan(x) = (x - x*x*x/9 + x*x*x*x*x/945) / (1 - 4*x*x/9 + x*x*x*x/63);
  //             reduced formula  tan(x) = (x - x*x*x/9) / (1 - 4*x*x/9 + x*x*x*x/63);
  float4 values;
 
  x = x * 0.3183098861 + 0.5;
  values.x = frac( x ) * 3.141592653589 - 1.570796326794;
  values.y = values.x * values.x;
  values.z = values.y * values.x;
  values.w = values.z * values.x;
 
#ifdef FTAN_PERFECT
  float xp5 = values.w * values.x;
  return (values.x - values.z * 0.111111111111 + xp5 * 0.001058201058) / (1 - values.y * 0.44444444444 + values.w * 0.01587301587);
#else
  return (values.x - values.z * 0.111111111111) / (1.0 - values.y * 0.44444444444 + values.w * 0.01587301587);
#endif
}
float facos( float x ) {
  return 0.93 * (1.68902830838 + x * (x * x - 1.07526881720));
}
float fasin( float x ) {
  return 0.25f * x * (4.0 + x * x);
}
float ffmod( float x, float y ) {
  // NOTE : input range : x >= 0, y > 0
  return x - y * floor( x / y );
}
float fsqrt( float value ) {
        unsigned int tmp = (uint32(0x3f800000 << 1) + 0x3f800000 - *(uint32*)&x) >> 1;
  float y = *(float*)&tmp;
  return y * (1.47f - 0.47f * x * y * y);
}

Сообщение отредактировано: Sazabis - 12.01.06, 18:15

На вскидку сложно сказать. Пояснить немоного не мог бы? Это вообще, твои мысли или взял откуда?

Добавлено 12.01.06, 18:33
Вообще, круто - ни корней, ни циклов.

Добавлено 12.01.06, 18:36
Если б не экзамен завтра сдавать, опопробовал бы.

P. S. Вообще, заранее всем спасибо, за проявленный интерес к теме. Очень многое для себя открыл.

это формулы быстрого вычисления для чисел с плавающей точкой, накопленные в процессе разработки игр. взято с www.gamedev.ru или исходников Cg . Есть еще более быстрые вычисления для чисел с фиксированной точкой (самопальное представление числа с плавающей точкой в неком unsigned старшие биты отвечают за целую часть, младшие за дробную или наоборот, непомню ) Но для современных процессоров имеющих регистры под числа с плавающей точкой они не актуальны, да и код весь пришлось бы перелопачивать. Как конкретно работают эту формулы надо у математиков спрашивать, мне эти коэфициэнты мало о чем говорят.
Могу только добавить что вычисление синуса и косинуса может быть быстрее табличного, это связано с тем что на некоторых процах таблицы могут к кэш не залазить и конвеерное вычисление осуществляется быстрее.

Да, это объясняет, почему табличка на 1300 элементов, большая часть из которых (99%) не используется, не дает ожидаемых результатов.

В Абрамовиче есть неплохие аппроксимации для аркфункций - используются полиномы третьей и пятой степеней, достигается точность порядка 10^-5. Как я понимаю, в вашем случае более высокая точность и не нужна.

Формулы для арксинуса можно использовать без изменений, для арктангенса требуется предварительно привести аргумент к диапазону [-1, 1] (несложное преобразование).

Скриншот страницы с коэффициентами аппроксимации прилагается.

Прикреплённая картинка

Обещанный пример сравнения скорости выполнения п/п, вызываемых из .dll(вложение)

Насколько мне известно, в настоящее время наиболее быстрый код генерируют компиляторы C++ и Fortran от Intel при устоновке соответствующих ключей. Если писать на Delphi, то наиболее критичные участки надо писать на ассемблере, т.к. генерируемый код очень далек от оптимального. Примеры сравнений - в теме по перемножению матриц.

Сообщение отредактировано: PAV - 13.01.06, 08:54

Прикреплённый файл_Sample.zip (14.03 Кбайт, скачиваний: 110)

Цитата Sazabis @ 12.01.06, 19:04

это формулы быстрого вычисления для чисел с плавающей точкой, накопленные в процессе разработки игр. взято с www.gamedev.ru или исходников Cg . Есть еще более быстрые вычисления для чисел с фиксированной точкой (самопальное представление числа с плавающей точкой в неком unsigned старшие биты отвечают за целую часть, младшие за дробную или наоборот, непомню ) Но для современных процессоров имеющих регистры под числа с плавающей точкой они не актуальны, да и код весь пришлось бы перелопачивать.

Меня как раз интересует для фиксированной точки (поскольку на процессорах АРМ не то, что плавающей точки, даже деления нет)
"Есть еще более быстрые..." - это у тебя есть или вообще есть? Если это CORDIC, то скорее всего не пригодится (слишком много операций, хоть и очень простых)

Цитата Beeblebrox @ 26.07.07, 15:36

это у тебя есть или вообще есть

вообще есть, у меня в коде точно нету. Техника описана, например, у Андре Ламота.
почитай тут
http://www.gamedev.ru/code/forum/?id=50927