1

I have tried to create a program to calculate the base-10 logarithm based on the Taylor series-based algorithm described in "The Mathematical-Function Computation Handbook" (I found an online copy via my University's library).

A similar algorithm is given on another question on StackOverflow for which I cannot find the link right now.

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).

For this I wrote this program in 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

I compared the results to the math.log10 function, and the results are not as expected. My biggest issue when debugging is understanding the algorithm and why it works.

A.G.
  • 37
  • 8
  • 1
    I see some minor errors. Can n become 0? Does D*n*sum look like the correct expression in that case? And 2*log10(e) is not the number you have in the code. – Arndt Jonasson Mar 18 '18 at 23:54
  • If you're trying to get x into the range `(sqrt(1/10), 1)`, shouldn't you need a while on the lower bound, just like on the upper bound? Also, the algorithm you quoted specifies bounds of `(sqrt(1/10), sqrt(10)]`; you are you using different ones? – abarnert Mar 19 '18 at 00:00
  • And that final step of `D*n*sum` can't be the right use of `n`. That would mean that the log of every value in `(sqrt(1/10), 1)` is 0. Presumably you're supposed to add it, not multiply it? – abarnert Mar 19 '18 at 00:00
  • @ArndtJonasson Thank you for pointing out the incorrect value of `2*log10(e)`. I must have made a typo when calculating it in Excel. I have changed it to `0.868588964`. @abarnert You are absolutely right. The `n` term should be added. I have changed the return statement to `return D*sum+n`. With respect to your comment on the lower bound, I believe the issue of the `f` variable being too small is already handled by the `if(x<=0.316227766016838)` statement. I am debugging the program right now, but it looks much better with those two fixes. Thanks! – A.G. Mar 19 '18 at 00:16
  • @abarnert I now understand what you meant. Changed the `if` statement with a `while` statement. – A.G. Mar 19 '18 at 00:30
  • @dreamerboy there are also another options to compute this see [Building a logarithm function in C without using float type](https://stackoverflow.com/a/42108287/2521214) – Spektre Mar 19 '18 at 11:34
  • @Spektre Thanks. – A.G. Mar 20 '18 at 00:19

1 Answers1

1

Here is my source code after the suggested corrections (changed return statement to D*sum+n fixed the value of D, and changed if(x<=0.316227766016838) to while(x<=0.316227766016838). I added some if statements to handle exceptional cases.

The code below works well within my target precision of 6 digits (I tested it with very small input, large input).

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
A.G.
  • 37
  • 8
  • Your code as shown has the obvious indent error after the first line. – Ancora Imparo Mar 19 '18 at 00:28
  • @AncoraImparo I don't get an indent error when running it. – A.G. Mar 19 '18 at 00:32
  • Yes, you have re-edited. Your original answer didn't have the indents. – Ancora Imparo Mar 19 '18 at 00:50
  • 1
    So, does this work (within the error limits expected for your approximation) on a range of values? If so, congrats. Now edit the answer to explain that, and to explain what you did and why (people aren't expected to read the comments to understand an answer, especially since comments occasionally get deleted). Then you can accept your own answer, and you'll help future seekers and get some rep points. – abarnert Mar 19 '18 at 02:46
  • @abarnert great. Updated answer. – A.G. Mar 20 '18 at 00:12