Вычисление синуса и косинуса

В предложенной программе синус и косинус вычисляются через тангенс половинного аргумента. Тангенс половинного аргумента вычисляется в интервале [0,pi/8] с помощью быстро сходящейся цепной дроби.

Тангенс вычисляется на интервале [0,pi/8], для этого используется цепная дробь:

tan(x) = x/(1-x²/(3-x²/(5-x²/(7-x²/...)))).

Для одинарной точности требуется 4 итерации, для двойной - 6-7.
Затем вычисляются синус и косинус на интервале [0,pi/4] через тангенс половинного аргумента: sin(x₁)=2*t/(1+t²), nbsp;   cos(x₁)=(1-t²)/(1+t²). Потом, применяя алгоритм к x₂=(pi/2-x₁), можно будет расширить интервал до [0,pi/2]. Применяя известные формулы тригонометрии, интервал расширяется до полного цикла: [-pi,pi].

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

 

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

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

/* Вычисление синуса и косинуса без мат. библиотеки. 
    Оптимизировано для чисел одинарной точности. 

   Вычисление синуса и косинуса на [0, PI/4]:
   void _Sico(float arg,float *sine,float *cosi);

   Вычисление синуса и косинуса на [0, PI/2]:
   void Sico(float arg,float *sine,float *cosi);

   Вычисление синуса и косинуса на [-PI, PI]:
   void Sico1p(float arg,float *sine,float *cosi);
*/

#define M_PI ((float)3.141592653589793)
#define M_PI4 (M_PI*0.25F)
#define M_PI2 (M_PI*0.5F)

/* Вычисление синуса и косинуса на 0-PI/4. MFRAC=4 для одинарной точности */
#define MFRAC 4	
void _Sico(float arg, float *sine, float *cosi) {
  int n,n2;
  float arg2,t;
  /* вычисление тангенса непрерывным делением */
  t=0.; arg*=0.5F; arg2=arg*arg; n=MFRAC-1; n2=(n<<1)+1;
  for(;n>=0;n--) {
    if(n>0) t=arg2/(n2-t);
    else t=arg/(1.F-t);
    n2--; n2--;
  }
  /* синус и косинус */
  arg=t*t; arg2=arg+1.F; arg2=1.F/arg2;
  *sine=t*arg2; *sine+=(*sine);
  *cosi=1.F-arg; *cosi*=arg2;
}

/* аргумент 0-PI/2 */
void Sico(float arg, float *sine, float *cosi) {
  if(arg<=M_PI4) _Sico(arg,sine,cosi);
  else _Sico(M_PI2-arg,cosi,sine);
}

/* первый период: -PI<=arg<=PI */
void Sico1p(float arg, float *sine, float *cosi) {
  int s=0;
  if(arg<0.F) {arg=-arg;s=1;}
  if(arg<=M_PI2) Sico(arg,sine,cosi);
  else {
    Sico(arg-M_PI2,cosi,sine);
    *cosi=-(*cosi);
  }
  if(s) *sine=-(*sine);
}