tag:blogger.com,1999:blog-6568744196982634289.post7881611169286097799..comments2012-09-06T19:25:36.717-07:00Comments on Ondřej Čertík: When double precision is not enoughOndřej Čertíkhttps://plus.google.com/104039945248245758823noreply@blogger.comBlogger8125tag:blogger.com,1999:blog-6568744196982634289.post-17074024394295632152012-03-07T02:47:34.073-08:002012-03-07T02:47:34.073-08:00Pavel, thanks a lot for the feedback and links. In...Pavel, thanks a lot for the feedback and links. Indeed, I think that the correct C++ approach is to use templates and classes for scalar types and also arrays. So compared to Fortran, it is more heavyweight, but it does the job and also it is more general (i.e. as you wrote, you can then easily use MPFR --- while in Fortran, few more modifications would be needed).Ondřej Čertíkhttp://www.blogger.com/profile/02845032202161204018noreply@blogger.comtag:blogger.com,1999:blog-6568744196982634289.post-31734358825274500022012-03-06T23:37:29.009-08:002012-03-06T23:37:29.009-08:00Properly designed C++ mathematical code would use ...Properly designed C++ mathematical code would use templates to be invariant from scalar types (along with type traits and other meta-magic).<br /><br />Then change of scalar type (float-double-whatever) wouldn't be difficult operation.<br /><br />Thanks to this approach I was able to convert <a href="http://eigen.tuxfamily.org/index.php?title=Main_Page" rel="nofollow">Eigen: C++ matrix library</a> to arbitrary precision using my custom made multiprecision floating-point scalar type - <a href="http://www.holoborodko.com/pavel/mpfr/" rel="nofollow">mpreal</a> (C++ wrapper around MPFR/MPIR).<br /><br />I use similar OOP techniques to enable Matlab to run computations in arbitrary precision without changes in source code:<br /><a href="http://www.advanpix.com/" rel="nofollow">Multiprecision Computing Toolbox for MATLAB</a>Pavel Holoborodkohttp://www.blogger.com/profile/07669265710142706893noreply@blogger.comtag:blogger.com,1999:blog-6568744196982634289.post-57662378583363524862012-01-27T11:43:33.927-08:002012-01-27T11:43:33.927-08:00CarlK --- I didn't try his double-double and q...CarlK --- I didn't try his double-double and quad-double, but I tried his MPFUN90, and compared it with MPFR. On the sin(pi/10)**2 + cos(pi/10)**2 example, MPFUN90 takes about 6s on my computer, and MPFR is *instant*. I used the same precision for both, I forgot which, I think it was 10000, maybe it was too much.Ondřej Čertíkhttp://www.blogger.com/profile/02845032202161204018noreply@blogger.comtag:blogger.com,1999:blog-6568744196982634289.post-38703746151133186952012-01-27T11:39:25.707-08:002012-01-27T11:39:25.707-08:00sojrhthu, that is a great idea. But I don't kn...sojrhthu, that is a great idea. But I don't know how to force LAPACK to use it. For my own code that will work well.<br /><br />Markus, thanks for he tip. This looks really cool. And I guess you can also hook in MPFR automatically to get any precision. I like Fortran, because it makes things easy for scientific computing, but I also like C++, because almost anything is possible to do, if one really wants to.Ondřej Čertíkhttp://www.blogger.com/profile/02845032202161204018noreply@blogger.comtag:blogger.com,1999:blog-6568744196982634289.post-85115211257109336082012-01-27T06:00:26.621-08:002012-01-27T06:00:26.621-08:00Does anyone knows how David H. Bailey's double...Does anyone knows how David H. Bailey's double-double and quad-double package (based on IEEE double precision pairs) compares to libquadmath?. See: http://crd-legacy.lbl.gov/~dhbailey/mpdist/CarlKhttp://www.blogger.com/profile/00917295480785698315noreply@blogger.comtag:blogger.com,1999:blog-6568744196982634289.post-34998295641103490672012-01-27T05:53:53.986-08:002012-01-27T05:53:53.986-08:00On x86/amd64, gfortran also supports the 80-bit &q...On x86/amd64, gfortran also supports the 80-bit "extended precision" of the processor with real(kind=10). It's a little bit slower than the standard double precision real(kind=8), but much faster than software simulated quadruple precision real(kind=16). So if 80 bit are enough in your case, you might want to try this.sojrhthuhttp://www.blogger.com/profile/12110257998930253503noreply@blogger.comtag:blogger.com,1999:blog-6568744196982634289.post-72106357096533212132012-01-27T01:13:26.321-08:002012-01-27T01:13:26.321-08:00I suppose the cleanest solution is to use a typede...I suppose the cleanest solution is to use a typedef. If you want to get into trouble, you may use #define double __float128.<br /><br />The C++-ish way to define something similar to Fortran's real type would be a simple class template. But it feels a lot more clumsy.<br /><br />You can use static assertions to get better error messages in case something goes wrong and other techniques to improve compatibility. Furthermore, you can use preprocessor macros to work around different compilers, implementations, etc.<br /><br />I don't know if this has been done already.<br /><br />typedef typename real< FLOAT_BYTE_SIZE >::t my_real;<br /><br />template < int byte_size ><br />struct real;<br />template <> struct real<4><br />{<br /> typedef float t;<br />};<br />template <> struct real<8><br />{<br /> typedef double t;<br />};<br />template <> struct real<10><br />{<br /> STATIC_ASSERT(sizeof(long double) == 10);<br /> typedef long double t;<br />};<br />#ifdef MACRO_TO_CHECK_WHETHER_YOU_ARE_USING_A_GNU_COMPILER_ETC<br />template <> struct real<16><br />{<br /> typedef __float128 t;<br />};<br />#endifMarkushttp://www.blogger.com/profile/10750652068584028947noreply@blogger.comtag:blogger.com,1999:blog-6568744196982634289.post-57768150197070166982012-01-26T18:58:32.677-08:002012-01-26T18:58:32.677-08:00Could you be sneaky and #define double __float128?...Could you be sneaky and #define double __float128?Marius Gedminashttp://www.blogger.com/profile/15155998626202067226noreply@blogger.com