How to disable python rounding?

Question

I'm looking for a way to disable this:

print(0+1e-20) returns 1e-20, but print(1+1e-20) returns 1.0

I want it to return something like 1+1e-20.

I need it because of this problem:

from numpy import sqrt


def f1(x):
  return 1/((x+1)*sqrt(x))


def f2(x):
  return 1/((x+2)*sqrt(x+1))


def f3(x):
  return f2(x-1)


print(f1(1e-6))
print(f3(1e-6))
print(f1(1e-20))
print(f3(1e-20))

returns

999.9990000010001
999.998999986622
10000000000.0
main.py:10: RuntimeWarning: divide by zero encountered in double_scalars
  return 1/((x+2)*sqrt(x+1))
inf

f1 is the original function, f2 is f1 shifted by 1 to the left and f3 is f2 moved back by 1 to the right. By this logic, f1 and f3 should be equal to eachother, but this is not the case.

I know about decimal. Decimal and it doesn't work, because decimal doesn't support some functions including sin. If you could somehow make Decimal work for all functions, I'd like to know how.

This is the limit for double precision floating point numbers. It has only 15 to 16 significant digits, while your request exceeds 20. For the operation of single value, please use `math` instead of `numpy`. The functions in `math` support decimal. — Mechanic Pig, Oct 05 '22 at 14:03
@MechanicPig While `math.sin` can take a `decimal.Decimal` as an argument, the result is still `float`, so it may not be sufficient if one requires heightened precision to be maintained. OP, I have added a sample in `mpmath` to my answer. — Amadan, Oct 05 '22 at 14:34
note that you might want to read about the field of "numerical stability". for some easier reading: https://www.johndcook.com/blog/2010/06/07/math-library-functions-that-seem-unnecessary/ describes C math functions that also show up in Python's `math` module for the same reason as described in the blog post — Sam Mason, Oct 05 '22 at 22:25

Amadan · Answer 1 · 2022-10-05T14:37:06.973

Can't be done. There is no rounding - that would imply that the exact result ever existed, which it did not. Floating-point numbers have limited precision. The easiest way to envision this is an analogy with art and computer screens. Imagine someone making a fabulously detailed painting, and all you have is a 1024x768 screen to view it through. If a microscopic dot is added to the painting, the image on the screen might not change at all. Maybe you need a 4K screen instead.

In Python, the closest representable number after 1.0 is 1.0000000000000002 (*), according to math.nextafter(1.0, math.inf) (Python 3.9+ required for math.nextafter).1e-20 and 1 are too different in magnitude, so the result of their addition cannot be represented by a Python floating-point number, which are precise to up to about 16 digits.

See Is floating point math broken? for an in-depth explanation of the cause.

As this answer suggests, there are libraries like mpmath that implement arbitrary-precision arithmetics:

from mpmath import mp, mpf
mp.dps = 25 # set precision to 25 decimal digits
mp.sin(1)
# => mpf('0.8414709848078965066525023183')
mp.sin(1 + mpf('1e-20'))  # mpf is a constructor for mpmath floats
# => mpf('0.8414709848078965066579053457')

mpmath floats are sticky; if you add up an int and a mpf you get a mpf, so I did not have to write mp.mpf(1). The result is still not precise, but you can select what precision is sufficient for your needs. Also, note that the difference between these two results is, again, too small to be representable by Python's floating point numbers, so if the difference is meaningful to you, you have to keep it in the mpmath land:

float(mpf('0.8414709848078965066525023183')) == float(mpf('0.8414709848078965066579053457'))
# => True

(*) This is actually a white lie. The next number after 1.0 is 0x1.0000000000001, or 0b1.0000000000000000000000000000000000000000000000001, but Python doesn't like hexadecimal or binary float literals. 1.0000000000000002 is Python's approximation of that number for your decimal convenience.

`float.fromhex('0x1.0000000000001')` works... – Sam Mason Oct 05 '22 at 22:21 — Sam Mason, Oct 05 '22 at 22:21
@SamMason It does; but that's not a literal :) – Amadan Oct 06 '22 at 00:01 — Amadan, Oct 06 '22 at 00:01
oops, somehow missed that word! :( – Sam Mason Oct 06 '22 at 09:40 — Sam Mason, Oct 06 '22 at 09:40

score 1 · Accepted Answer · answered Oct 05 '22 at 22:39

As others have stated, in general this can't be done (due to how computers commonly represent numbers).

It's common to work with the precision you've got, ensuring that algorithms are numerically stable can be awkward.

In this case I'd redefine f1 to work on the logarithms of numbers, e.g.:

from numpy as sqrt, log, log1p, 

def f1(x):
    prod = log1p(x) + log(x) / 2
    return exp(-prod)

You might need to alter other parts of the code to work in log space as well depending on what you need to do. Note that most stats algorithms work with log-probabilities because it's much more compatible with how computers represent numbers.

f3 is a bit more work due to the subtraction.

How to disable python rounding?

2 Answers2