C++学习------cmath头文件的源码学习01

桑榆晚
875 阅读5分钟

携手创作，共同成长！这是我参与「掘金日新计划 · 8 月更文挑战」的第23天，点击查看活动详情

引言

cmath是C++对math.h头文件的封装，里面定义了一系列的数学函数，用来进行通用的数学计算和转换。我们来看看他的源码实现。

math.h

参考代码： www.aospxref.com/android-12.…

类型别名

这里将double和float类型定义了别名，方便后续的使用。

25  typedef double __double_t;
26  typedef __double_t double_t;
27  typedef float __float_t;
28  typedef __float_t float_t;

宏常量定义

通过builtin方法获取最值的情况，定义边界值，包括INFINITY和NAN

30  #define HUGE_VAL __builtin_huge_val()
31  #define HUGE_VALF __builtin_huge_valf()
32  #define HUGE_VALL __builtin_huge_vall()
33  
34  #define INFINITY __builtin_inff()
35  
36  #define NAN __builtin_nanf("")

通过直接定义的方式，定义一些浮点数运算的标志值，后续计算函数当中会使用到。

38  #define FP_INFINITE 0x01
39  #define FP_NAN 0x02
40  #define FP_NORMAL 0x04
41  #define FP_SUBNORMAL 0x08
42  #define FP_ZERO 0x10
43  
44  #if defined(__FP_FAST_FMA)
45  #define FP_FAST_FMA 1
46  #endif
47  #if defined(__FP_FAST_FMAF)
48  #define FP_FAST_FMAF 1
49  #endif
50  #if defined(__FP_FAST_FMAL)
51  #define FP_FAST_FMAL 1
52  #endif
53  
54  #define FP_ILOGB0 (-INT_MAX)
55  #define FP_ILOGBNAN INT_MAX
56  
57  #define MATH_ERRNO 1
58  #define MATH_ERREXCEPT 2
59  #define math_errhandling MATH_ERREXCEPT

POSIX扩展的常量，包括常用的自然常数，圆周率，对数等，便于后面的快速计算。

337  #define M_E		2.7182818284590452354	/* e */
338  #define M_LOG2E		1.4426950408889634074	/* log 2e */
339  #define M_LOG10E	0.43429448190325182765	/* log 10e */
340  #define M_LN2		0.69314718055994530942	/* log e2 */
341  #define M_LN10		2.30258509299404568402	/* log e10 */
342  #define M_PI		3.14159265358979323846	/* pi */
343  #define M_PI_2		1.57079632679489661923	/* pi/2 */
344  #define M_PI_4		0.78539816339744830962	/* pi/4 */
345  #define M_1_PI		0.31830988618379067154	/* 1/pi */
346  #define M_2_PI		0.63661977236758134308	/* 2/pi */
347  #define M_2_SQRTPI	1.12837916709551257390	/* 2/sqrt(pi) */
348  #define M_SQRT2		1.41421356237309504880	/* sqrt(2) */
349  #define M_SQRT1_2	0.70710678118654752440	/* 1/sqrt(2) */
350  
351  #define MAXFLOAT	((float)3.40282346638528860e+38)

同样的，GNU的扩展中也有类似的值

390  /* GNU extensions. */
391  
392  #if defined(__USE_GNU)
393  #define M_El            2.718281828459045235360287471352662498L /* e */
394  #define M_LOG2El        1.442695040888963407359924681001892137L /* log 2e */
395  #define M_LOG10El       0.434294481903251827651128918916605082L /* log 10e */
396  #define M_LN2l          0.693147180559945309417232121458176568L /* log e2 */
397  #define M_LN10l         2.302585092994045684017991454684364208L /* log e10 */
398  #define M_PIl           3.141592653589793238462643383279502884L /* pi */
399  #define M_PI_2l         1.570796326794896619231321691639751442L /* pi/2 */
400  #define M_PI_4l         0.785398163397448309615660845819875721L /* pi/4 */
401  #define M_1_PIl         0.318309886183790671537767526745028724L /* 1/pi */
402  #define M_2_PIl         0.636619772367581343075535053490057448L /* 2/pi */
403  #define M_2_SQRTPIl     1.128379167095512573896158903121545172L /* 2/sqrt(pi) */
404  #define M_SQRT2l        1.414213562373095048801688724209698079L /* sqrt(2) */
405  #define M_SQRT1_2l      0.707106781186547524400844362104849039L /* 1/sqrt(2) */
406  #endif

宏函数定义---分类函数

61  #define fpclassify(x) __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x)
62  
63  #define isfinite(x) __builtin_isfinite(x)
64  
65  #define isinf(x) __builtin_isinf(x)
66  
67  #define isnan(x) __builtin_isnan(x)
68  
69  #define isnormal(x) __builtin_isnormal(x)
70  
71  #define signbit(x) \
72      ((sizeof(x) == sizeof(float)) ? __builtin_signbitf(x) \
73      : (sizeof(x) == sizeof(double)) ? __builtin_signbit(x) \
74      : __builtin_signbitl(x))

fpclassify：返回输入数x的类型

  • FP_INFINITE：Positive or negative infinity (overflow)
  • FP_NAN：Not-A-Number
  • FP_ZERO：Value of zero
  • FP_SUBNORMAL：Sub-normal value (underflow)
  • FP_NORMAL：Normal value (none of the above) 这里我们很难找到通用的代码逻辑，因为这个操作是与具体的平台相关的，我们参考glibc库中的实现来看看这个函数的实现： (glibc的源代码在这里可以下载：www.gnu.org/software/li…)
 957 /* Return number of classification appropriate for X.  */
 958 # if ((__GNUC_PREREQ (4,4) && !defined __SUPPORT_SNAN__)              \
 959       || __glibc_clang_prereq (2,8))                          \
 960      && (!defined __OPTIMIZE_SIZE__ || defined __cplusplus)
 961      /* The check for __cplusplus allows the use of the builtin, even
 962     when optimization for size is on.  This is provided for
 963     libstdc++, only to let its configure test work when it is built
 964     with -Os.  No further use of this definition of fpclassify is
 965     expected in C++ mode, since libstdc++ provides its own version
 966     of fpclassify in cmath (which undefines fpclassify).  */
 967 #  define fpclassify(x) __builtin_fpclassify (FP_NAN, FP_INFINITE,        \
 968      FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x)
 969 # else
 970 #  define fpclassify(x) __MATH_TG ((x), __fpclassify, (x))
 971 # endif

同样的，我们来看看第二个实现，可以看到，这里使用宏连接的方式，判断参数的类型然后确定调用__fpclassifyf函数，还是__fpclassifyl函数（这种宏的拼接手法在glibc代码中很常见，同时这也导致了代码跳转非常多，变得晦涩难懂）

 922 # define __MATH_TG(TG_ARG, FUNC, ARGS)      \
 923   (sizeof (TG_ARG) == sizeof (float)        \
 924    ? FUNC ## f ARGS             \
 925    : sizeof (TG_ARG) == sizeof (double)     \
 926    ? FUNC ARGS                  \
 927    : FUNC ## l ARGS)
 928 #endif

继续往下追对应的函数__fpclassifyf，对应该函数，有两种实现，ieee754和riscv，这是两种不同的标准，我们优先看最为通用的ieee754

image.png

//glibc/sysdeps/ieee754/flt-32/s_fpclassifyf.c
 24 int
 25 __fpclassifyf (float x)
 26 {
 27   uint32_t wx;
 28   int retval = FP_NORMAL;
 29 
 30   GET_FLOAT_WORD (wx, x);
 31   wx &= 0x7fffffff;
 32   if (wx == 0)
 33     retval = FP_ZERO;
 34   else if (wx < 0x800000)
 35     retval = FP_SUBNORMAL;
 36   else if (wx >= 0x7f800000)
 37     retval = wx > 0x7f800000 ? FP_NAN : FP_INFINITE;
 38 
 39   return retval;
 40 }

这下一看，逻辑就比较清晰了，首先默认是FP_NORMAL，然后调用GET_FLOAT_WORD宏，作用如下： 定义了一个union结构体，用于将float数转换为uint32的字解释，一个非常精妙的方法。

//glibc/include/math.h
63 typedef union
64 {
65   float value;
66   uint32_t word;
67 } ieee_float_shape_type;
68 
69 /* Get a 32 bit int from a float.  */
70 #ifndef GET_FLOAT_WORD
71 # define GET_FLOAT_WORD(i,d)                    \
72 do {                                \
73   ieee_float_shape_type gf_u;                   \
74   gf_u.value = (d);                     \
75   (i) = gf_u.word;                      \
76 } while (0)
77 #endif

然后做了一个与操作，0x7fffffff，为什么是这个数呢？ 因为32位浮点数的表示中，第31位表示符号位，第30位到23位表示指数域，第22位到第0位表示小数域。 0x7fffffff截断了除符号位的其余数，即指数域+小数域；

  • 如果此时wx为0，那说明指数域+小数域都为0，也就返回FP_ZERO；
  • 如果wx < 0x800000（第23位为1），说明只有第0位到第22位有值，即小数域有值，但实际上指数域的范围是（1 – 254），不可能为0，所以返回FP_SUBNORMAL；
  • 如果wx == 0x7f800000，说明指数域占满了1,小数域全为0，这也是不可能的，标记为溢出FP_INFINITE
  • 如果wx > 0x7f800000，说明指数域全为1,小数域有数值，因为8位指数的取值范围为0255（0000000011111111），但在浮点数的阶码中，00000000与11111111被保留用作特殊情况，所以阶码可用范围只有1~254，总共有254个值，这样的数字是没法进行表示的，所以标识为FP_NAN（非数字）。 riscv的实现与其指令集相关，我们就不过多分析了。 __fpclassifyl函数的实现与__fpclassifyf类似，区别在于 64位浮点数1位符号位S，11位指数位E 和52位尾数位M

image.png

后续函数的分析见下一节

avatar
文章
阅读
粉丝