[Ann] Proposed Math library additions (complex number inverse trig functions).
As part of the project to produce a Boost TR1 implementation I've implemented the additional algorithms in <complex> that are part of the TR: fabs asin asinh acos acosh atan atanh It occurs to me that maybe these should really be spun out as a separate library rather than "slipped in" as part of another submission (comments welcome on that). In the mean time these algorithms can be downloaded as part of the TR1 library here: http://boost-sandbox.sourceforge.net/vault/index.php?&direction=0&order=&directory=tr1 All the inverse trig functions are defined for the full complex plain, and do the right thing even for exceptional input (no overflow or underflow during the computation if the result is actually representable). They also do the right thing for infinities and nan's, as specified by C99. The asin/asinh/acos/acosh algorithms have a theoretical max error bound of 9.5e (experimental bounds are much lower, about 4e), there is no theoretical work on atan/atanh that I know of, but experimentally the error bound seems lower than for asin/acos. This is some of the most heavily commented code I've written, so refer to the header for more information on the literature sources used etc. During the implementation I found that I needed a few other things that might make useful Boost components (part implemented, part wish list as noted below): log1p: ~~~~ This algorithm is part of C99, but by no means all compilers support it yet, I used a Taylor series expansion for small x: I'm aware that there are much more efficient methods, but optimizing compilers completely trash the logic of these (Intel C++ proved to be particularly bad). Series Summation (Kahan Method): ~~~~~~~~~~~~~~~~~~~~~~~~ I wrote a small routine to apply the Kahan summation method to a nullary-functor that generates successive terms in a series. Summation stops when the next term is below a particular threshold. I've found the Kahan method to be significantly more accurate for certain series, but again some highly optimizing compilers can cause small errors to creep in (I believe that the issue is double-rounding of intermediate terms on machines with extended-double registers). Probably my current interface to this one needs some more thought before it could be anything other than a detail, and if we're talking about wish lists, infinite product, and continued fraction evaluation would be useful additions too.... Floating point testing: ~~~~~~~~~~~~~~ Currently just a version of C99's isnan, and only works for quiet nan's with IEEE-conforming compilers. More of a wish list really, I could really use all the C99 fp-testing macros (as C++ functions obviously). Numeric constants: ~~~~~~~~~~~~~ Currently I've embedded the numeric constants used inside the implementation, a numeric constants library really would have simplified things... just a pity we could never agree on an interface last time this came up. Whew, I think that's all, thanks for reading this far! Regards, John.
Somewhere in the E.U., le 05/04/2005 Bonjour In article <007601c537a0$696a8210$51eb1b52@fuji>, "John Maddock" <john@johnmaddock.co.uk> wrote:
As part of the project to produce a Boost TR1 implementation I've implemented the additional algorithms in <complex> that are part of the TR:
fabs asin asinh acos acosh atan atanh
It occurs to me that maybe these should really be spun out as a separate library rather than "slipped in" as part of another submission (comments welcome on that).
FWIW, asinh, acosh and atanh are already in Boost! If you wish, I can retire them (though I planed to add valarray support, I even printed a few days ago the discussion which had taken place on this two years (ahem) ago...).
In the mean time these algorithms can be downloaded as part of the TR1 library here: http://boost-sandbox.sourceforge.net/vault/index.php?&direction=0&order=&direc tory=tr1
(not showing up yet...)
All the inverse trig functions are defined for the full complex plain, and do the right thing even for exceptional input (no overflow or underflow during the computation if the result is actually representable). They also do the right thing for infinities and nan's, as specified by C99.
The asin/asinh/acos/acosh algorithms have a theoretical max error bound of 9.5e (experimental bounds are much lower, about 4e), there is no theoretical work on atan/atanh that I know of, but experimentally the error bound seems lower than for asin/acos. This is some of the most heavily commented code I've written, so refer to the header for more information on the literature sources used etc.
During the implementation I found that I needed a few other things that might make useful Boost components (part implemented, part wish list as noted below):
log1p: ~~~~
This algorithm is part of C99, but by no means all compilers support it yet, I used a Taylor series expansion for small x: I'm aware that there are much more efficient methods, but optimizing compilers completely trash the logic of these (Intel C++ proved to be particularly bad).
I had planed to work on it (and a few related such as expm1...) something liketwo years ago and had planed to add it to the same special functions library. They really are important.
Series Summation (Kahan Method): ~~~~~~~~~~~~~~~~~~~~~~~~
I wrote a small routine to apply the Kahan summation method to a nullary-functor that generates successive terms in a series. Summation stops when the next term is below a particular threshold. I've found the Kahan method to be significantly more accurate for certain series, but again some highly optimizing compilers can cause small errors to creep in (I believe that the issue is double-rounding of intermediate terms on machines with extended-double registers). Probably my current interface to this one needs some more thought before it could be anything other than a detail, and if we're talking about wish lists, infinite product, and continued fraction evaluation would be useful additions too....
Floating point testing: ~~~~~~~~~~~~~~
Currently just a version of C99's isnan, and only works for quiet nan's with IEEE-conforming compilers. More of a wish list really, I could really use all the C99 fp-testing macros (as C++ functions obviously).
Numeric constants: ~~~~~~~~~~~~~
Currently I've embedded the numeric constants used inside the implementation, a numeric constants library really would have simplified things... just a pity we could never agree on an interface last time this came up.
Whew, I think that's all, thanks for reading this far!
Regards,
John.
Thanks for the work! I'll look into details when they show up. Merci Hubert Holin
FWIW, asinh, acosh and atanh are already in Boost! If you wish, I can retire them (though I planed to add valarray support, I even printed a few days ago the discussion which had taken place on this two years (ahem) ago...).
I wouldn't dream of retiring them: however your implementations make assumptions that are only true for scalar types, and not complex numbers (like calling numeric_limits<T>::epsilon() and abs), my versions are strictly complex-number only. In other words they should complement each other. The other issue, is that calculating the real and imaginary parts separately is more accurate than treating the complex number "like a scalar" and doing it the obvious way.
In the mean time these algorithms can be downloaded as part of the TR1 library here: http://boost-sandbox.sourceforge.net/vault/index.php?&direction=0&order=&direc tory=tr1
(not showing up yet...)
Works for me... try going to http://boost-sandbox.sourceforge.net/vault/ and then navigating down into the tr1 directory.
This algorithm is part of C99, but by no means all compilers support it yet, I used a Taylor series expansion for small x: I'm aware that there are much more efficient methods, but optimizing compilers completely trash the logic of these (Intel C++ proved to be particularly bad).
I had planed to work on it (and a few related such as expm1...) something liketwo years ago and had planed to add it to the same special functions library. They really are important.
Agreed, it's just a pity that some compilers mangle the floating point logic so badly :-( Thanks for the interest, John.
Somewhere in the E.U., le 07/04/2005 Bonjour In article <00a701c539f1$91514090$55270d52@fuji>, "John Maddock" <john@johnmaddock.co.uk> wrote:
FWIW, asinh, acosh and atanh are already in Boost! If you wish, I can retire them (though I planed to add valarray support, I even printed a few days ago the discussion which had taken place on this two years (ahem) ago...).
I wouldn't dream of retiring them: however your implementations make assumptions that are only true for scalar types, and not complex numbers
Arrgh! I only coded the complex versions of sinc and sinhc!
(like calling numeric_limits<T>::epsilon() and abs), my versions are strictly complex-number only. In other words they should complement each other. The other issue, is that calculating the real and imaginary parts separately is more accurate than treating the complex number "like a scalar" and doing it the obvious way.
I will look at your implementation tomorrow, but I feel it would perhaps be simpler to replace my files with yours or merge the implementations (specialization for scalars and specialization for complexes), and update the library documentation. I'd also like to add valarray support, at least for the scalars (pending the language decisions concerning layout guaranties), at some point, though rereading the discussion of two years ago this may prove contentious. At any rate, I believe this is an enrichment of an existing, accepted, library, and as such does not warrant a separate review. Adding valarray support might warrant one, though, since there has been controversy (and thus the opportunity for improvement thru review).
In the mean time these algorithms can be downloaded as part of the TR1 library here: http://boost-sandbox.sourceforge.net/vault/index.php?&direction=0&order=&di rec tory=tr1
(not showing up yet...)
Works for me... try going to http://boost-sandbox.sourceforge.net/vault/ and then navigating down into the tr1 directory.
Got it! The problem I encountered (and which is still here) is that chose the "Sandbox CVS" link from the page at "http://www.boost.org", and then the "via the web" link, which lead me to a page at "http://www.boost.org/more/mailing_lists.htm#sandbox", which, going from there to "boost" and then "tr1" only locates "tuple". Choosing the "Sandbox Files" link from the page at "http://www.boost.org" works correctly however.
This algorithm is part of C99, but by no means all compilers support it yet, I used a Taylor series expansion for small x: I'm aware that there are much more efficient methods, but optimizing compilers completely trash the logic of these (Intel C++ proved to be particularly bad).
I had planed to work on it (and a few related such as expm1...) something liketwo years ago and had planed to add it to the same special functions library. They really are important.
Agreed, it's just a pity that some compilers mangle the floating point logic so badly :-(
Thanks for the interest,
John.
Hubert
Arrgh! I only coded the complex versions of sinc and sinhc!
:-) It happens.
I will look at your implementation tomorrow, but I feel it would perhaps be simpler to replace my files with yours or merge the implementations (specialization for scalars and specialization for complexes), and update the library documentation. I'd also like to add valarray support, at least for the scalars (pending the language decisions concerning layout guaranties), at some point, though rereading the discussion of two years ago this may prove contentious.
At any rate, I believe this is an enrichment of an existing, accepted, library, and as such does not warrant a separate review. Adding valarray support might warrant one, though, since there has been controversy (and thus the opportunity for improvement thru review).
OK, since my implementations are strictly complex number only, yours would still be needed for reals. BTW, I haven't implemented asinh/acosh as such: these are trivial variations of asin/acos (which are implemented). Just to be different, atan is a trivial variation of atanh: C99 expressed the relationship this way around, so I just followed.
I do not have institutional free access from my "beloved workplace" (grumble...).
No nor do I, but I've used this authors work before, and the current exchange rate makes ACM access not too painful...
Any others? You've got my curiosity roused...
In the stopgap series of implementations that my "special functions library" is, I would have simply used a succession of truncated series expansion, at least near zero (like I did for the various other functions), and then rely on whatever exponential function was coded in the user's platform's standard library.
Yep, that's easy enough: take a look at my log1p: the series evaluation is all boilerplate, so doing the same for expm1 is just a case of writing a short function object that calculates successive terms each time it's called. John.
log1p: ~~~~
This algorithm is part of C99, but by no means all compilers support it yet, I used a Taylor series expansion for small x: I'm aware that there are much more efficient methods, but optimizing compilers completely trash the logic of these (Intel C++ proved to be particularly bad).
I had planed to work on it (and a few related such as expm1...) something liketwo years ago and had planed to add it to the same special functions library. They really are important.
Agreed: how were you planning to implement expm1? The obvious series expansion: expm1(z) = SUM[k=1, INF] (z^k / k!) would converge in less than B steps for a B bit floating point type and |z| < 0.5 (in other words reasonable but not astounding performance), in fact it would be trivial to plug that into my existing series summation code. Alternatively there's a table based method (http://portal.acm.org/citation.cfm?doid=146847.146928) which claims to be particularly accurate. Any others? You've got my curiosity roused... John.
Somewhere in the E.U., le 07/04/2005 Bonjour In article <001101c53a90$8591e580$ec4c0e52@fuji>, "John Maddock" <john@johnmaddock.co.uk> wrote:
log1p: ~~~~
This algorithm is part of C99, but by no means all compilers support it yet, I used a Taylor series expansion for small x: I'm aware that there are much more efficient methods, but optimizing compilers completely trash the logic of these (Intel C++ proved to be particularly bad).
I had planed to work on it (and a few related such as expm1...) something liketwo years ago and had planed to add it to the same special functions library. They really are important.
Agreed: how were you planning to implement expm1? The obvious series expansion:
expm1(z) = SUM[k=1, INF] (z^k / k!)
would converge in less than B steps for a B bit floating point type and |z| < 0.5 (in other words reasonable but not astounding performance), in fact it would be trivial to plug that into my existing series summation code.
Alternatively there's a table based method (http://portal.acm.org/citation.cfm?doid=146847.146928) which claims to be particularly accurate.
I do not have institutional free access from my "beloved workplace" (grumble...).
Any others? You've got my curiosity roused...
John.
In the stopgap series of implementations that my "special functions library" is, I would have simply used a succession of truncated series expansion, at least near zero (like I did for the various other functions), and then rely on whatever exponential function was coded in the user's platform's standard library. Of course, one can only dream of writing a C++ implementation of NIST's DLMF (http://dlmf.nist.gov/), ready any decade now... (OK, officially this year...). On a related point, I recently discovered that, at least for the Sun platform where tests are carried out for Boost, "long double" has numeric limits set correctly (and with values different than those of "double"), but some standard library functions (such as the exponential, apparently) have no "long double" overload, and seemingly truncate to "double". I have to adapt my test to declare the library to pass (rather than fail, such as is now the case), with warnings of QOI. This problem might bite your implementation as well, and IMHO does militate for the above dream to come true! Hubert
participants (2)
-
Hubert Holin -
John Maddock