float128 / longdouble on PPC - is it broken?
Hi, I just ran into this on a PPC machine: In [1]: import numpy as np In [2]: np.__version__ Out[2]: '2.0.0.dev-4daf949' In [3]: res = np.longdouble(2)**64 In [4]: res Out[4]: 18446744073709551616.0 In [5]: 2**64 Out[5]: 18446744073709551616L In [6]: res-1 Out[6]: 36893488147419103231.0 Same for numpy 1.4.1. I don't have a SPARC to test on but I believe it's the same double-double type? See you, Matthew
25.10.2011 06:59, Matthew Brett kirjoitti:
res = np.longdouble(2)**64 res-1 36893488147419103231.0
Can you check if long double works properly (not a given) in C on that platform: long double x; x = powl(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x); or, in case the platform doesn't have powl: long double x; x = pow(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x);
Hi, On Tue, Oct 25, 2011 at 2:43 AM, Pauli Virtanen <pav@iki.fi> wrote:
25.10.2011 06:59, Matthew Brett kirjoitti:
res = np.longdouble(2)**64 res-1 36893488147419103231.0
Can you check if long double works properly (not a given) in C on that platform:
long double x; x = powl(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x);
or, in case the platform doesn't have powl:
long double x; x = pow(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x);
Both the same as numpy: [mb312@jerry ~]$ gcc test.c test.c: In function 'main': test.c:5: warning: incompatible implicit declaration of built-in function 'powl' [mb312@jerry ~]$ ./a.out 1.84467e+19 3.68935e+19 Thanks, Matthew
On Tue, Oct 25, 2011 at 11:45 AM, Matthew Brett <matthew.brett@gmail.com>wrote:
Hi,
On Tue, Oct 25, 2011 at 2:43 AM, Pauli Virtanen <pav@iki.fi> wrote:
25.10.2011 06:59, Matthew Brett kirjoitti:
res = np.longdouble(2)**64 res-1 36893488147419103231.0
Can you check if long double works properly (not a given) in C on that platform:
long double x; x = powl(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x);
or, in case the platform doesn't have powl:
long double x; x = pow(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x);
Both the same as numpy:
[mb312@jerry ~]$ gcc test.c test.c: In function 'main': test.c:5: warning: incompatible implicit declaration of built-in function 'powl'
I think implicit here means that that the arguments and the return values are treated as integers. Did you #include <math.h>?
[mb312@jerry ~]$ ./a.out 1.84467e+19 3.68935e+19
Chuck
Hi, On Tue, Oct 25, 2011 at 10:52 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Tue, Oct 25, 2011 at 11:45 AM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Tue, Oct 25, 2011 at 2:43 AM, Pauli Virtanen <pav@iki.fi> wrote:
25.10.2011 06:59, Matthew Brett kirjoitti:
res = np.longdouble(2)**64 res-1 36893488147419103231.0
Can you check if long double works properly (not a given) in C on that platform:
long double x; x = powl(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x);
or, in case the platform doesn't have powl:
long double x; x = pow(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x);
Both the same as numpy:
[mb312@jerry ~]$ gcc test.c test.c: In function 'main': test.c:5: warning: incompatible implicit declaration of built-in function 'powl'
I think implicit here means that that the arguments and the return values are treated as integers. Did you #include <math.h>?
Ah - you've detected my severe ignorance of c. But with math.h, the result is the same, #include <stdio.h> #include <math.h> int main(int argc, char* argv) { long double x; x = pow(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x); } See you, Matthew
On Tue, Oct 25, 2011 at 11:05 AM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Tue, Oct 25, 2011 at 10:52 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Tue, Oct 25, 2011 at 11:45 AM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Tue, Oct 25, 2011 at 2:43 AM, Pauli Virtanen <pav@iki.fi> wrote:
25.10.2011 06:59, Matthew Brett kirjoitti:
res = np.longdouble(2)**64 res-1 36893488147419103231.0
Can you check if long double works properly (not a given) in C on that platform:
long double x; x = powl(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x);
or, in case the platform doesn't have powl:
long double x; x = pow(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x);
Both the same as numpy:
[mb312@jerry ~]$ gcc test.c test.c: In function 'main': test.c:5: warning: incompatible implicit declaration of built-in function 'powl'
I think implicit here means that that the arguments and the return values are treated as integers. Did you #include <math.h>?
Ah - you've detected my severe ignorance of c. But with math.h, the result is the same,
#include <stdio.h> #include <math.h>
int main(int argc, char* argv) { long double x; x = pow(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x); }
By the way - if you want a login to this machine, let me know - it's always on and we're using it as a buildslave already. Matthew
On 25 Oct 2011, at 20:05, Matthew Brett wrote:
Both the same as numpy:
[mb312@jerry ~]$ gcc test.c test.c: In function 'main': test.c:5: warning: incompatible implicit declaration of built-in function 'powl'
I think implicit here means that that the arguments and the return values are treated as integers. Did you #include <math.h>?
Ah - you've detected my severe ignorance of c. But with math.h, the result is the same,
#include <stdio.h> #include <math.h>
int main(int argc, char* argv) { long double x; x = pow(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x); }
What system/compiler is this? I am getting ./ldouble 1.84467e+19 1.84467e+19 and
res = np.longdouble(2)**64 res 18446744073709551616.0 2**64 18446744073709551616L res-1 18446744073709551615.0 np.__version__ '1.6.1'
as well as with
np.__version__ '2.0.0.dev-3d06f02' [yes, not very up to date]
and for all gcc versions /usr/bin/gcc -v Using built-in specs. Target: powerpc-apple-darwin9 Configured with: /var/tmp/gcc/gcc-5493~1/src/configure --disable-checking -enable-werror --prefix=/usr --mandir=/share/man --enable-languages=c,objc,c++,obj-c++ --program-transform-name=/^[cg][^.-]*$/s/$/-4.0/ --with-gxx-include-dir=/include/c++/4.0.0 --with-slibdir=/usr/lib --build=i686-apple-darwin9 --program-prefix= --host=powerpc-apple-darwin9 --target=powerpc-apple-darwin9 Thread model: posix gcc version 4.0.1 (Apple Inc. build 5493) to /sw/bin/gcc-fsf-4.6 -v Using built-in specs. COLLECT_GCC=/sw/bin/gcc-fsf-4.6 COLLECT_LTO_WRAPPER=/sw/lib/gcc4.6/libexec/gcc/powerpc-apple-darwin9.8.0/4.6.1/lto-wrapper Target: powerpc-apple-darwin9.8.0 Configured with: ../gcc-4.6.1/configure --prefix=/sw --prefix=/sw/lib/gcc4.6 --mandir=/sw/share/man --infodir=/sw/lib/gcc4.6/info --enable-languages=c,c++,fortran,lto,objc,obj-c++,java --with-gmp=/sw --with-libiconv-prefix=/sw --with-ppl=/sw --with-cloog=/sw --with-mpc=/sw --with-system-zlib --x-includes=/usr/X11R6/include --x-libraries=/usr/X11R6/lib --program-suffix=-fsf-4.6 --enable-cloog-backend=isl --disable-libjava-multilib --disable-libquadmath Thread model: posix gcc version 4.6.1 (GCC) uname -a Darwin osiris.astro.physik.uni-goettingen.de 9.8.0 Darwin Kernel Version 9.8.0: Wed Jul 15 16:57:01 PDT 2009; root:xnu-1228.15.4~1/RELEASE_PPC Power Macintosh Cheers, Derek
Hi, On Tue, Oct 25, 2011 at 12:01 PM, Derek Homeier <derek@astro.physik.uni-goettingen.de> wrote:
On 25 Oct 2011, at 20:05, Matthew Brett wrote:
Both the same as numpy:
[mb312@jerry ~]$ gcc test.c test.c: In function 'main': test.c:5: warning: incompatible implicit declaration of built-in function 'powl'
I think implicit here means that that the arguments and the return values are treated as integers. Did you #include <math.h>?
Ah - you've detected my severe ignorance of c. But with math.h, the result is the same,
#include <stdio.h> #include <math.h>
int main(int argc, char* argv) { long double x; x = pow(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x); }
What system/compiler is this? I am getting ./ldouble 1.84467e+19 1.84467e+19
and
res = np.longdouble(2)**64 res 18446744073709551616.0 2**64 18446744073709551616L res-1 18446744073709551615.0 np.__version__ '1.6.1'
as well as with
np.__version__ '2.0.0.dev-3d06f02' [yes, not very up to date]
and for all gcc versions /usr/bin/gcc -v Using built-in specs. Target: powerpc-apple-darwin9 Configured with: /var/tmp/gcc/gcc-5493~1/src/configure --disable-checking -enable-werror --prefix=/usr --mandir=/share/man --enable-languages=c,objc,c++,obj-c++ --program-transform-name=/^[cg][^.-]*$/s/$/-4.0/ --with-gxx-include-dir=/include/c++/4.0.0 --with-slibdir=/usr/lib --build=i686-apple-darwin9 --program-prefix= --host=powerpc-apple-darwin9 --target=powerpc-apple-darwin9 Thread model: posix gcc version 4.0.1 (Apple Inc. build 5493)
to
/sw/bin/gcc-fsf-4.6 -v Using built-in specs. COLLECT_GCC=/sw/bin/gcc-fsf-4.6 COLLECT_LTO_WRAPPER=/sw/lib/gcc4.6/libexec/gcc/powerpc-apple-darwin9.8.0/4.6.1/lto-wrapper Target: powerpc-apple-darwin9.8.0 Configured with: ../gcc-4.6.1/configure --prefix=/sw --prefix=/sw/lib/gcc4.6 --mandir=/sw/share/man --infodir=/sw/lib/gcc4.6/info --enable-languages=c,c++,fortran,lto,objc,obj-c++,java --with-gmp=/sw --with-libiconv-prefix=/sw --with-ppl=/sw --with-cloog=/sw --with-mpc=/sw --with-system-zlib --x-includes=/usr/X11R6/include --x-libraries=/usr/X11R6/lib --program-suffix=-fsf-4.6 --enable-cloog-backend=isl --disable-libjava-multilib --disable-libquadmath Thread model: posix gcc version 4.6.1 (GCC)
uname -a Darwin osiris.astro.physik.uni-goettingen.de 9.8.0 Darwin Kernel Version 9.8.0: Wed Jul 15 16:57:01 PDT 2009; root:xnu-1228.15.4~1/RELEASE_PPC Power Macintosh
mb312@jerry ~]$ gcc -v Using built-in specs. Target: powerpc-apple-darwin8 Configured with: /var/tmp/gcc/gcc-5370~2/src/configure --disable-checking -enable-werror --prefix=/usr --mandir=/share/man --enable-languages=c,objc,c++,obj-c++ --program-transform-name=/^[cg][^.-]*$/s/$/-4.0/ --with-gxx-include-dir=/include/c++/4.0.0 --with-slibdir=/usr/lib --build=powerpc-apple-darwin8 --host=powerpc-apple-darwin8 --target=powerpc-apple-darwin8 Thread model: posix gcc version 4.0.1 (Apple Computer, Inc. build 5370) [mb312@jerry ~]$ uname -a Darwin jerry.bic.berkeley.edu 8.11.0 Darwin Kernel Version 8.11.0: Wed Oct 10 18:26:00 PDT 2007; root:xnu-792.24.17~1/RELEASE_PPC Power Macintosh powerpc Best, Matthew
25.10.2011 19:45, Matthew Brett kirjoitti: [clip]
or, in case the platform doesn't have powl:
long double x; x = pow(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x);
Both the same as numpy:
[mb312@jerry ~]$ gcc test.c test.c: In function 'main': test.c:5: warning: incompatible implicit declaration of built-in function 'powl' [mb312@jerry ~]$ ./a.out 1.84467e+19 3.68935e+19
This result may indicate that it's the *printing* of long doubles that's broken. Note how the value cast as double prints the correct result, whereas the %Lg format code gives something wrong. Can you try to check this by doing something like: - do some set of calculations using np.longdouble in Numpy (that requires the extra accuracy) - at the end, cast the result back to double -- Pauli Virtanen
Hi, On Tue, Oct 25, 2011 at 11:14 AM, Pauli Virtanen <pav@iki.fi> wrote:
25.10.2011 19:45, Matthew Brett kirjoitti: [clip]
or, in case the platform doesn't have powl:
long double x; x = pow(2, 64); x -= 1; printf("%g %Lg\n", (double)x, x);
Both the same as numpy:
[mb312@jerry ~]$ gcc test.c test.c: In function 'main': test.c:5: warning: incompatible implicit declaration of built-in function 'powl' [mb312@jerry ~]$ ./a.out 1.84467e+19 3.68935e+19
This result may indicate that it's the *printing* of long doubles that's broken. Note how the value cast as double prints the correct result, whereas the %Lg format code gives something wrong.
Ah - sorry - I see now what you were trying to do.
Can you try to check this by doing something like:
- do some set of calculations using np.longdouble in Numpy (that requires the extra accuracy)
- at the end, cast the result back to double
In [1]: import numpy as np In [2]: res = np.longdouble(2)**64 In [6]: res / 2**32 Out[6]: 4294967296.0 In [7]: (res-1) / 2**32 Out[7]: 8589934591.9999999998 In [8]: np.float((res-1) / 2**32) Out[8]: 4294967296.0 In [9]: np.float((res) / 2**32) Out[9]: 4294967296.0 Thanks, Matthew
25.10.2011 20:29, Matthew Brett kirjoitti: [clip]
In [7]: (res-1) / 2**32 Out[7]: 8589934591.9999999998
In [8]: np.float((res-1) / 2**32) Out[8]: 4294967296.0
Looks like a bug in the C library installed on the machine, then. It's either in wontfix territory for us, or in the "cast to doubles before formatting" one. In the latter case, one would have to maintain a list of broken C libraries (ugh). -- Pauli Virtanen
Hi, On 25 Oct 2011, at 21:14, Pauli Virtanen wrote:
25.10.2011 20:29, Matthew Brett kirjoitti: [clip]
In [7]: (res-1) / 2**32 Out[7]: 8589934591.9999999998
In [8]: np.float((res-1) / 2**32) Out[8]: 4294967296.0
Looks like a bug in the C library installed on the machine, then.
It's either in wontfix territory for us, or in the "cast to doubles before formatting" one. In the latter case, one would have to maintain a list of broken C libraries (ugh).
as it appears to be a "Tiger-only" problem, probably the former? On 25 Oct 2011, at 21:13, Matthew Brett wrote:
[mb312@jerry ~]$ uname -a Darwin jerry.bic.berkeley.edu 8.11.0 Darwin Kernel Version 8.11.0: Wed Oct 10 18:26:00 PDT 2007; root:xnu-792.24.17~1/RELEASE_PPC Power Macintosh powerpc
Cheers, Derek
Hi, On Tue, Oct 25, 2011 at 12:14 PM, Pauli Virtanen <pav@iki.fi> wrote:
25.10.2011 20:29, Matthew Brett kirjoitti: [clip]
In [7]: (res-1) / 2**32 Out[7]: 8589934591.9999999998
In [8]: np.float((res-1) / 2**32) Out[8]: 4294967296.0
Looks like a bug in the C library installed on the machine, then.
It's either in wontfix territory for us, or in the "cast to doubles before formatting" one. In the latter case, one would have to maintain a list of broken C libraries (ugh).
How about a check at import time and a warning when printing? Is that hard to do? See you, Matthew
On Tue, Oct 25, 2011 at 8:22 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Tue, Oct 25, 2011 at 12:14 PM, Pauli Virtanen <pav@iki.fi> wrote:
25.10.2011 20:29, Matthew Brett kirjoitti: [clip]
In [7]: (res-1) / 2**32 Out[7]: 8589934591.9999999998
In [8]: np.float((res-1) / 2**32) Out[8]: 4294967296.0
Looks like a bug in the C library installed on the machine, then.
It's either in wontfix territory for us, or in the "cast to doubles before formatting" one. In the latter case, one would have to maintain a list of broken C libraries (ugh).
How about a check at import time and a warning when printing? Is that hard to do?
That's fragile IMO. I think that Chuck summed it well: long double are not portable, don't use them unless you have to or you can rely on platform-specificities. I would rather spend some time on implementing/integrating portable quad precision in software, cheers, David
Hi, On Tue, Oct 25, 2011 at 2:58 PM, David Cournapeau <cournape@gmail.com> wrote:
On Tue, Oct 25, 2011 at 8:22 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Tue, Oct 25, 2011 at 12:14 PM, Pauli Virtanen <pav@iki.fi> wrote:
25.10.2011 20:29, Matthew Brett kirjoitti: [clip]
In [7]: (res-1) / 2**32 Out[7]: 8589934591.9999999998
In [8]: np.float((res-1) / 2**32) Out[8]: 4294967296.0
Looks like a bug in the C library installed on the machine, then.
It's either in wontfix territory for us, or in the "cast to doubles before formatting" one. In the latter case, one would have to maintain a list of broken C libraries (ugh).
How about a check at import time and a warning when printing? Is that hard to do?
That's fragile IMO. I think that Chuck summed it well: long double are not portable, don't use them unless you have to or you can rely on platform-specificities.
That reminds me of the old joke about the Irishman giving directions - "If I were you, I wouldn't start from here".
I would rather spend some time on implementing/integrating portable quad precision in software,
I guess from your answer that such a warning would be complicated to implement, and if that's the case, I can imagine it would be low priority. See you, Matthew
On Wed, Oct 26, 2011 at 12:49 AM, Matthew Brett <matthew.brett@gmail.com> wrote:
That reminds me of the old joke about the Irishman giving directions - "If I were you, I wouldn't start from here".
Sounds about accurate `1
I would rather spend some time on implementing/integrating portable quad precision in software,
I guess from your answer that such a warning would be complicated to implement, and if that's the case, I can imagine it would be low priority.
No, it would not be hard to implement, but I don't think it is a nice solution, because it is very platform dependent, so it will likely bitrot as it won't appear on usual platforms. David
On Tue, Oct 25, 2011 at 4:49 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
I guess from your answer that such a warning would be complicated to implement, and if that's the case, I can imagine it would be low priority.
I assume the problem is more that it would be a weirdo check that becomes a maintenance burden ("what is this doing here? Do we still need it? who knows?") than that it would be hard to do. You can easily do it yourself as a workaround... if not str(np.longdouble(2)**64 - 1).startswith("1844"): warn("Printing of longdoubles is fubared! Beware! Beware!") -- Nathaniel
Hi, On Wed, Oct 26, 2011 at 1:07 AM, Nathaniel Smith <njs@pobox.com> wrote:
On Tue, Oct 25, 2011 at 4:49 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
I guess from your answer that such a warning would be complicated to implement, and if that's the case, I can imagine it would be low priority.
I assume the problem is more that it would be a weirdo check that becomes a maintenance burden ("what is this doing here? Do we still need it? who knows?") than that it would be hard to do.
You can easily do it yourself as a workaround...
if not str(np.longdouble(2)**64 - 1).startswith("1844"): warn("Printing of longdoubles is fubared! Beware! Beware!")
Thanks - yes - I was only thinking of someone like me getting confused and thinking badly of us if they run into this. See you, Matthew
26.10.2011 10:07, Nathaniel Smith kirjoitti:
On Tue, Oct 25, 2011 at 4:49 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
I guess from your answer that such a warning would be complicated to implement, and if that's the case, I can imagine it would be low priority.
I assume the problem is more that it would be a weirdo check that becomes a maintenance burden ("what is this doing here? Do we still need it? who knows?") than that it would be hard to do.
This check should be done compile-time, and set a flag like BROKEN_LONG_DOUBLE_FORMATTING for compilation. IIRC, we already have some related workarounds in the code (for Windows). The main argument is, (i) it's not our bug, (ii) adding workarounds for broken platforms should be weighed having the commonness of the issue in mind. But I'll admit that someone could have got into work and implemented the workaround in the same time it has taken to write mails to this thread :) Pauli
On Mon, Oct 24, 2011 at 10:59 PM, Matthew Brett <matthew.brett@gmail.com>wrote:
Hi,
I just ran into this on a PPC machine:
In [1]: import numpy as np
In [2]: np.__version__ Out[2]: '2.0.0.dev-4daf949'
In [3]: res = np.longdouble(2)**64
In [4]: res Out[4]: 18446744073709551616.0
In [5]: 2**64 Out[5]: 18446744073709551616L
In [6]: res-1 Out[6]: 36893488147419103231.0
Same for numpy 1.4.1.
I don't have a SPARC to test on but I believe it's the same double-double type?
The PPC uses two doubles to represent long doubles, the SPARC uses software emulation of ieee quad precision for long doubles, very different. The subtraction of 1 working like multiplication by two is strange, perhaps the one is getting subtracted from the exponent somehow? It would be interesting to see if the same problem happens in pure c. As a work around, can I ask what you are trying to do with the long doubles? Chuck
Hi, On Tue, Oct 25, 2011 at 7:31 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Oct 24, 2011 at 10:59 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
I just ran into this on a PPC machine:
In [1]: import numpy as np
In [2]: np.__version__ Out[2]: '2.0.0.dev-4daf949'
In [3]: res = np.longdouble(2)**64
In [4]: res Out[4]: 18446744073709551616.0
In [5]: 2**64 Out[5]: 18446744073709551616L
In [6]: res-1 Out[6]: 36893488147419103231.0
Same for numpy 1.4.1.
I don't have a SPARC to test on but I believe it's the same double-double type?
The PPC uses two doubles to represent long doubles, the SPARC uses software emulation of ieee quad precision for long doubles, very different.
Yes, thanks - I read more after my post. I guess from this: http://publib.boulder.ibm.com/infocenter/aix/v6r1/index.jsp?topic=/com.ibm.a... that AIX does use double-double.
The subtraction of 1 working like multiplication by two is strange, perhaps the one is getting subtracted from the exponent somehow? It would be interesting to see if the same problem happens in pure c.
As a work around, can I ask what you are trying to do with the long doubles?
I was trying to use them as an intermediate format for high-precision floating point calculations, before converting to integers. Best, Matthew
participants (6)
-
Charles R Harris -
David Cournapeau -
Derek Homeier -
Matthew Brett -
Nathaniel Smith -
Pauli Virtanen