math.log(243, 3) should result in 5.0, but instead results in 4.9999...

Can be worked around via base conversion (math.log10(243) / math.log10(3))

Issue seems to be with all power of 3 divisible by 5 (3**5, 3**10, etc).

Powers of 7 also seem to have the issue, where math.log(7**5, 7) = 5.0000000001

According to https://docs.python.org/3/library/math.html#math.log :

"""With two arguments, return the logarithm of x to the given base, calculated as log(x)/log(base)."""

Indeed, this is consistent with that:

>>> from math import log
>>> log(243)
5.493061443340548
>>> log(3)
1.0986122886681098
>>> 5.493061443340548/1.0986122886681098
4.999999999999999

math.log is a floating-point operation, not an integer operation, and this is the nature of floating point operations: there can be rounding errors whenever there are intermediate steps of the computation. Strictly speaking, I would say this is not a bug, and results of floating point operations should generally be compared with math.isclose(), not with ==. Just like you can't expect (1/49)*49 == 1.0.

On the other hand, it may be theoretically possible to do some clever floating point tricks to get more accurate results out of math.log(a, b), like maybe some kind of correction factor involving the platform libm pow() function. I don't know if those are typically any more reliable than the status quo quotient-of-logs.

It turns out that the current implementation is already very good, within 1ulp in 99.85% of cases and perfect (as good as floats can get) in 65% of cases. So my thinking would be to leave the implementation as is.

#################################################

from decimal import Decimal as D, getcontext
from math import log, nextafter
from random import uniform

getcontext().prec = 200

N = 10_000
perfect = 0
good = 0

for i in range(N):
    x, base = uniform(1, 10**6), uniform(2, 100)
    d_result = float(D(x).ln() / D(base).ln())
    f_result = log(x, base)

    if d_result == f_result:
        perfect += 1
        good += 1
    elif nextafter(d_result, f_result) == f_result:
        good += 1

print(f"{perfect/N = }")
print(f"{good/N = }")

# results:
# perfect/N = 0.6503
# good/N = 0.9985

CPython's log() builds on the platform C's libm facilities, and C simply doesn't define primitives capable of producing a worst-case < 1 ulp error 2-argument log in reasonable time. Instead we have to build it out of two separate log operations, and a division. Each of those 3 operations can suffer its own rounding errors, which may, overall, cancel out, or build on each other. There's no error bound we can guarantee, although "< 2 ulp worst-case error" seems likely under modern libms.

Which is actually quite good. Doing better than that is out of scope for CPython's implementation. The easy way to get < 1 ulp is to use decimal to compute the intermediate results with excess precision. But that's also "too slow" for general use.

What Dennis did in his little test driver was fine for that, but we don't actually need to increase decimal's default precision at all to get < 1 ulp error in a converted-back-to-float-via-round-nearest result here.

Just an example (doesn't "prove" anything - but points at how to go about proving things):

>>> decimal.Decimal(3**8).ln() / decimal.Decimal(3).ln()
Decimal('7.999999999999999999999999999')
>>> float(_)
8.0