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?

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,In Python 3.2.3, this outputs the first two lines

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.