计算机如何计算对数函数 |
您所在的位置:网站首页 › 计算器如何计算ln的 › 计算机如何计算对数函数 |
计算机是如何计算对数函数 log a x \log_a{x} logax 的呢? 本文并不是为了探究各种编程语言中数学计算函数库对 log \log log 函数的具体实现,而是提供一种计算思路,来帮助理解 log \log log 函数的计算原理,并使用 C 语言来实现这个算法。 计算原理 泰勒展开提起对数函数 log a x \log_a{x} logax,大家最熟悉的对数函数应该是以自然底数 e e e 为底的对数函数: ln x \ln{x} lnx,即 log e x \log_e{x} logex。 ln x \ln{x} lnx 函数的图像如下所示: 对于不是以 e e e 为底的对数函数 log a x \log_a{x} logax ,可以使用换底公式,将其转换为以 e e e 为底的函数形式。 换底公式: log a b = log c b log c a 换底公式:\log_a{b} = \frac{\log_c{b}}{\log_c{a}} \qquad\qquad\qquad 换底公式:logab=logcalogcb 根据换底公式: log a x = log e x log e a = ln x ln a (1) \log_a{x} = \frac{\log_e{x}}{\log_e{a}} = \frac{\ln{x}}{\ln{a}} \tag{1} \qquad\qquad\qquad logax=logealogex=lnalnx(1) 于是,计算对数函数 log a x \log_a{x} logax 的问题就转换为计算 ln x \ln{x} lnx 的问题。 对于计算 ln x \ln{x} lnx,可能我们第一时间想到的就是泰勒展开公式: ln ( x ) = ∑ n = 0 ∞ ( − 1 ) n n + 1 ( x − 1 ) n + 1 = ( x − 1 ) − 1 2 ( x − 1 ) 2 + 1 3 ( x − 1 ) 3 − 1 4 ( x − 1 ) 4 + 1 5 ( x − 1 ) 5 + ⋅ ⋅ ⋅ , x ∈ ( 0 , 2 ] (2) \ln{(x)} = \sum_{n = 0}^{\infty} \frac{(-1)^n}{n+1}{(x-1)}^{n+1} = (x - 1) - \frac{1}{2}{(x-1)}^2 + \frac{1}{3}{(x-1)}^3 - \frac{1}{4}{(x-1)}^4 + \frac{1}{5}{(x-1)}^5 + ···, \quad x \in (0, 2] \tag{2} ln(x)=n=0∑∞n+1(−1)n(x−1)n+1=(x−1)−21(x−1)2+31(x−1)3−41(x−1)4+51(x−1)5+⋅⋅⋅,x∈(0,2](2) 但是,使用泰勒公式对 ln x \ln{x} lnx 展开是有条件的,即 x x x 所属的区间为 x ∈ ( 0 , 2 ] x \in (0, 2] x∈(0,2] 。而且,展开的项数一定是固定的,毕竟不能无限次的进行展开计算。 要知道的是,使用泰勒公式对 ln x \ln{x} lnx 进行展开,始终是一种近似,而非绝对的等于。决定展开后的误差的因素有两点: 一: x x x 与 1.0 1.0 1.0 的距离。在展开项数相同的情况下, x x x 越接近 1.0 1.0 1.0,误差越小; 二:展开的项数。对于同一个 x x x,展开项数越多,误差越小。 上图就是 ln x \ln{x} lnx 与不同项数的泰勒展开式的图像。 作为对比,下图列出了上述几个展开式与 ln x \ln{x} lnx 之间的误差: [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-LqjbDJGO-1676008929754)(计算机如何计算对数函数/ln与展开式误差.svg)] 结合上面两张图,可以看出, x x x 越趋向于 1.0 1.0 1.0 ,展开项数越多,精度越高。因此,要想办法把 x x x 转化到 1 1 1 的附近,就能够对 x x x 使用泰勒展开的方式进行计算了。 根据公式 ( 1 ) (1) (1) ,我们成功将 log a x \log_a{x} logax 转化为求 ln x \ln{x} lnx 的问题,对于 ln x \ln{x} lnx ,我们希望能够使用泰勒公式对其进行展开。但是由于不能确定 x x x 的值,不能直接使用泰勒公式进行展开,那么我们就要想办法处理 x x x,让其满足泰勒公式展开的条件,并且尽可能的提高计算精度。 单精度存储在处理 x x x 之前,我们要先弄清楚,数据在计算机中是如何存储的。 在计算机中,使用浮点数来表示一个带有小数的数(整数也可以使用浮点数来表示)。浮点数使用 IEEE 格式存储,分为单精度浮点数(float)、双精度浮点数(double)。 对于单精度浮点数而言,占 4 个字节,即 32 位。其中使用 1 位表示符号位,8 位表示二进制形式的阶码,剩余 23 位表示尾巴数。单精度浮点数的表示范围大约在 3. 4 − 38 − 3. 4 38 3.4^{-38} - 3.4^{38} 3.4−38−3.438 之间。 其中,最高位为符号位 f f f,符号位 f = 0 f = 0 f=0 代表该数是正数, f = 1 f = 1 f=1 代表该数是负数。 中间 8 位表示阶码 j j j ,阶码部分为纯整数,使用移码表示,真实值为 j − 127 j-127 j−127。 低 23 位表示尾数 m m m,在构造浮点数据的时候需要进行规格化,即要保证尾数的整数部分为 1,而这个 1 是不储存的,所以经过规格化后 m ∈ ( 1 , 2 ) m \in(1,2) m∈(1,2)。 因此采用 IEEE 格式保存的数据的值为: ( − 1 ) f ⋅ m ⋅ 2 j − 127 (-1)^f \cdot m \cdot 2^{j-127} (−1)f⋅m⋅2j−127 。 以 13.75 为例: 1、将十进制数转换为二进制数: 13.75 = 1101.11 B 13.75 = 1101.11B 13.75=1101.11B,即 13 = 8 + 4 + 1 = 2 3 + 2 2 + 2 0 = 1101 B 13 = 8 + 4 + 1 = 2^3 + 2^2 + 2^0 = 1101B 13=8+4+1=23+22+20=1101B, 0.75 = 0.5 + 0.25 = 2 − 1 + 2 − 2 = 0.11 B 0.75 = 0.5 + 0.25 = 2^{-1} + 2^{-2} = 0.11B 0.75=0.5+0.25=2−1+2−2=0.11B,因此 13.75 = 13 + 0.75 = 1101 B + 0.11 B = 1101.11 B 13.75 = 13 + 0.75 = 1101B + 0.11B= 1101.11B 13.75=13+0.75=1101B+0.11B=1101.11B; 2、将该二进制数规格化: 1101.11 = 1.10111 × 2 3 1101.11 = 1.10111 \times 2^3 1101.11=1.10111×23,尾数的整数位是隐藏的,因此尾数部分存储的是: 10111 10111 10111 。 3、将阶码转换为移码表示: 3 + 127 = 130 = 11 B + 01111111 B = 10000010 B 3 + 127 = 130 = 11B + 01111111B = 10000010B 3+127=130=11B+01111111B=10000010B,所以阶码中存储的数据为: 100000010 100000010 100000010。 因此,在计算机中, 13.75 存储为: 双精度浮点数(double)存储形式与单精度浮点数(float)相似,只不过双精度浮点数占 8 个字节,64 位,最高位代表符号位,中间 11 位为阶码,剩余 52 位表示尾数。 归约化借助 float 数据的存储方式,我们将 x x x 进行转换: x = ( − 1 ) f ⋅ m ⋅ 2 j − 127 x = (-1)^f \cdot m \cdot 2^{j-127} x=(−1)f⋅m⋅2j−127,由于 ln x \ln{x} lnx 函数的定义域为 ( 0 , + ∞ ) (0,+\infty) (0,+∞) ,因此阶码始终为 0,于是,化简为 x = m ⋅ 2 j − 127 x = m \cdot 2^{j-127} x=m⋅2j−127。 借助对数函数的另外两个性质: log a ( b ⋅ c ) = log a ( b ) + log a ( c ) log a ( b j ) = j ⋅ log a ( b ) (3) \log_a{(b \cdot c)} = \log_a{(b)} + \log_a{(c)} \\ \log_a{(b^j)} = j \cdot \log_a{(b)} \tag{3} loga(b⋅c)=loga(b)+loga(c)loga(bj)=j⋅loga(b)(3) 可得: ln x = ln [ ( 1. m ) ⋅ 2 j − 127 ] = ln ( 2 j − 127 ) + ln ( m ) = ( j − 127 ) ⋅ ln 2 + ln m (4) \ln{x} = \ln{[(1.m) \cdot 2^{j-127}]} = \ln{(2^{j-127})} + \ln{(m)} = (j-127) \cdot \ln{2} + \ln{m} \tag{4} lnx=ln[(1.m)⋅2j−127]=ln(2j−127)+ln(m)=(j−127)⋅ln2+lnm(4) 在 (4) 式中, ( j − 127 ) (j-127) (j−127) 为 x x x 的阶码,可以想办法取; ln 2 ≈ 0.693147 \ln2 \approx 0.693147 ln2≈0.693147 为常数,因此,只需要计算 ln m \ln{m} lnm 的值即可。 又因为 m m m 为 x x x 的尾数,根据规格化的规则, m ∈ ( 1 , 2 ) m \in(1,2) m∈(1,2),符合对其进行泰勒展开的要求。 为了展开的精度更加精确,我们应该尽量让 m m m 趋向于 1 1 1,根据 (3) 式: ln ( m ) = ln ( 2 ⋅ 2 2 ⋅ m ) = ln 2 + ln ( 2 2 m ) = 1 2 ln 2 + ln ( 2 2 m ) (5) \ln{(m)} = \ln{(\sqrt{2} \cdot\frac{\sqrt{2}}{2} \cdot m)} = \ln{\sqrt{2}} + \ln{(\frac{\sqrt{2}}{2} m)} = \frac{1}{2}\ln{2} + \ln{(\frac{\sqrt{2}}{2} m)} \tag{5} ln(m)=ln(2 ⋅22 ⋅m)=ln2 +ln(22 m)=21ln2+ln(22 m)(5) 因为 m ∈ ( 1 , 2 ) m \in(1,2) m∈(1,2),所以 2 2 m ∈ ( 2 2 , 2 ) \frac{\sqrt{2}}{2} m \in(\frac{\sqrt{2}}{2},\sqrt{2}) 22 m∈(22 ,2 ),即 m m m 的大致范围为 ( 0.707106 , 1.414213 ) (0.707106, 1.414213) (0.707106,1.414213),更加接近于 1 1 1。 于是,我们最终的式子就是: ln x = ( j − 127 ) ⋅ ln 2 + 1 2 ln 2 + ln ( 2 2 m ) x = m ⋅ 2 j − 127 (6) \ln{x} = (j-127) \cdot \ln{2} + \frac{1}{2}\ln{2} + \ln{(\frac{\sqrt{2}}{2} m)} \\ x = m \cdot 2^{j-127} \tag{6} lnx=(j−127)⋅ln2+21ln2+ln(22 m)x=m⋅2j−127(6) 对于 ln ( 2 2 m ) \ln{(\frac{\sqrt{2}}{2}m)} ln(22 m),使用泰勒公式将其展开: 令 t = ( 2 2 m ) 则 ln ( 2 2 m ) = ln ( t ) = ( t − 1 ) − 1 2 ( t − 1 ) 2 + 1 3 ( t − 1 ) 3 − 1 4 ( t − 1 ) 4 + 1 5 ( t − 1 ) 5 − 1 6 ( t − 1 ) 6 + 1 7 ( t − 1 ) 7 + o ( t − 1 ) 7 (7) 令 \quad t = (\frac{\sqrt{2}}{2}m) \\ 则 \quad \ln{(\frac{\sqrt{2}}{2}m)} = \ln{(t)} = (t-1) - \frac{1}{2}(t-1)^2 + \frac{1}{3}(t-1)^3 -\frac{1}{4}(t-1)^4 + \frac{1}{5}(t-1)^5 - \frac{1}{6}(t-1)^6 + \frac{1}{7}(t-1)^7 + o(t-1)^7 \tag{7} 令t=(22 m)则ln(22 m)=ln(t)=(t−1)−21(t−1)2+31(t−1)3−41(t−1)4+51(t−1)5−61(t−1)6+71(t−1)7+o(t−1)7(7) 实现思路根据 (6) 式,要想完成 ln x \ln{x} lnx 的计算,我们要先知道 x x x 的阶码 j j j、尾数 m m m 和 ln 2 \ln2 ln2 的值。其中 ln 2 \ln2 ln2 为常数,约等于 0.693147 0.693147 0.693147,那么剩下就是要通过 x x x 取得 j j j 和 m m m。 我们已经知道了浮点数 x x x 在内存中的二进制存放方式,因此只需要通过简单的位运算,就可以分别获取 j j j 和 m m m 了。 首先,计算机是不支持 float 型数据的位运算的,因此我们要将其内存中的数据以 int 型数据的存放方式读取出来,然后再对其进行位运算。 以 int 型数据格式读取 float 型二进制数据的方式如下: float x = 13.75; int i_x = *(int *)&x; // i_x = 0x415c0000其中 x 是我们想要计算的数,设其值为 13.75; i_x 是以 int 型数据读取的出的 x 的二进制数据,其十六进制的值为 0x415c0000,即: 为了取阶码,我们构造一个二进制 int 型数据 0x7f800000,其二进制表示为: 然后与 i_x 作 与 运算: int i_j = i_x & 0x7f800000; // i_j = 0x41000000将 i_j 右移 23 位: i_j = i_j >> 23; // i_j = 0x00000082 = 130i_j 减去移码的 127,就得到 x 的阶码: i_j = i_j - 127; // i_j = 3获取尾数的操作类似,我们先构造一个二进制 int 型数据 0x007fffff,其二进制表示为: 将其与 i_x 作 与 运算: int i_m = i_x & 0x007fffff; // i_m = 0x005c0000取到了尾数后,通过前面知道,浮点型的数据 x = m ⋅ 2 j − 127 x = m \cdot 2^{j-127} x=m⋅2j−127,而我们现在只取到尾数的二进制数据,需要将其构造出来,那么还需要将其阶码补齐。在此构造一个二进制 int 型数据 0x3f800000,其二进制表示如下: 将其与 i_m 作 或 运算: i_m = i_m | 0x3f800000; //i_m = 0x7fdc0000再以 float 的格式将其读取出来: float x_m = *(float *)&i_m;后续就是将其进行缩放与泰勒展开,在此不再赘述。 代码实现 源码 #include #include float my_log(float x){ float ln_2 = 0.693147; // ln(2) = 0.693147 int i_x = *(int*)&x; int x_j = ((i_x & 0x7f800000) >> 23) - 127; int x_m = (i_x & 0x007fffff) | 0x3f800000; float m = *(float*)&x_m; m *= 0.707106; // sqrt(2) / 2 = 0.707106 m -= 1.0; float t_m = m * m; float k = 1.0; float result = m; for (int i = 2; i |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |