携手创作，共同成长！这是我参与「掘金日新计划 · 8 月更文挑战」的第27天，点击查看活动详情
续接上文：C++学习------cmath头文件的源码学习04
函数族定义---三角函数
- cos---计算cos值
- sin---计算sin值
- tan---计算tan值
- acos---计算arccos值
- asin---计算arcsin值
- atan---计算arctan值
- atan2---计算输入参数y/x的arctan值 因为三角函数和反三角函数的计算在底层做了很多计算上的简化模拟，充分使用了数值分析的原理，运用各种近似，查表等方法得到三角函数的近似值，我们来看其中一个具体的计算，其它函数的计算方式类似，可以参考对应的代码资料。 代码位于：glibc/sysdeps/ieee754/dbl-64/s_sin.c
193 /*******************************************************************/
194 /* An ultimate sin routine. Given an IEEE double machine number x */
195 /* it computes the rounded value of sin(x). */
196 /*******************************************************************/
197 #ifndef IN_SINCOS
198 double
199 SECTION
200 __sin (double x)
201 {
202 double t, a, da;
203 mynumber u;
204 int4 k, m, n;
205 double retval = 0;
206
207 SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
208
209 u.x = x;
210 m = u.i[HIGH_HALF];
211 k = 0x7fffffff & m; /* no sign */
212 if (k < 0x3e500000) /* if x->0 =>sin(x)=x */
213 {
214 math_check_force_underflow (x);
215 retval = x;
216 }
217 /*--------------------------- 2^-26<|x|< 0.855469---------------------- */
218 else if (k < 0x3feb6000)
219 {
220 /* Max ULP is 0.548. */
221 retval = do_sin (x, 0);
222 } /* else if (k < 0x3feb6000) */
223
224 /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/
225 else if (k < 0x400368fd)
226 {
227 t = hp0 - fabs (x);
228 /* Max ULP is 0.51. */
229 retval = copysign (do_cos (t, hp1), x);
230 } /* else if (k < 0x400368fd) */
231
232 /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
233 else if (k < 0x419921FB)
234 {
235 n = reduce_sincos (x, &a, &da);
236 retval = do_sincos (a, da, n);
237 } /* else if (k < 0x419921FB ) */
238
239 /* --------------------105414350 <|x| <2^1024------------------------------*/
240 else if (k < 0x7ff00000)
241 {
242 n = __branred (x, &a, &da);
243 retval = do_sincos (a, da, n);
244 }
245 /*--------------------- |x| > 2^1024 ----------------------------------*/
246 else
247 {
248 if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
249 __set_errno (EDOM);
250 retval = x / x;
251 }
252
253 return retval;
254 }
我们来看几个比较关键的过程：
- 数值转换 这一步是将输入的double x中的信息做分析
typedef int int4;
typedef union { int4 i[2]; double x; double d; } mynumber;
将64位的double信息转换为两个32为int，其中m是根据大小端序取的最高位的数据，然后截断符号位，赋值给k
209 u.x = x;
210 m = u.i[HIGH_HALF];
211 k = 0x7fffffff & m; /* no sign */
这里HIGH_HALF的定义也是值得我们观察的 glibc/include/endian.h，即如果大端序，则高位在前面表示，HIGH_HALF就表示int数组的第0位。
3 #if defined _LIBC && !defined _ISOMAC
4 # if __FLOAT_WORD_ORDER == __BIG_ENDIAN
5 # define BIG_ENDI 1
6 # undef LITTLE_ENDI
7 # define HIGH_HALF 0
8 # define LOW_HALF 1
9 # else
10 # if __FLOAT_WORD_ORDER == __LITTLE_ENDIAN
11 # undef BIG_ENDI
12 # define LITTLE_ENDI 1
13 # define HIGH_HALF 1
14 # define LOW_HALF 0
15 # endif
16 # endif
17 #endif
- 根据k的范围，做不同的处理 这里还是要强调一点，64位浮点数，第63位是符号位，62到52位是指数域，51到0位是小数域
（1）k < 0x3e500000，中间的指数域正好是0x3e5(997),这样计算出来的数就是1.M*2E(-26)，这种情况下，x已经足够小，x->0 =>sin(x)=x，所以我们可以使用x作为sin(x)的值返回，但是要做underflow检查;
（2）k < 0x3feb6000，指数域0x3f6(1014),即2^-26<|x|< 0.855469时，使用do_sin (x, 0)方式获取结果；
（3）k < 0x400368fd，即0.855469 <|x|<2.426265，使用如下的方式获取结果：
t = hp0 - fabs (x);
retval = copysign (do_cos (t, hp1), x);
（4）k < 0x419921FB，即2.426265<|x|< 105414350，使用如下的方式获取结果：
n = reduce_sincos (x, &a, &da);
retval = do_sincos (a, da, n);
（5）k < 0x7ff00000,指数域0x7ff(2047),即105414350 <|x| <2^1024，使用如下的方式获取结果：
n = __branred (x, &a, &da);
retval = do_sincos (a, da, n);
（6）|x| > 2^1024，先考虑输入数溢出的情况(即指数域=2047且小数域全为0)
if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
__set_errno (EDOM);
否则返回x/x即1。 中间的每一步的具体计算都是使用了具体的数值分析知识做的近似估计计算，后续另开文章讲解。
