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

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

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

1. Вначале проверить знак x, изменить знак, сделав аргумент неотрицательным.
2. Если x>1, обратить его: x₁=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+x²) +0.60310579-0.05160454*x²)
   
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-x²)),   где sqrt(x) - квадратный корень.
Арккосинус связан с арксинусом соотношением acos(x)=pi/2-asin(x). Для его вычисления следует сначала воспользоваться вычислением арксинуса. Для обоих случаев необходимо отслеживать случаи x=1, x=-1 (где требуется находить арктангенс бесконечности), а также случаи выхода x из области определения, где следует генерировать ошибку. В программе arcsin.c, приведенной ниже, используются литералы NORMAL для указания нормальной работы и EDOM для завершения с ошибкой.

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

Ниже приведена программа вычисления арктангенса по указанному алгоритму. Программа вычисления арксинуса и арккосинуса, приведенная еще ниже, использует как внешние ее и программу для квадратного корня.

/* вычисление арктангенса без мат. библиотеки */

   float Arctan(float x);
*/

#define M_PI ((float)3.141592653589793)
#define M_PI12 (M_PI/12.F)
#define M_PI6 (M_PI/6.F)
#define M_PI2 (M_PI/2.F)
/* квадратный корень из 3 */
#define SQRT3 ((float)1.732050807569)

float Arctan(float x) {
  int sta=0,sp=0;
  float x2,a;
  /* проверка смены знака */
  if(x<0.F) {x=-x;sta|=1;}
  /* проверка инверсии */
  if(x>1.F) {x=1.F/x;sta|=2;}
  /* сжатие области, пока не x<PI/12 */
  while(x>M_PI12) {
    sp++; a=x+SQRT3; a=1.F/a; x*=SQRT3; x-=1.F; x*=a;
  }
  /* вычисления */
  x2=x*x; a=x2+1.4087812F; a=0.55913709F/a; a+=0.60310579F;
  a-=0.05160454F*x2; a*=x;
  /* пока не sp=0 */
  while(sp>0) {a+=M_PI6;sp--;}
  if(sta&2) a=M_PI2-a;
  if(sta&1) a=-a;
  return(a);
}

/* Вычисление арксинуса и арккосинуса. Использует Arctan() и Sqroot() как внешние вызовы. 

   int Arcsin(float *a,float x);
       x - аргумент,
       a - указатель на результат
   Возвращает код ошибки NORMAL или EDOM

   int Arccos(float *a,float x);
       x - аргумент,
       a - указатель на результат
   Возвращает код ошибки NORMAL или EDOM

*/

#define M_PI ((float)3.141592653589793)
#define M_PI2 (M_PI/2.F)

enum{NORMAL,EDOM};

extern float Arctan(float x);
extern float Sqroot(float x);

int Arcsin(float *a,float x) {
 /* проверка исключений */
 if(x<-1.F) {*a=-M_PI2;return(EDOM);}
 if(x==-1.F) {*a=-M_PI2;return(NORMAL);}
 if(x>1.F) {*a=M_PI2;return(EDOM);}
 if(x==1.F) {*a=M_PI2;return(NORMAL);}
 /* преобразование аргумента */
 x/=Sqroot(1.F-x*x);
 *a=Arctan(x);
// положительный результат
 return(NORMAL);
}

int Arccos(float *a,float x) {
 int re=Arcsin(a,x);
 *a=M_PI2-(*a);
 return(re);
}