|
Copyright © Nikitine Valeri F. 2000
Синус и косинус вычисляются через тангенс половинного аргумента. Последний
вычисляется в интервале [0,pi/8] с помощью быстро сходящейся цепной дроби.
На интервале [0,pi/8] вычисляется тангенс, для чего используется
цепная дробь:
tan(x) =
x/(1-x2/(3-x2/(5-x2/(7-x2/...)))).
Для single precision требуется 4 итерации, для double -
6-7.
Затем синус и косинус на интервале [0,pi/4] вычисляются
через тангенс половинного аргумента:
sin(x1)=2*t/(1+t2), nbsp;
cos(x1)=(1-t2)/(1+t2). Затем,
применяя алгоритм к x2=(pi/2-x1), можно расширить
интервал до [0,pi/2]. Применяя известные формулы тригонометрии, интервал
расширяется до полного цикла: [-pi,pi]. Дальнейшее расширение достаточно
тривиально, но здесь не сделано -- сами пробуйте если хотите.
Имеется альтернативный подход, требующий большего числа операций, но удобный
при представлении аргумента в виде fixed point. В этом случае алгоритм связан с
проведением серии комплексных умножений на табличные данные. Будет приведен в
дальнейшем -- ждите добавления.
Ниже приведена программа, реализующая вычисление синуса и косинуса через
тангенс.
/* Sine and cosine without mathematic library. Optimized for floating
point single precision.
Copyright (c) Nikitin V.F. 2000
Calculate sine and cosine within [0, PI/4]:
void _Sico(float arg,float *sine,float *cosi);
Calculate sine and cosine within [0, PI/2]:
void Sico(float arg,float *sine,float *cosi);
Calculate sine and cosine within one period [-PI, PI]:
void Sico1p(float arg,float *sine,float *cosi);
No argument domain check is performed: insert yourself.
*/
#define M_PI ((float)3.141592653589793)
#define M_PI4 (M_PI*0.25F)
#define M_PI2 (M_PI*0.5F)
/* sine and cosine within 0-PI/4. MFRAC=4 optimized for single precision */
#define MFRAC 4
void _Sico(float arg, float *sine, float *cosi) {
int n,n2;
float arg2,t;
/* calculate tangent by continuous fraction */
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--;
}
/* sine and cosine */
arg=t*t; arg2=arg+1.F; arg2=1.F/arg2;
*sine=t*arg2; *sine+=(*sine);
*cosi=1.F-arg; *cosi*=arg2;
}
/* argument 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);
}
/* first period: -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);
}
|
