Asked  7 Months ago    Answers:  5   Viewed   19 times

Just for fun and because it was really easy, I've written a short program to generate Grafting numbers, but because of floating point precision issues it's not finding some of the larger examples.

def isGrafting(a):
  for i in xrange(1, int(ceil(log10(a))) + 2):
    if a == floor((sqrt(a) * 10**(i-1)) % 10**int(ceil(log10(a)))):
      return 1

a = 0
while(1):
  if (isGrafting(a)):
    print "%d %.15f" % (a, sqrt(a))
  a += 1

This code misses at least one known Grafting number. 9999999998 => 99999.99998999999999949999999994999999999374999999912... It seems to drop extra precision after multiplying by 10**5.

>>> a = 9999999998
>>> sqrt(a)
99999.99999
>>> a == floor((sqrt(a) * 10**(5)) % 10**int(ceil(log10(a))))
False
>>> floor((sqrt(a) * 10**(5)) % 10**int(ceil(log10(a))))
9999999999.0
>>> print "%.15f" % sqrt(a)
99999.999989999996615
>>> print "%.15f" % (sqrt(a) * 10**5)
9999999999.000000000000000

So I wrote a short C++ program to see if it was my CPU truncating the floating point number or python somehow.

#include <cstdio>
#include <cmath>
#include <stdint.h>

int main()
{
  uint64_t a = 9999999998;
  printf("%ld %.15f %.15f %.15f %.15fn", a, sqrt((double)a), sqrt((double)a)*1e4, sqrt((double)a)*1e5, sqrt((double)a)*1e6);
  a = 999999999998;
  printf("%ld %.15f %.15f %.15f %.15fn", a, sqrt((double)a), sqrt((double)a)*1e5, sqrt((double)a)*1e6, sqrt((double)a)*1e7);
  a = 99999999999998;
  printf("%ld %.15f %.15f %.15f %.15fn", a, sqrt((double)a), sqrt((double)a)*1e6, sqrt((double)a)*1e7, sqrt((double)a)*1e8);
  return 0;
}

Which outputs:

9999999998 99999.999989999996615 999999999.899999976158142 9999999999.000000000000000 99999999990.000000000000000
999999999998 999999.999998999992386 99999999999.899993896484375 999999999999.000000000000000 9999999999990.000000000000000
99999999999998 9999999.999999899417162 9999999999999.900390625000000 99999999999999.000000000000000 999999999999990.000000000000000

So it looks like I'm running up hard against the limits of floating point precision and the CPU is chopping off the remaining bits because it thinks that the remaining difference is floating point error. Is there a way to work around this under Python? Or do I need to move to C and use GMP or something?

 Answers

19

In the standard library, the decimal module may be what you're looking for. Also, I have found mpmath to be quite helpful. The documentation has many great examples as well (unfortunately my office computer does not have mpmath installed; otherwise I would verify a few examples and post them).

One caveat about the decimal module, though. The module contains several in-built functions for simple mathematical operations (e.g. sqrt), but the results from these functions may not always match the corresponding function in math or other modules at higher precisions (although they may be more accurate). For example,

from decimal import *
import math

getcontext().prec = 30
num = Decimal(1) / Decimal(7)

print("   math.sqrt: {0}".format(Decimal(math.sqrt(num))))
print("decimal.sqrt: {0}".format(num.sqrt()))

In Python 3.2.3, this outputs the first two lines

   math.sqrt: 0.37796447300922719758631274089566431939601898193359375
decimal.sqrt: 0.377964473009227227214516536234
actual value: 0.3779644730092272272145165362341800608157513118689214

which as stated, isn't exactly what you would expect, and you can see that the higher the precision, the less the results match. Note that the decimal module does have more accuracy in this example, since it more closely matches the actual value.

Tuesday, June 1, 2021
 
edorian
answered 7 Months ago
90

str(0.47000000000000003) give '0.47' and float('0.47') can be 0.46999999999999997. This is due to the way floating point number are represented (see this wikipedia article)

Note: float(repr(0.47000000000000003)) or eval(repr(0.47000000000000003)) will give you the expected result, but you should use Decimal if you need precision.

Tuesday, June 1, 2021
 
CodeCaster
answered 7 Months ago
42

Check out ?0 (number) in Wikipedia

Basically IEEE does actually define a negative zero.

And by this definition for all purposes:

-0.0 == +0.0 == 0

I agree with aaronasterling that -0.0 and +0.0 are different objects. Making them equal (equality operator) makes sure that subtle bugs are not introduced in the code.
Think of a * b == c * d

>>> a = 3.4
>>> b =4.4
>>> c = -0.0
>>> d = +0.0
>>> a*c
-0.0
>>> b*d
0.0
>>> a*c == b*d
True
>>> 

[Edit: More info based on comments]

When I said for all practical purposes, I had chosen the word rather hastily. I meant standard equality comparison.

As the reference says, the IEEE standard defines comparison so that +0 = -0, rather than -0 < +0. Although it would be possible always to ignore the sign of zero, the IEEE standard does not do so. When a multiplication or division involves a signed zero, the usual sign rules apply in computing the sign of the answer.

Operations like divmod and atan2 exhibit this behavior. In fact, atan2 complies with the IEEE definition as does the underlying "C" lib.

>>> divmod(-0.0,100)
(-0.0, 0.0)
>>> divmod(+0.0,100)
(0.0, 0.0)

>>> math.atan2(0.0, 0.0) == math.atan2(-0.0, 0.0)
True 
>>> math.atan2(0.0, -0.0) == math.atan2(-0.0, -0.0)
False

One way is to find out through the documentation, if the implementation complies with IEEE behavior . It also seems from the discussion that there are subtle platform variations too.

However this aspect (IEEE definition compliance) has not been respected everywhere. See the rejection of PEP 754 due to disinterest! I am not sure if this was picked up later.

See also What Every Computer Scientist Should Know About Floating-Point Arithmetic.

Wednesday, June 2, 2021
 
Sidarta
answered 7 Months ago
86

The number 2.01 represented in binary is:

b10.00000010100011111100001010001111110000101000111111000010100011111100...

The computer uses only a finite number of digits to store floating-point values, but the binary representation of 2.01 requires infinitely many digits; as a result, it is rounded to the closest representable value:

b10.000000101000111111000010100011111100001010001111110

Expressed in decimal, this number is exactly:

2.0099999999999997868371792719699442386627197265625

When you print it out, it is rounded a second time to seventeen decimal digits, giving:

2.0099999999999998
Saturday, July 31, 2021
 
Abdel
answered 4 Months ago
11

There is the GWT-MATH library at http://code.google.com/p/gwt-math/.

However, I warn you, it's a GWT jsni overlay of a java->javascript automated conversion of java.BigDecimal (actually the old com.ibm.math.BigDecimal).

It works, but speedy it is not. (Nor lean. It will pad on a good 70k into your project).

At my workplace, we are working on a fixed point simple decimal, but nothing worth releasing yet. :(

Thursday, August 5, 2021
 
Gersom
answered 4 Months ago
Only authorized users can answer the question. Please sign in first, or register a free account.
Not the answer you're looking for? Browse other questions tagged :  
Share