在编程中我们总要进行一些数学运算以及数字处理,尤其是浮点数的运算和处理,这篇文章主要介绍C语言下的数学库。而其他语言中的数学库函数的定义以及最终实现也是通过对C数学库的调用来完成的,其内容大同小异,因此就不在这里介绍了。
C语言标准库中的math.h定义了非常多的数学运算和数字处理函数。这些函数大部分都是在C89标准中定义的,而有些C99标准下的函数我会特殊的说明,同时因为不同的编译器下的C标准库中有些函数的定义有差别,我也会分别的说明。
数字的范围
整型
整型用来存储整数数值,它按存储的字节长短分为:字符型,短整型,整型,长整型。 所有类型的存储长度都是定长的。既然类型是定长的就有一个最大最小可表示的范围,对于整型来说各种类型的最大最小的定义可以在limits.h中找到。下面表格列出了不同类型的存储长度和最大最小值:
类型 | 字节数 | 最小值 | 宏定义 | 最大值 | 宏定义 | 备注 |
---|---|---|---|---|---|---|
char | 1 | -2^7 | SCHAR_MIN | 2^7-1 | SCHAR_MAX | |
unsigned char | 1 | 0 | UCHAR_MIN | 2^8-1 | UCHAR_MAX | |
short | 2 | -2^15 | SHRT_MIN | 2^15-1 | SHRT_MAX | |
unsigned short | 2 | 0 | USHRT_MIN | 2^16-1 | USHRT_MAX | |
int | 4? | -2^31 | INT_MIN | 2^31-1 | INT_MAX | |
unsinged int | 4? | 0 | UINT_MIN | 2^32-1 | UINT_MAX | |
long | 4? | -2^31 | LONG_MIN | 2^31-1 | LONG_MAX | |
unsigned long | 4? | 0 | ULONG_MIN | 2^32-1 | ULONG_MAX | |
long long | 8 | -2^63 | LLONG_MIN | 2^63-1 | LLONG_MAX | C99 |
unsigned long long | 8 | 0 | ULLONG_MIN | 2^64-1 | ULLONG_MAX | C99 |
对于int和long类型来说,二者的长度是依赖于操作系统的字长或者机器的字长。因此如果我们要编写跨平台或跨系统的程序就应该尽量减少对这两个类型变量的直接定义。 下面表格列出了int和long两种类型在不同操作系统字长下的长度。
类型 | 16位系统/字节 | 32位系统/字节 | 64位系统/字节 |
---|---|---|---|
int | 2 | 4 | 4 |
long | 4 | 4 | 8 |
在很多系统中都对32位的整型以及64位的整型进行特殊的定义,比如Windows中的DWORD,UINT32,INT64等等。
浮点型
浮点型用来存储浮点数值。它按精度分为:单精度浮点型,双精度浮点型,扩展双精度浮点型。 浮点数是连续并且无限的,但是计算机并不能表达出所有连续的值。因此对浮点数定义了最小规格化值和最大规格化值,这些定义可以在float.h中找到。下面表格列出了不同类型的存储长度和最值:
类型 | 字节数 | 最小规格化值 | 宏定义 | 最大规格化值 | 宏定义 | 备注 |
---|---|---|---|---|---|---|
float | 4 | 1.175494351e-38 | FLT_MIN | 3.402823466e+38 | FLT_MAX | |
double | 8 | 2.2250738585072014e-308 | DBL_MIN | 1.7976931348623158e+308 | DBL_MAX | |
long double | 8? | 2.2250738585072014e-308 | LDBL_MIN | 1.7976931348623158e+308 | LDBL_MAX | C99 |
这里的FLT_MIN,DBL_MIN,LDBL_MIN并不是指最小可表示的浮点数,而是最小规格化浮点值,具体我会在下面详细介绍。
对 long double 的定义,取决于编译器和机器字长,所以对于不同平台可能有不同的实现,有的是8字节,有的是10字节,有的是12字节或16字节。
为了和数学中的无穷∞对应,标准库中定义了一个宏:INFINITY来表示无穷大。比如1.0/0.0等于INFINITY,-1.0/0.0等于-INFINITY。无穷大可以进行加减乘除操作,比如1.0/INFINITY == 0。
为了和数学中的非法数字对应,标准库中定义了一个宏:NAN来表示非法数字。比如负数开方、负数求对数、0.0/0.0、0.0* INFINITY、INFINITY/INFINITY、INFINITY-INFINITY这些操作都会得到NAN。注意:如果是整数0/0会产生操作异常
浮点数的存储结构
浮点数不像整数那样离散值,而是连续的值。但是用计算机来描述一个浮点数时就不可能完全实现其精度和连续性,现在的浮点型的存储和描述普遍都是遵循IEEE754标准。如果您想详细的了解关于浮点数的存储格式那么您可以花费一点时间来阅读:https://wenku.baidu.com/view/d02978d8d15abe23482f4dac.html 这篇文章。
简单来说浮点数的存储由:S(sign)符号位、E(exponent)指数位、M(mantissa 或significand)尾数位三个部分组成。我们以一个32位的float类型举例来说,一个浮点数N的从高位到低位的存储结构如下:
浮点数的存储结构
也就是一个32位的浮点数由1个符号位,8个指数位,23个尾数位组成。 而为了表示不同类型的浮点数,根据存储格式对浮点数进行了如下分类:
如果一个浮点数中指数位部分全为1,而尾数位部分全为0则这个浮点数表示为无穷大** INFINITY **,如果符号位为0表示正无穷大,否则就是负无穷大。
如果一个浮点数中指数位部分全为1,而尾数位部分不全为0则这个浮点数表示为非法数字NAN。因此可以看出非法数字并非一个数字而是一类数字。在下面介绍nan函数时我会更加深入的介绍NAN
如果一个浮点数中除符号位外全部都是0,那么这个浮点数就是0
如果一个浮点数中指数位部分全为0,而尾数位部分不全为0则这个浮点数称为非规格化浮点数,英文称为:subnormal number 或 denormal number 或 denormalized number。非规格化浮点数常用来表示一个非常接近于0的浮点数。
如果一个浮点数中的指数位部分即非全1又非全0。那么这个浮点数称之为规格化浮点数,英文称之为:normal number。我们上面定义的FLT_MIN, DBL_MIN 指的就是最小的规格化浮点数。
我们把规格化浮点数和非规格化浮点数合称为可表示的浮点数,英文称之为:machine representable number
一个规格化浮点数N的值可以用如下公式算出:
规格化浮点数计算公式
从上面的公式中可以看出对于一个32位浮点数来说,指数位占8位,最小值是1(全0为非常规浮点),而最大值是254(全1为无穷或者非法浮点),而减去127则表示指数部分的最小值为-126,最大值为127;同时我们发现除了23位尾数外,还有一个隐藏的1作为尾数的头部。因此我们就很容易得出:
**FLT_MIN = 1.0 * 2^-126 = 1.175494351e-38 **
FLT_MAX = (1.11111111111111111111111)b * 2^127 = 3.402823466e+38
一个非规格化浮点数N的值的可以用如下公式算出:
非规格化浮点数计算公式
从上面的公式中可以看出对于一个32位的浮点数来说,我们发现虽然非规格化浮点的指数位部分全0,但是这里并不是0-127,而是1-127,同时发现尾数位部分并没有使用隐藏的1作为尾数的头部,而是将头部的1移到了指数部分,这样做的目的是为了保持浮点数字的连续性。我们可以看出当一个浮点数小于FLT_MIN时,他就变为了一个非规格化浮点。我们知道FLT_MIN的值是1.0 * 2^-126,而一个比FLT_MIN小的值就应该是:(0.11111111111111111111111)b * 2^-126,而一个比0大的值就是:(0.00000000000000000000001)b * 2^-126。如果非规格化浮点数以-127作为指数,而继续使用1作为尾数的头部时,那么这种数字连续性将会被打破。这也是为什么要定义规格化浮点数和非规格化浮点数的意义所在。可以看出浮点数的这种存储设计的精妙之处!!。
从上面两种类型的浮点数中可以总结出浮点数的计算公式可以表示为:
** N = 符号 * 尾数 * 2^指数 **
数学函数
数字判断函数或宏
//如果x是正无穷大返回1,负无穷大返回-1,否则返回0int isinf(x)//如果x是无穷大返回0int isfinite(x)//如果x是一个规格化浮点数则返回非0int isnormal(x)//如果x是一个非法的数字返回非0int isnan(x)//如果x是负数返回非0int signbit(x) /** *返回浮点数的分类: FP_INFINITE: x是无穷大或者无穷小 FP_NAN:x是一个非法数字 FP_NORMAL:x是一个规格化浮点数 FP_SUBNORMAL:x是一个非规格化浮点数 FP_ZERO:x是0 */int fpclassify(x)
三角函数
1. 反余弦函数: y = arccos(x)
extern float acosf(float x);extern double acos(double x);extern long double acosl(long double x);
2. 反正弦函数:y = arcsin(x)
extern float asinf(float x);extern double asin(double x);extern long double asinl(long double x);
3. 反正切函数:* y = arctan(x)*
extern float atanf(float x);extern double atan(double x);extern long double atanl(long double x);
4. 2个参数的反正切函数:z = arctan(y/x)
extern float atan2f(float y, float x);extern double atan2(double y, double x);extern long double atan2l(long double y, long double x);
因为arctan的定义域是在(-∞, +∞),而值域是在(-?/2, ?/2)之间。因此 :atan2f(-1.0, 0.0) == -?/2; atan2f(1.0, 0.0) == ?/2;
这个函数提供的另外一个意义在于tan函数的值其实就是对边除以邻边的结果,因此当知道对边和邻边时就可以直接用这个逆三角函数来求得对应的弧度值。假如特殊情况下对边和邻边的值都是0.0,那么如果你调用atan(0.0/0.0)得到的值将是NAN而不是0。因为0.0/0.0的值是NAN,而对NAN调用atan函数返回的也是NAN,但是对atan2(0.0,0.0)调用返回的结果就是正确值0。
5. 余弦函数: y = cos(x)
extern float cosf(float x);extern double cos(double x);extern long double cosl(long double x);
6. 正弦函数:y = sin(x)
extern float sinf(float x);extern double sin(double x);extern long double sinl(long double x);
7. 正切函数:y = tan(x)
extern float tanf(float x);extern double tan(double x);extern long double tanl(long double x);
双曲函数
1. 反双曲余弦函数:y = arccosh(x)
extern float acoshf(float x);extern double acosh(double x);extern long double acoshl(long double x);
2. 反双曲正弦函数:y = arcsinh(x)
extern float asinhf(float x);extern double asinh(double x);extern long double asinhl(long double x);
3. 反双曲正切函数:y = arctanh(x)
extern float atanhf(float x);extern double atanh(double x);extern long double atanhl(long double x);
4. 双曲余弦函数:y = cosh(x)
extern float coshf(float x);extern double cosh(double x);extern long double coshl(long double x);
5. 双曲正弦函数:y = sinh(x)
extern float sinhf(float x);extern double sinh(double x);extern long double sinhl(long double x);
6. 双曲正切函数: y = tanh(x)
extern float tanhf(float x);extern double tanh(double x);extern long double tanhl(long double x);
指数函数
1. 自然常数e为基数的指数函数:y = e^x
extern float expf(float x);extern double exp(double x);extern long double expl(long double x);
2. 自然常数e为基数的指数减1:y = e^x - 1
extern float expm1f(float x);extern double expm1(double x); extern long double expm1l(long double x);
我们既然定义了exp函数,那么按理说要实现e^x-1
就很简单,为什么要单独定义这个函数呢?先看下面两个输出:
double o1 = exp(1.0e-13) - 1.0; double o2 = expm1(1.0e-13); printf("o1 = %e, o2 = %e", o1, o2);//output: o1 = 9.992007e-14, o2 = 1.000000e-13
从上面的例子中发现当用exp函数时出现了有效数字损失而expm1则没有。出现这种问题的原因就是浮点加减运算本身机制的问题,在浮点运算中下面两种类型的运算都有可能出现损失有效数字的情况:
两个相近的数相减
两个数量级相差很大的数字相加减
我们可以做一个实验,分别在调试器中查看a1,a2和b1,b2的结果:
double a1 = 5.37-5.36; double a2 = (5.37*100 - 5.36*100)/100;double b1 = 100.0-0.01; double b2 = (100.0/0.01 - 0.01/0.01)*0.01;//我们发现a1的值是0.0099999999999997868,而a2的值就是0.01//我们发现b1的值是99.989999999999994而b2的值是99.990000000000009
从上面的例子中可以看出当浮点数相近或者差异很大时加减运算出现了有效数字损失的情况,同时上面的例子也给出了一个减少这种损失的简易解决方案。再回到上面exp函数的场景中,因为exp(1.0e-13)的值和1.0是非常接近,因此当对这两个数做减法时就会出现有效数字损失的情况。我们再来考察expm1函数,这个函数主要用于当x接近于0时的场景。我们知道函数 y = e^x - 1 当x趋近于0时的极限是0,因此我们可以用泰勒级数来展开他:
e^x-1泰勒级数展开
可以看出这个级数收敛的很快,因此可以肯定的是expm1函数的内部实现就是通过上面的泰勒级数的方法来实现求值的。下面这段函数使用手册的文档也给出了用expm1代替exp函数的例子和说明:
Note that computations numerically equivalent to exp(x) - 1.0 are often hidden in more complicated expressions; some amount of algebraic manipu- lation may be necessary to take advantage of the expm1() function. Con- sider the following example, abstracted from a developer's actual produc- tion code in a bug report: double z = exp(-x/y)*(x*x/y/y + 2*x/y + 2) - 2 When x is small relative to y, this expression is approximately equal to: double z = 2*(exp(-x/y) - 1) and all precision of the result is lost in the computation due to cata- strophic cancellation. The developer was aware that they were losing precision, but didn't know what to do about it. To remedy the situation, we do a little algebra and re-write the expression to take advantage of the expm1() function: exp(-x/y)*(x*x/y/y + 2*x/y + 2) - 2 = (2*exp(-x/y) - 2) + exp(-x/y)*((x*x)/(y*y) + 2*x/y) This transformation allows the result to be computed to a high degree of accuracy as follows: const double r = x/y; const double emrm1 = expm1(-r); double z = 2.0*emrm1 + (1.0 + emrm1)*(2.0 + r)*r; It is not always easy to spot such opportunities for improvement; if an expression involving exp() seems to be suffering from an undue loss of accuracy, try a few simple algebraic operations to see if you can iden- tify a factor with the form exp(x) - 1.0, and substitute expm1(x) in its place.
3. 2为基数的指数函数:y = 2^x
extern float exp2f(float x);extern double exp2(double x); extern long double exp2l(long double x);
4. 浮点数构造函数:* y = x * 2^n*
extern float ldexpf(float x, int n);extern double ldexp(double x, int n);extern long double ldexpl(long double x, int n);
既然上面已经存在了一个exp函数,如果我们要实现相同的功能按理来只要:x*exp(n)
就好了,为什么还要单独提供一个新的ldexp函数呢?原因就是ldexp函数其实是一个用来构造浮点数的函数,我们知道浮点数的格式定义在IEEE754中,具体的结构为:符号*尾数*2^指数,刚好和ldexp所实现的功能是一致的,这里的x用来指定符号*尾数,而n则指定为指数。因此我们就可以借助这个函数来实现浮点数的构造。
5. 以FLT_RADIX基数的浮点数构造函数:y = x* FLT_RADIX^n
extern float scalbnf(float x, int n);extern double scalbn(double x, int n);extern long double scalbnl(long double x, int n);extern float scalblnf(float x, long int n);extern double scalbln(double x, long int n);extern long double scalblnl(long double x, long int n);
这里的FLT_RADIX是浮点数存储里面的基数(在float.h中有定义这个宏),一般情况下是2,这时候这个函数就和ldexp函数是一致的。但是有些系统的浮点数存储并不是以2为基数(比如IBM 360的机器)。因此如果你要构造一个和机器相关的浮点数时就用这个函数。
对数函数
1. 自然常数e为基数的对数函数:y = ln(x)
extern float logf(float x);extern double log(double x);extern long double logl(long double x);
2. 自然常数e为基数的对数函数: y = ln(x + 1)
extern float log1pf(float x);extern double log1p(double x);extern long double log1pl(long double x);
这个函数的使用场景主要用于当x趋近于0的情况,上面曾经描述过当两个浮点数之间的数量值相差很大时数字的加减会存在有效位丢失的情况。因此如果我们用log函数来计算时当x趋近于0的ln(x+1)时就会存在有效位的损失情况。比如下面的例子:
double o1 = log(1.0e-13 + 1); double o2 = log1p(1.0e-13); printf("o1 = %e, o2 = %e", o1, o2); //output: o1 = 9.992007e-14, o2 = 1.000000e-13
可以看出函数log1p主要用于当x接近于0时的场景。我们知道函数 y = ln(x+1) 当x趋近于0时的极限是0,因此我们可以用泰勒级数来展开他:
ln(x+1)的泰勒级数展开
可以看出这个级数收敛的很快,因此可以肯定的是log1p函数的内部实现就是通过上面的泰勒级数的方法来实现求值的。
3. 10为基数的对数函数:y = log10(x)
extern float log10f(float x);extern double log10(double x);extern long double log10l(long double x);
4. 2为基数的对数函数1:y = log2(x)
extern float log2f(float x);extern double log2(double x);extern long double log2l(long double x);
5. FLT_RADIX为基数的对数函数并取整:y = floor(log2(x))
extern float logbf(float x);extern double logb(double x);extern long double logbl(long double x);
函数返回的是一个小于等于真实指数的最大整数,也就是对返回的值进行了floor操作,具体floor函数的定义见下面。这里的FLT_RADIX是浮点数的基数,大部分系统定义为2。下面是这个函数的一些例子:
logb(2.5) == floor(log2(2.5)) == 1; logb(4.0) == floor(log2(4.0)) == 2; logb(4.1) == floor(log2(4.1)) == 2; logb(7) == floor(log2(7)) == 2; logb(7.9999) == floor(log2(7.9999)) == 2; logb(8.0) == floor(log2(8.0)) == 3;
6. FLT_RADIX为基数的对数函数并取整:y = floor(log2(x))
extern int ilogbf(float x);extern int ilogb(double x);extern int ilogbl(long double x);
函数返回的是一个小于等于真实指数的最大整数,也就是对返回的值进行了floor操作,具体floor函数的定义见下面。需要注意的是这里返回的类型是整型,因此不可能存在返回NAN或者** INFINITY**的情况。下面是当x是0或者负数时返回的特殊值:
FP_ILOGB0: 当x是0时返回这个特殊值。 FP_ILOGBNAN:当x是负数时返回这个特殊值。
这里区分一下log2,logb,ilogb 这三个函数的差异:
logb,ilogb是以FLT_RADIX为基数的对数,而log2则是以2为基数的对数,虽然大部分系统中FLT_RADIX默认是定义为2。
log2,logb返回的都是浮点型,因此有可能返回INFINITY和NAN这两个特殊值;而ilogb则返回的是整型,因此如果x是特殊的话那么将会返回FP_ILOGB0和FP_ILOGBNAN两个值。
log2返回的是有可能带小数的指数,而logb和ilogb则返回的是一个不大于实际指数的整数。
绝对值函数
1. 取绝对值函数:y = |x|
extern float fabsf(float);extern double fabs(double);extern long double fabsl(long double);
幂函数
1. 平方根函数:y = √x
extern float sqrtf(float x);extern double sqrt(double x);extern long double sqrtl(long double x);
2. 立方根函数: y = ∛x
extern float cbrtf(float x);extern double cbrt(double x);extern long double cbrtl(long double x);
3. 幂函数:z = x ^ y
extern float powf(float x, float y);extern double pow(double x, double y);extern long double powl(long double x, long double y);
4. 欧几里得距离函数: *d =√x2+y2 *
extern float hypotf(float x, float y);extern double hypot(double x, double y);extern long double hypotl(long double x, long double y);
这个函数可以用来求直角三角形的斜边长度。
误差函数
误差函数主要用于概率论和偏微分方程中使用,具体参考误差函数
1. 误差函数
extern float erff(float x);extern double erf(double x);extern long double erfl(long double x);
2. 互补误差函数
extern float erfcf(float x);extern double erfc(double x);extern long double erfcl(long double x);
伽玛函数
1. 伽玛函数 :y = ?(x)
extern float lgammaf(float x);extern double lgamma(double x);extern long double lgammal(long double x);
2. 阶乘函数:y = (x-1)!
extern float tgammaf(float x);extern double tgamma(double x);extern long double tgammal(long double x);
伽玛函数其实就是阶乘在实数上的扩展,一般我们知道3! = 3*2*1 = 8。那么我们要求2.5!怎么办,这时候就可以用这个函数来实现。这个函数也可以用来进行阶乘计算。 注意这里是x-1后再计算的。
取整函数
1. 返回一个大于等于x的最小整数
extern float ceilf(float x);extern double ceil(double x);extern long double ceill(long double x);
举例来说我们要对于一个负浮点数按0.5进行四舍五入处理:即当某个负数的小数部分大于等于0并且小于0.5时则舍弃掉小数部分,而当小数部分大于等于0.5并且小于1时则等于0.5。我们就可以用ceil函数来实现如下:
double y = ceil(x*0.5)/0.5;
2. 返回一个小于等于x的最大整数
extern float floorf(float x);extern double floor(double x);extern long double floorl(long double x);
举例来说我们要对于一个正浮点数按0.5进行四舍五入处理:即当某个正数的小数部分大于等于0并且小于0.5时则舍弃掉小数部分,而当小数部分大于等于0.5并且小于1时则等于0.5。我们就可以用floor函数来实现如下:
double y = floor(x*0.5)/0.5;
3. 返回一个最接近x的整数
extern float nearbyintf(float x);extern double nearbyint(double x);extern long double nearbyintl(long double x);extern float rintf(float x);extern double rint(double x);extern long double rintl(long double x);//下面三个函数返回的是整数。extern long int lrintf(float x);extern long int lrint(double x);extern long int lrintl(long double x);//下面三个函数是C99或者gnu99中的函数。extern long long int llrintf(float x);extern long long int llrint(double x);extern long long int llrintl(long double x);
上述各函数的区别请参考:http://zh.cppreference.com/w/c/numeric/math/rint
4. 对x进行四舍五入取整
extern float roundf(float x);extern double round(double x);extern long double roundl(long double x);extern long int lroundf(float x);extern long int lround(double x);extern long int lroundl(long double x);//下面三个函数是C99或者gnu99中的函数。extern long long int llroundf(float x);extern long long int llround(double x);extern long long int llroundl(long double x);
如果x是正数,那么当小数部分小于0.5则返回的整数小于浮点数,如果小数部分大于等于0.5则返回的整数大于浮点数;如果x是负数,那么当小数部分小于0.5则返回的整数大于浮点数,如果小数部分大于等于0.5则返回的整数小于浮点数。
** 如果我们要实现保留N位小数的四舍五入时。我们可以用如下的方法实现:**
double y = round(x * pow(10, N)) / pow(10, N)
数字拆分
1. 返回浮点数x的整数部分
extern float truncf(float x);extern double trunc(double x);extern long double truncl(long double x);
这个函数和floor函数的区别主要体现在负数上,对一个负数求floor则会返回一个小于等于负数的负整数,而对一个负数求trunc则会返回一个大于等于负数的负整数。
如果我们要实现保留N位小数的截取时。我们可以用如下的方法实现:
double y = trunc(x * pow(10, N)) / pow(10, N)
2. 返回x/y的余数1: z = mod(x, y)
extern float fmodf(float x, float y);extern double fmod(double x, double y);extern long double fmodl(long double x, long double y);
函数返回值r = x - n*y, 其中n等于x/y的值截取的整数。
3. 返回x/y的余数2: z = mod(x, y)
extern float remainderf(float x, float y);extern double remainder(double x, double y);extern long double remainderl(long double x, long double y);
函数返回值r = x - n*y, 其中n等于x/y的值取最接近的整数,如果有两个数都接近x/y,那么n就取偶数。比如我们要求remainder(7,2)
。因为7/2是3.5,按上面规则n就取4,因此最后的结果是r = 7 - 4*2 = -1。同样我们可以得出remainder(7,3) == 7-2\*3 == 1
。
从上面的描述可以看出fmod和remainder的区别主要在于x/y的整数部分的处理不一样:前者是取x/y的整数来算余数,而后者则取最接近x/y的整数来算余数。
4. 返回x/y的余数和整数商
extern float remquof(float x, float y , int *quo);extern double remquo(double x, double y, int *quo);extern long double remquol(long double x, long double y, int * quo);
这个函数和** remainder**函数一样,只不过会将整数商也返回给quo,也就是说r = x - n *y这个等式中,r作为函数的返回,而n则返回给quo。
5. 分解出x的整数和小数部分
extern float modff(float x, float p*);extern double modf(double x, double p*);extern long double modfl(long double x, long double p*);
函数返回小数部分,整数部分存储在p中。这里面返回值和p都和x具有相同的符号。
6. 分解出x的指数和尾数部分
extern float frexpf(float x, int * p);extern double frexp(double x, int * p);extern long double frexpl(long double x, int * p);
函数返回尾数*符号部分,指数部分存储在p中。需要明确的是如果浮点数x为0或者非规格化浮点数时按浮点数的定义格式返回尾数和指数,而当x为规格化浮点数那么返回的值的区间是[0.5, 1)。这里的返回值和指数值p和上面介绍的规格化浮点数格式:** 符号 * (1.尾数) * 2^指数 有差异。因为按照定义返回的尾数部分应该是1.xxx,但是这里的返回值却是[0.5, 1)。其实这并不矛盾,只是函数对返回的值做了特殊处理:因为一个正浮点数可以表示为:1.m * 2^e ==> (2^0 + 0.m) * 2^e ==> (2^0 / 2 + 0.m / 2) *2^(e+1) =>(0.5 + 0.m/2) *2^(e+1)。因此frexp函数返回的真实值是: 尾数除以2,而p存储的是:指数+1**
下面函数使用的一些例子:
int p1 = 0; double y1 = frexp(16.0, &p); //y1=0.5, p= 5 int p2 = 0; double y2 = frexp(1.0, &p); //y2=0.5, p = 1 int p3 = 0; double y3 = frexp(0.0, &p); //y3=0, p = 0
这个函数和上面的ldexp函数为互逆函数。要详细的了解浮点数存储格式请参考IEEE754
符号改变
1. 将y的符号赋值给x并返回具有和y相同符号的x值
extern float copysignf(float x, float y);extern double copysign(double x, double y);extern long double copysignl(long double x, long double y);
举例如下:
copysign(10.0, 9.0) == 10; copysign(-10.0, -9.0) == -10; copysign(-10.0, 9.0) == 10; copysign(10.0, -9.0) == -10;
这个函数的作用是实现符号的赋值,有就是将y的符号赋值给x。
无效数字定义
1.生成一个quient NAN浮点数
extern float nanf(const char *tagp);extern double nan(const char *tagp);extern long double nanl(const char *tagp);
前面我有介绍了浮点数里面有两个特殊的值:无穷INFINITY和非法NAN,既然这两个数字都可以用浮点数来描述,那么他就肯定也有对应的存储格式。我们知道浮点数的格式为:符号*尾数*2^指数。在IEEE754标准中就对无穷和非法这两种特殊的数进行了定义:
当浮点数中的指数部分的二进制位全为1。而尾数部分的二进制位全为0时则表示的浮点数是无穷INFINITY,如果符号位为0则表示正无穷大,而符号位为1则表示负无穷大。
当浮点数中的指数部分的二进制位全为1。而尾数部分的二进制位不全为0时则表示的浮点数是非法数字NAN,或者表示为未定义的数字。
从上面的对NAN的定义可以得出非法数字并不是一个具体的数字而是一类数字,因此对两个为NAN的浮点数字并不能用等号来比较。以32位IEEE单精度浮点数的NAN为例,按位表示即:S111 1111 1AXX XXXX XXXX XXXX XXXX XXXX,其中的S是符号位,而符号位后面的指数位为8个1表示这个数字是一个特殊的浮点数,剩余的A和X则组成为了尾数部分,因为是NAN 所以我们要求A和X这些位中至少有一个是1。在IEEE 754-2008标准中,又对NAN的类型进行了细分:
如果A = 1,则该数是quiet NAN。也就是quiet NAN中尾数的最高位为1。
如果A为零、其余X部分非零,则是signaling NAN。
区分两种NAN的目的是为了更好的对浮点数进行处理。一般我们将signaling NAN来表示为某个数字未初始化,而将quiet NAN则用来表示浮点运算的结果出现了某类异常,比如0除异常,比如负数开根异常等等。既然quiet NAN可以用来对无效数字进行分类,也就是说我们可以构建出一个有类别标志的quiet NAN。因此nan函数就是一个专门构建具有无效类别的NAN函数(绕了这么多终于说到点子上了)。nan函数中的tagp参数就是用来指定非法数字中的类别,虽然参数类型是字符串,但是要求里面的值必须是整数或者空字符串,而且系统在构造一个quiet NAN时会将tagp所表示的整数放在除A外的其他尾数位上。下面是使用nan函数的例子:
float f1 = NAN; //0b01111111110000000000000000000000 float f2 = nanf(""); //0b01111111110000000000000000000000 float f3 = nanf("123"); //0b01111111110000000000000001111011 float f4 = nanf("456"); //0b01111111110000000000000111001000 float f5 = nanf("abc"); //0b01111111110000000000000000000000
具体操作时我们可以用如下来方法来处理各种异常情况:
//定义部分:float testfn(){ //有异常时根据不同的情况返回不同的nan。 if (异常1) return nan("100"); else if (异常2) return nan("200");else return 正常数字; }//调用部分:float ret = testfn();if (isnan(ret)) { //取非法数字的错误标志部分 int exceptionType = ret & 0x3FFFFF; if (exceptionType == 100) { } else if (exceptionType == 200) { } }else{ //正常处理。}
有一个地方疑惑的是为什么NAN定义默认值是一个quiet NAN而不是signaling NAN
?递增函数
1. 返回x在y方向上的下一个可表示的浮点数。
extern float nextafterf(float x, float y);extern double nextafter(double x, double y);extern long double nextafterl(long double x, long double y);extern double nexttoward(double x, long double y);extern float nexttowardf(float x, long double y);extern long double nexttowardl(long double x, long double y);
如果x等于y则返回x。这个函数主要用来实现那些需要高精度增量循环的处理逻辑。也就是说如果对浮点数进行for循环处理时,这个函数可以用来实现最小的浮点数可表示的数字的增量。比如下面的代码:
for (double x = 0.1; x < 0.2; x=nextafter(x,0.2)) { //... }
注意这里是下一个可表示的浮点数,也就是说当x为0而y为1时,那么返回的值将是最小的非常规浮点数;而如果x为1而y为2时,那么返回的值将是1+DBL_MIN(or FLT_MIN). 下面是具体的示例代码:
// 0.0f == 0b00000000000000000000000000000000 float a = nextafterf(0.0f, 1.0f); //a == 0b00000000000000000000000000000001 // FLT_MIN == 0b00000000100000000000000000000000 float b = nextafterf(FLT_MIN, 1.0f); // b = 0b00000000100000000000000000000001 // 1.0f == 0b00111111100000000000000000000001 float c = nextafterf(1.0f, 1.1f); // c = 0b00111111100000000000000000000001
?比较函数
1. 返回x减去y的差如果x>y,否则返回0
extern float fdimf(float x, float y);extern double fdim(double x, double y);extern long double fdiml(long double x, long double y);
这个函数可以用来求两个数的差,并且保证不会出现负数。下面是使用的例子:
double a = fdim(5.0, 3.0); //2.0 double b = fdim(5.0, 5.0); //0.0 double c = fdim(5.0, 6.0); //0.0
2. 返回x和y中大的数字: z = max(x,y)
extern float fmaxf(float x, float x);extern double fmax(double x, double x);extern long double fmaxl(long double x, long double x);
3. 返回x和y中小的数字: z = min(x,y)
extern float fminf(float x, float y);extern double fmin(double x, double y);extern long double fminl(long double x, long double y);
浮点乘加运算
1. 浮点乘加运算:w = x*y + z
extern float fmaf(float x, float y, float z);extern double fma(double x, double y, double z);extern long double fmal(long double x, long double y, long double z);
这个函数返回x*y+z的结果,而且会保证中间计算不会丢失精度。这个函数会比直接用x*y+z要快,因为CPU中专门提供了一个用于浮点数乘加的指令FMA。具体情况请参考关于浮点乘加器方面的资料和应用。
作者:欧阳大哥2013
共同学习,写下你的评论
评论加载中...
作者其他优质文章