Оптимизация sin() на BASM

[InternetStranger]

Уважаемые, Кодеры! Помогите пожалуйста ускорить функцию вычисления sin(x) через полином в форме Горнера (Delphi).

История проблемы:
Для вычисления sin(x) при x=[-Pi/2,Pi/2] вполне достаточно учесть 3 члена ряда. Сия замечательная функция реализуется вот так:

function ASin(fAngle: Single): Single;
const
  sin1: Double = 7.61e-03;
  sin2: Double = -1.6605e-01;
  sin3: Double = 1;
asm
    { Result := fAngle * (sin3 + sqr(fAngle) * (sin2 + sin1 * sqr(fAngle) )) }
    fld fAngle.Single 
    fmul fAngle.Single 
    fld sin1.Double    
    fmul st(0),st(1)    
    fld sin2.Double    
    faddp //st(1),st(0)
    fmulp //st(1),st(0) 
    fadd sin3.Double    
    fmul fAngle       
end;

Однако в моей задаче аргумент x меняется в пределах [-2,8], поэтому приходится учитывать гораздо больше членов ряда:

function mysin_11(x: real): real;
var
  x2: real;
const
  okr = Pi / 4; // Окрестность точки Pi/4
  A3 = -2.66592780638819;
  A5 = 2.12775693434532;
  A7 = -0.79794217405921;
  A9 = 0.16264977471553;
  A11 = -0.01573479830529;
begin
  x := 0.25 * x;
  x2 := sqr(x - okr);
  Result := -4 * (x - okr) * (1 + x2 * (A3 + x2 * (A5 + x2 * (A7 + x2 * (A9 + x2* A11)))));
end;

Ежу понятно, что глупо столько раз перемножать, но это эффективнее стандартной реализации sin() - ходят слухи, что она съедает > 200 тактов. (применительно к моей задаче ускорение c mysin_11() ~ 15%)
Вроде решение очевидно -> использовать симметрию функции, но нельзя использовать условные конструкции (в перспективе алгоритм будет распараллеливаться -> ветвления губительны)

Собственно задача:
Существует отличная формула приведения, позволяющая свести вычисление синуса произвольного угла к [-Pi/2, Pi/2] :

sin( n*Pi+y ) = (-1)^n * sin( y )

Нужно дополнить функцию Asin(x). Ломаю голову как эффективно выделить n и y и возвести (-1)^n и затратить при этом меньше тактов, чем для mysin_11(x).
Вообще есть способ избежать операции деления?

Зы: Везде двойная точность

[Mikl___]

InternetStranger
не очень понял, но попробую помочь
1) как эффективно возвести (-1)^n и затратить при этом меньше тактов
проверить последний бит у n если 0, то (-1)^n=1 если 1, то (-1)^n=-1
2) Вообще есть способ избежать операции деления?
умножение занимает меньше времени, чем деление, поэтому, если делят на константу х, то предварительно находят обратную ей 1/х и заменяют деление на умножение
3) а зачем искать синус рядами используя команды fpu? если в fpu есть команды fsin и fsincos, а для вычисления arcsin(x) используют fpu-шную команду fpatan
arcsin(x)=arctg(tg(x))=arctg(sin(x)/sqrt(1-sin(x)^2))

Однако в моей задаче аргумент x меняется в пределах [-2,8]

для начала sin(-x)=-sin(x)
x=n*Pi + y
y = x mod Pi (остаток от деления х на Pi)
n = x div Pi (целочисленное деление х на Pi)
последний бит(n) = n and 1 или последний бит(n) = n mod 2

[InternetStranger]

Спасибо за советы ) Половину проблемы ты мне уже разрешил! ))

1) Проверяя последний бит у n, я так понимаю, мы проверяем на четность? ( ааа. сообразил -> это проверка кратности двум))
Даже если это так, интересно конечно, но мне нельзя использовать условные переходы -> иначе не смогу распараллеливать потом алгоритм. Хотя...
С учетом твоего совета в голову вот что пришло:

(-1)^n = 1 + (n and 1)*(-2) = 1 - (n and 1) shr 1

Всего одна вещественных операция + по мелочи (shr - выполняется вроде за один такт, and наверное тоже).
Вот про это я и говорил, спасибо )))

2) Действительно, вообще сглупил. Pi - это же константа. Пора завязывать ночами напролет кодить
Но этот вопрос снимается уже - я просто нормировал (сделал замену переменной) все уравнение вместе с входящим в него синусом. Теперь осталось реализовать собственно операцию:
n= trunc (x), y = frac(x)

(Чтобы воспользоваться формулой приведения, предварительно нужно вычислить n и y ).

Или я плохо знаю матчасть, но mod и div в качестве операндов принимают Целые числа , а не Double? Я как-то всю жизнь считал не применимыми к вещественной арифметике.

3) А производительность в тактах этой команды известна? С тем учетом, что считает она до 15го знака после запятой и на всей области -> мне столько не нужно ))
Ряд до определенной длины сосчитать быстрее, чем вызвать эту команду (Ну как-то так я посчитал). Это если при вычислении ряда воспользуешься самыми базовыми операциями - * + - shl shr.

