Hi, On Wed, Oct 12, 2011 at 11:24 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Tue, Oct 11, 2011 at 12:17 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
While struggling with floating point precision, I ran into this:
In [52]: a = 2**54+3
In [53]: a Out[53]: 18014398509481987L
In [54]: np.float128(a) Out[54]: 18014398509481988.0
In [55]: np.float128(a)-1 Out[55]: 18014398509481987.0
The line above tells us that float128 can exactly represent 2**54+3, but the line above that says that np.float128(2**54+3) rounds upwards as if it were a float64:
In [59]: np.float64(a) Out[59]: 18014398509481988.0
In [60]: np.float64(a)-1 Out[60]: 18014398509481988.0
Similarly:
In [66]: np.float128('1e308') Out[66]: 1.000000000000000011e+308
In [67]: np.float128('1e309') Out[67]: inf
Is it possible that float64 is being used somewhere in float128 casting?
The problem is probably in specifying the values. Python doesn't support long double and I expect python integers to be converted to doubles, then cast to long double.
Presumably our (numpy) casting function receives the python integer, and therefore its us who are doing the conversion? If so, surely that is a bug?
The only way to get around this is probably using string representations of the numbers, and I don't know how well/consistently numpy does that at the moment. If it calls python to do the job, then double is probably what is returned. It doesn't help on my system:
In [1]: float128("18014398509481987.0") Out[1]: 18014398509481988.0
Note though that:
In [66]: np.float128('1e308') Out[66]: 1.000000000000000011e+308
In [67]: np.float128('1e309') Out[67]: inf
and I infer that we (numpy) are using float64 for converting the strings; that seems to me to be the likely explanation of both phenomena - do you agree? See you, Matthew