在 Python 中计算以 10 为底的对数的算法

Question

我已经尝试创建一个程序来根据 "The Mathematical-Function Computation Handbook" 中描述的基于泰勒级数的算法来计算以 10 为底的对数（我通过我大学的图书馆找到了一个在线副本）。

Whosebug 上的另一个问题给出了类似的算法，我现在找不到 link。

10.3.2 Computing logarithms in a decimal base

For a decimal base, the base-10 logarithm is the natural choice, and the decomposition of the argument into an exponent and a fraction gives us a decimal representation: x = (−1)^s × f × 10^n, either f = 0 exactly, or f is in [1/10, 1).

If f ≤√1/10, set f = 10 × f and n = n − 1, so that f is now in the interval (√1/10,√10]. Then introduce a change of variable, a Taylor-series expansion, and a polynomial representation of that expansion:

z = (f − 1)/( f + 1),

f = (1 + z)/(1 − z),

D = 2 log10(e)

= 2/ log(10)

log10( f) = D × (z + z3/3 + z5/5 + z7/7 + z9/9 + z11/11 + · · · )

≈ D × z + z3Q(z2), polynomial fit incorporates D in Q(z2).

For f in (√1/10,√10], we have z in roughly [−0.5195,+0.5195]. The wider range of z requires longer polynomials compared to the binary case, and also makes the correction term z3Q(z2) relatively larger. Its magnitude does not exceed |0.35z|, so it provides barely one extra decimal digit of precision, instead of two. Accurate computation of z is easier than in the binary case: just set z = fl(fl(f−12)−12)/fl(f+1).

为此，我在 Python 中编写了这个程序：

def log10(x):

n = 0.0 #Start exponent of base 10

while (x >= 1.0):
    x = x/10.0
    n+=1


# if x <= sqrt(1/10)
if(x<=0.316227766016838):
    x = x*10.0
    n = n-1

#Produce a change of variable
z = (x-1.0)/(x+1.0)
D = 4.60517018598809 #2*log10(e)

sum = z
for k in range(3,111,2):
    sum+=(z**k)/k

return D*n*sum

我将结果与math.log10函数进行了比较，结果并不如预期。调试时我最大的问题是理解算法及其工作原理。

Answer 1

这是我在建议更正后的源代码（将 return 语句更改为 D*sum+n 固定了 D 的值，并将 if(x<=0.316227766016838) 更改为 while(x<=0.316227766016838). 我添加了一些 if 语句来处理异常情况。

下面的代码在我的 6 位数目标精度内运行良好（我用非常小的输入、大的输入测试了它）。

def log10(x):

    # Handle exceptional cases
    if (x == 1):
        return 0
    if (x == 0):
        return float('-Inf')
    if (x < 0):
        return float('nan')

    n = 0 #Start exponent of base 10

    while (x >= 1.0):
        x = x/10.0
        n+=1

    # if x <= sqrt(1/10)
    while(x<=0.316227766016838):
        x = x*10.0
        n = n-1

    #Produce a change of variable
    z = (x-1.0)/(x+1.0)
    D = 0.868588964 #2*log10(e)

    #Taylor series
    sum = z
    for k in range(3,23,2):
        sum+=(z**k)/k

    return D*sum+n

在 Python 中计算以 10 为底的对数的算法

Algorithm for computing base-10 logarithm in Python

python

algorithm

math

elementary-functions