Алгоритм вычисления арктангенса

ОГЛАВЛЕНИЕ

Арктангенс вычисляется методом сужения области определения до [0,pi/12] и аппроксимации в этой области. Алгоритм оптимизирован для чисел с плавающей точкой одинарной точности. Арксинус легко вычисляется через арктангенс (с дополнительным использованием вычисления квадратного корня), а арккосинус - через арксинус.

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

  1. Вначале проверить знак x, изменить знак, сделав аргумент неотрицательным.
  2. Если x>1, обратить его: x1=1/x.
  3. Сокращаем область определения, используя формулу: atan(x)=pi/6+atan((x*sqrt(3)-1)/(x+sqrt(3))). Здесь sqrt(3)3. При этом необходимо запомнить число шагов (возможно, ноль).
  4. Арктангенс на интервале [0,pi/12] аппроксимируется формулой (для одинарной точности, в случае двойной точности формула должна быть улучшена!): atan(x) = x*(0.55913709/(1.4087812+x2) +0.60310579-0.05160454*x2)
  5. К полученному результату добавляется столько pi/6, сколько было шагов сокращения области определения.
  6. В случае обращения, аргумента, результат вычитается из pi/2.
  7. Если была смена знака, у результата меняем знак.

Для повышенной точности, формулу на участке [0,pi/12] следует брать в виде:
atan(x) = x*(m0+n0*(x*x)+k0/(m1+n1*(x*x)+k1/(m2+n2*(x*x)+k2/(...)))),
то есть в том же виде цепной дроби, как и для single precision, только с некоторыми другими значениями m0,n0,k0;m1,n1,k1;... Определению подобных значений будет посвящена задача, представляющая из себя частный случай алгоритма минимизации функции нескольких переменных.

Для вычисления арксинуса следует воспользоваться соотношением:
asin(x)=atan(x/sqrt(1-x2)),   где sqrt(x) - квадратный корень.
Арккосинус связан с арксинусом соотношением acos(x)=pi/2-asin(x). Для его вычисления следует сначала воспользоваться вычислением арксинуса. Для обоих случаев необходимо отслеживать случаи x=1, x=-1 (где требуется находить арктангенс бесконечности), а также случаи выхода x из области определения, где следует генерировать ошибку. В программе arcsin.c, приведенной ниже, используются литералы NORMAL для указания нормальной работы и EDOM для завершения с ошибкой.