How to understand inaccurate results while using np.log?-CodePudding

I wish to calculate the natural logarithm of a value that is very close to 1, but not exactly one.

For example, np.log(1 1e-22) is 0 and not some non-zero value. However, np.log(1 1e-13) is not zero, and calculated to be 9.992007221625909e-14.

How can I understand this tradeoff of precision while using a numpy function vs. defining a numpy array with dtype as np.longdouble?

Floating precision information of numpy (v1.22.2) on the system I am using:

>>> np.finfo(np.longdouble)
finfo(resolution=1e-15, min=-1.7976931348623157e 308, max=1.7976931348623157e 308, dtype=float64)

>>> 1   np.finfo(np.longdouble).eps
1.0000000000000002

CodePudding user response：

For practical usage, take a look at np.log1p(x), which computes log(1 x) without the roundoff error that comes from 1 x. From the docs:

For real-valued input, log1p is accurate also for x so small that 1 x == 1 in floating-point accuracy.

Even the seemingly working example of log(1 1e-13) differs from the true value at the 3rd decimal place with 64-bit floats, and at the 6th with 128-bit floats (true value is from WolframAlpha):

>>> (1   1e-13) - 1
9.992007221626409e-14
>>> np.log(1   1e-13)
9.992007221625909e-14
>>> np.log(1   np.array(1e-13, dtype=np.float128))
9.999997791637127032e-14
>>> np.log1p(1e-13)
9.9999999999995e-14
>>> 9.999999999999500000000000033333333333330833333333333533333*10**-14
9.9999999999995e-14

CodePudding user response：

To complete the good solution of @yut23 using Numpy. If you need to deal with very small floats that does not fit in native types defined by Numpy or numbers close to 1 with a precision more than ~10 digits, then you can use the decimal package. It is slower than native float but it can gives you an arbitrary precision. The thing is it does not support the natural logarithm (ie. log) function directly since it is based on the transcendental Euler number (ie. e) which can hardly be computed with an arbitrary precision (at least not when the precision is huge). Fortunately, you can compute the natural logarithm from the based-10 logarithm and a precomputed Euler number based on existing number databases like this one (I guess 10000 digits should be enough for your needs ;) ). Here is an example:

import decimal
decimal.setcontext(decimal.Context(prec=100)) # 100 digits of precision
e = Decimal('2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516643')
result = (Decimal("1") Decimal("1e-13")).log10() / e.log10()

# result = 9.999999999999500000000000033333333333330833333333333533333333333316666666666668095238095237970238087E-14

The precision of result is of 99 digits (only the last one is not correct).