By OmnipotentEntity


2012-07-17 12:54:48 8 Comments

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 %.15f\n", a, sqrt((double)a), sqrt((double)a)*1e4, sqrt((double)a)*1e5, sqrt((double)a)*1e6);
  a = 999999999998;
  printf("%ld %.15f %.15f %.15f %.15f\n", a, sqrt((double)a), sqrt((double)a)*1e5, sqrt((double)a)*1e6, sqrt((double)a)*1e7);
  a = 99999999999998;
  printf("%ld %.15f %.15f %.15f %.15f\n", 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?

5 comments

@U10-Forward 2018-09-18 01:29:41

use decimal, (here is a clearer example):

>>> 2.3-2.2
0.09999999999999964
>>> from decimal import Decimal
>>> Decimal('2.3')-Decimal('2.2')
Decimal('0.1')
>>> float(Decimal('2.3')-Decimal('2.2'))
0.1
>>> 

@Ricardo Altamirano 2012-07-17 13:05:25

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.

@DSM 2012-07-17 13:08:13

+1 for mpmath. The problem with using Decimal numbers is that you can't do much in the way of math functions on Decimal objects, so if you're just playing around it's pretty limiting.

@Ricardo Altamirano 2012-07-17 13:10:40

@DSM I agree. I've only used mpmath for fairly simple problems, but nevertheless, I've found it a good, solid package.

@senderle 2012-07-17 13:30:58

Just to be clear -- I'm pretty sure that in your test of math.sqrt vs. Decimal.sqrt(), the result produced by math.sqrt is less correct, because of binary-to-decimal conversion. Consider the output of decimal.Decimal(math.sqrt(num) ** 2) * 7 vs. the output of decimal.Decimal(num.sqrt() ** 2) * 7.

@Ricardo Altamirano 2012-07-17 13:34:18

@senderle Good point. I didn't specify clearly enough which is more correct. +1 if you post that as an answer. (edit: or edit it into your existing answer, which just popped up). Otherwise I'm happy to have it added to my answer if you want.

@Ricardo Altamirano 2012-07-17 13:41:43

@senderle It's definitely a good point that probably belongs in one answer. I'd say add it to yours to be safe (otherwise as I said, I'm happy to link to your comment and edit it into mine). Any more caveats that I'm missing?

@OmnipotentEntity 2012-07-17 13:41:46

Considering that the actual value of sqrt(1/7) is 0.3779644730092272272145165362341800608157513118689214543383‌​33494171581260461469‌​08968005612663922051‌​5802... it seems that the decimal sqrt function is more accurate.

@Ricardo Altamirano 2012-07-17 13:43:36

@OmnipotentEntity I was just adding that to my answer, loosely based on the discussion above.

@DSM 2012-07-17 13:46:51

Rather than Decimal(math.sqrt(num)), you simply want num.sqrt(). Decimal(math.sqrt(num)) builds a Decimal object from a low-precision math function, rather than doing a high-precision sqrt.

@Ricardo Altamirano 2012-07-17 13:47:52

@DSM Exactly. I included both in my answer for comparison (as OmnipotentEntity pointed out, decimal is more accurate).

@Ricardo Altamirano 2012-07-17 13:49:02

@senderle No worries, my example is different. Plus it gives me an excuse to link to Wolfram Alpha!

@太極者無極而生 2016-01-08 09:50:56

hm... i don't think you can write it as "actual value" if the "actual value" is actually longer than that

@senderle 2012-07-17 13:24:35

For this particular problem, decimal is a great way to go, because it stores the decimal digits as tuples!

>>> a = decimal.Decimal(9999999998)
>>> a.as_tuple()
DecimalTuple(sign=0, digits=(9, 9, 9, 9, 9, 9, 9, 9, 9, 8), exponent=0)

Since you're looking for a property that is most naturally expressed in decimal notation, it's a bit silly to use a binary representation. The wikipedia page you linked to didn't indicate how many "non-grafting digits" may appear before the "grafting digits" begin, so this lets you specify:

>>> def isGrafting(dec, max_offset=5):
...     dec_digits = dec.as_tuple().digits
...     sqrt_digits = dec.sqrt().as_tuple().digits
...     windows = [sqrt_digits[o:o + len(dec_digits)] for o in range(max_offset)]
...     return dec_digits in windows
... 
>>> isGrafting(decimal.Decimal(9999999998))
True
>>> isGrafting(decimal.Decimal(77))
True

I think there's a good chance the result of Decimal.sqrt() will be more accurate, at least for this, than the result of math.sqrt() because of the conversion between binary representation and decimal representation. Consider the following, for example:

>>> num = decimal.Decimal(1) / decimal.Decimal(7)
>>> decimal.Decimal(math.sqrt(num) ** 2) * 7
Decimal('0.9999999999999997501998194593')
>>> decimal.Decimal(num.sqrt() ** 2) * 7
Decimal('1.000000000000000000000000000')

@f p 2012-07-17 12:57:30

You can try with Decimal instead of floatingpoint.

@Ned Batchelder 2012-07-17 13:02:49

Python has no built-in arbitrary-precision floats, but there are 3rd-party Python packages that use GMP: gmpy and PyGMP.

Related Questions

Sponsored Content

58 Answered Questions

[SOLVED] Calling an external command in Python

31 Answered Questions

[SOLVED] How do I check if a string is a number (float)?

23 Answered Questions

[SOLVED] How do I parse a string to a float or int?

21 Answered Questions

[SOLVED] Does Python have a ternary conditional operator?

25 Answered Questions

[SOLVED] How can I safely create a nested directory?

30 Answered Questions

[SOLVED] Is floating point math broken?

24 Answered Questions

[SOLVED] Limiting floats to two decimal points

16 Answered Questions

[SOLVED] What are metaclasses in Python?

10 Answered Questions

[SOLVED] Does Python have a string 'contains' substring method?

16 Answered Questions

[SOLVED] Difference between decimal, float and double in .NET?

Sponsored Content