Итак, с учетом вышесказанного, что же мне еще осталось сделать (Для наглядности полу-оптимизированный Delphi - эквивалент ):

{Аргумент fAngle уже попадает сюда поделенным на Pi}

function ASin(fAngle: Double ): Double ;
const
  sin1: Double = 7.61e-03;
  sin2: Double = -1.6605e-01;
  sin3: Double = 1;
var n: Byte;
     y: Double;
begin
    n:= trunc(fAngle); // Вот это нужно сделать.
    y:= fAngle - n; // Наверняка это самый быстрый способ

    Result := fAngle * (sin3 + sqr(fAngle) * (sin2 + sin1 * sqr(fAngle) )); // Это уже реализовано
    Result := ( 1 - (n and 1) shr 1 ) * Result; // За это спасибо Mikl___ 
end;

ЗЫ: Эта функция лишь используется лишь для тестирования. Потом ее код будет тупо использоваться посреди уравнения, это же ряд, его хоть куда воткнуть можно. Поэтому это быстрее, чем вызывать всякие функции.

Т.е. теперь проблема формулируется так -> как числа двойной точности достать целую и дробные части?
Так-то я пробовал fist, но резобрался в какую сторону она округляет.

[Mikl___]

InternetStranger
Если требуется скорость для вычисления sin, то для чего тогда использовать fpu? Для больших скоростей используют целочисленную арифметику, а при вычислении тригонометрических функций используют табличное преобразование -- величина программы я думаю тебя не интересует, а значения угла как правило с определенным шагом получают с какого-нибудь датчика, погрешность измерения это и определяет точность вычисления

n:= trunc(fAngle); // Вот это нужно сделать.
y:= fAngle - n;

можно сделать так
fld Angle; от -2 до 8
fldpi; загрузили pi= 3,141593265358979323…
fdiv st(1),st(0)
fld st(0); дублируем результат
frndint; округляем до целого
fist n; сохраняем целую часть
fsub st(1),st(0); получаем остаток
Но всеравно, я считаю, что встроенная fsin считает быстрее, чем вычислять синус рядами, и это необходимо тебе проверить

Так-то я пробовал fist, но резобрался в какую сторону она округляет.

Команда FRNDINT округления вещественного числа до целого = ”Floating point ROUND to INTEGER”
Синтаксис команды: FRNDINT
Семантика команды: Команда округления до целого
Псевдокод команды: ST(0)<--RoundToIntegerValue(ST(0))
Эта команда округляет текущее содержимое вершины стека до целого числа. Биты RC (Rounding Control) в регистре управления (RCW) определяют способ округления результатов команд блока FPU до заданной точности
RC |Способ округления
00b| к ближайшему числу near
01b| к -∞ down
10b| к +∞ up
11b| к нулю chop/zero

;пример управления округлением
fstcw ax
and eax,1111xx1111111111b;11-10 bits rc
or eax,xx0000000000b
push eax
fldcw word ptr [esp]
pop eax
frndint

Команда инициализации FINIT/FNINIT (FPU INITIALIZE) сбрасывает FPU в состояние, в которое он попадает при включении питания. Рекомендуется использовать эту команду перед обращением к FPU. В результате FINIT:
rc=0 (режим округления к ближайшему числу);
флаги особых ситуаций в регистре SWR сброшены;
флаги запрета прерываний в регистре CWR по особым ситуациям установлены;
TOP=7.

[InternetStranger]

Спасибо за подробное разъяснение этой инструкции.
Буду разбираться )) Осталось немного до полноценной победы над синусом (а борюсь я с ним уже больше месяца)).

1) Целочисленная математика для меня не вариант, я уже не раз пытался свести мое уравнение к целым числам. все безуспешно.

2) Таблицы - тоже пройденный этап. У меня ж задача научная - все подсчитано. Для моей задачи (вернее необходимой точности) потребуется таблица ~ 4 мегабайта .
К тому же просто навскидку подумай, думаешь получится организовать быстрый доступ к значениям этой таблицы, когда у меня аргументы дробные. Это опять округления, отбрасывания вещественной части части и прочее....
По поводу таблиц ты конечно прав отчасти. Они используются:
2.1) там, где не требуется высокая точность. Например, при рисовании графики.
2.2) там, где требуется очень большая точность. тут таблицы - это отправные точки для итерационных методов ( вот эти алгоритмы кстати эффективнее полиномиального разложения для точности > 1e-10 )


Понимаю, выглядит на борьбу с тенью...
Вообще вопрос не столько программерский - больше философский: "Чтобы решить проблему, нужно прежде осознать, что она существует". А вычисление синуса - это действительно проблема.
Ведь большинство считает, что алгоритмы крупных компаний (aka AMD, Intel, nVidia...) наилучшие )) Ничего гениального там нет- пользуются они математикой, придуманной еще 17-19 веках. Ряд Тейлора/Маклорена вроде даже раньше ))
Я даже на 99% знаю какая аппаратная реализация используется в сопроцессорах ( см. пункт 2.2 ). А также ее сильные и слабые стороны. И вот за этой инструкцией "fsin" кроется куча не нужных и кое-каких бесполезных действий. Это всмысле для точности 1e-5 (которая нужна мне).

