计算机是如何计算对数函数 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(22m)=21ln2+ln(22m)(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})
22m∈(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(22m)x=m⋅2j−127(6)
对于
ln
(
2
2
m
)
\ln{(\frac{\sqrt{2}}{2}m)}
ln(22m),使用泰勒公式将其展开:
令
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=(22m)则ln(22m)=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 = 130
i_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<stdio.h>
#include <math.h>
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 < 8; i++){
k = -1.0 * k;
result += k * t_m / (1.0 * i);
t_m *= m;
}
return x_j * ln_2 + 0.5 * ln_2 + result;
}
int main(){
for (int i = 1; i <= 30; i++){
float x = i * 0.2;
printf("x = %2.6f, my log(x) = %2.6f, log(x) = %2.6f, err = %2.6f\n", x, my_log(x), log(x), fabs(my_log(x) - log(x)));
}
return 0;
}
误差
运行结果如下