引言
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 0x0139 #define FP_NAN 0x0240 #define FP_NORMAL 0x0441 #define FP_SUBNORMAL 0x0842 #define FP_ZERO 0x1043 44 #if defined(__FP_FAST_FMA)45 #define FP_FAST_FMA 146 #endif47 #if defined(__FP_FAST_FMAF)48 #define FP_FAST_FMAF 149 #endif50 #if defined(__FP_FAST_FMAL)51 #define FP_FAST_FMAL 152 #endif53 54 #define FP_ILOGB0 (-INT_MAX)55 #define FP_ILOGBNAN INT_MAX56 57 #define MATH_ERRNO 158 #define MATH_ERREXCEPT 259 #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
//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.h63 typedef union64 {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_WORD71 # 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
// glibc/sysdeps/ieee754/ldbl-128/s_fpclassifyl.c24 int 25 __fpclassifyl (_Float128 x) 26 { 27 uint64_t hx, lx; 28 int retval = FP_NORMAL; 29 30 GET_LDOUBLE_WORDS64 (hx, lx, x); 31 lx |= (hx & 0x0000ffffffffffffLL); 32 hx &= 0x7fff000000000000LL; 33 if ((hx | lx) == 0) 34 retval = FP_ZERO; 35 else if (hx == 0) 36 retval = FP_SUBNORMAL; 37 else if (hx == 0x7fff000000000000LL) 38 retval = lx != 0 ? FP_NAN : FP_INFINITE; 39 40 return retval; 41 }
复制代码
评论