Hi Anne, Long time no see ;) On Wed, Jun 5, 2013 at 10:07 AM, Anne Archibald <archibald@astron.nl> wrote:
Hi folks,
I recently came across an application I needed quad precision for (high-accuracy solution of a differential equation). I found a C++ library (odeint) that worked for the integration itself, but unfortunately it appears numpy is not able to work with quad precision arrays. For my application the quad precision is only really needed for integrator state, so I can manually convert my data to long doubles as I go to and from numpy, but it's a pain. So quad precision support would be nice.
There's a thread discussing quad support: http://mail.scipy.org/pipermail/numpy-discussion/2012-February/061080.html Essentially, there isn't any, but since gcc >= 4.6 supports them on Intel hardware (in software), it should be possible. (Then the thread got bogged down in bike-shedding about what to call them.)
What would it take to support quads in numpy? I looked into the numpy base dtype definitions, and it's a complex arrangement involving detection of platform support and templatized C code; in the end I couldn't really see where to start. But it seems to me all the basics are available: native C syntax for basic arithmetic, "qabs"-style versions of all the basic functions, NaNs and Infs. So how would one go about adding quad support?
There are some improvements for user types committed in 1.8-dev. Perhaps quad support could be added as a user type as it is still platform/compiler dependent. The rational type added to numpy could supply a template for adding the new type. Long term, we need to have quad support. If you run on SUN hardware I think it is already available as the extended precision type, but that doesn't help the majority of users at this point and I don't think LAPACK/BLAS supports it. Chuck