[InternetStranger]

Ну довольно лирики...

Но всеравно, я считаю, что встроенная fsin считает быстрее, чем вычислять синус рядами, и это необходимо тебе проверить

Вот результаты тестирования (дабы разубедить тебя во всесильности FPU). Знаю, что способ кривоватый, но этого вполне достаточно, чтобы продемонстрировать превосходства полиномиального разложения.

Итак, стоит задача рассчитать синус с точностью до 1e-5 в пределах [0..Pi]:

 
{ Пустая процедура }
function B_empty(x: Real): Real;   // Пустая процедура
const A0:Real= 1;    A2:Real= -0.49993908064694;    A4:Real=  0.04151144082871;      A6:Real=  -0.00127685584279; Pi_popolam:Real= -1.5707963267948966192;
asm  fld [x]
end;

{ Реализовано на Delphi }
function BSin(x: Real): Real;// [0,Pi], точность более 1e-5 
const A0:Real= 1;    A2:Real= -0.49993908064694;    A4:Real=  0.04151144082871;      A6:Real=  -0.00127685584279; Pi_popolam:Real= -1.5707963267948966192;
var x2 : real;
begin
     x2:= sqr(x+Pi_popolam); // компилятор эту строчку замечательно преобразовамал - все как надо
     Result := A0 + x2*( A2 + x2*( A4 + A6*x2 ) );
end;

{ Bsin - asm1 }
function BSin_asm1(x: Real): Real;//  (Моя реализация)
const A0:Real= 1;    A2:Real= -0.49993908064694;    A4:Real=  0.04151144082871;      A6:Real=  -0.00127685584279; Pi_popolam:Real= -1.5707963267948966192;
asm
   {x2:= sqr(x+Pi_popolam); }
   fld [x.Real]
   fadd [Pi_popolam.Real]
   fmul st(0),st(0)
   {Result := A0 + x2*( A2 + x2*( A4 + A6*x2 ) );     }
   fld [A6.Real]
   fmul st(0),st(1)
   fadd [A4.Real]
   fmul st(0),st(1)
   fadd [A2.Real]
   fmul st(0),st(1)
   fadd [A0.Real]
   ffree st(1)
end;

{ Bsin - asm2 }
function BSin_asm2(x: Real): Real;//   (Оптимизация от Борланда)
const A0:Real= 1;    A2:Real= -0.49993908064694;    A4:Real=  0.04151144082871;      A6:Real=  -0.00127685584279; Pi_popolam:Real= -1.5707963267948966192;
var x2 : real;
asm
   {x2:= sqr(x+Pi_popolam); }
   fld [x]
   fadd [Pi_popolam]
   fmul st(0),st(0)
   fstp  [x2]
   {Result := A0 + x2*( A2 + x2*( A4 + A6*x2 ) );     }
   fld  [A6]
   fmul [x2]
   fadd [A4]
   fmul [x2]
   fadd [A2]
   fmul [x2]
   fadd [A0]
end;

{ Bsin - asm3 - FPU_fsin }
function BSin_FPU_sin(x: Real): Real; ( Синус через команду FPU )
asm
   fld [x]
   fsin
end;

А вот и результаты тестов на моем AMD Athlon 7750 Dual-Core / 2.89GHz / 2Gb Hynix PC-6400 :

Empty = 31.81
----------------------------
Bsin - BorlandRTL (Standart) = 127.90
Bsin - BorlandRTL (Standart)(abs) = 96.08
----------------------------
Bsin - FPU_fsin = 133.50
Bsin - FPU_fsin (abs) = 101.69
----------------------------
Bsin - Delphi = 94.51
Bsin - Delphi (abs) = 62.69
----------------------------
Bsin - asm1 = 59.67
Bsin - asm1 (abs) = 27.85
----------------------------
Bsin - asm2 = 78.20
Bsin - asm2 (abs) = 46.39
----------------------------

(abs) - это число тактов, без учета вызова самой процедуры (для этого и замеряется время вызова процедуры Empty)

Прошу громко не пинать ( ассемблером владею всего пару дней, с даты первого поста этой темы ), а при необходимости указать где чего поменять для объективности результатов тестирования.

Эти значения имеют смысл только при сравнении между собой. Но 27.85 - это, кстати, тоже многовато )) Существуют и другие способы подсчета суммы полинома. Но это уже совсем другая история...

ЗЫ: в принципе необязательно было даже ухитряться и измерять. Достаточно просто заглянуть вот сюда: Intel®64 and IA-32 Architectures Optimization Reference Manual (ноябрь 2009), чтобы раз и навсегда отбило желание использовать fsin ))

Наиболее интересны страницы 569 и 570.