From: Paul A. Bristow (boost_at_[hidden])
Date: 2001-09-13 12:51:10
I have recently coded W Cody's algorithm, originally in FORTRAN,
ACM TOMS 715
using rational approximation, using C++,
as a sample test for submission for math functions.
Despite Eric Ford enthusiam for making things more complex,
I have for this kept it simple, so far, and only defined
double psi(double x);
(Cody noted that system libraries normally only provide one library
using double precision for IEEE conforming hardware, and I have
followed this. Any saving in space and time that a float version
would provide is small and in my view not worth the complexity.
At the other end of the scale, long double precision, the algorithms
are not available. Most users will find double ample precision.
All number theoretic workers are likely to use arbitrary precision,
and quite different algorithms.)
Several questions emerge from implementing this algorithm.
The function has some exciting behaviour around zero and negative integers,
when it zips to infinity and back.
For these integral values, it really should return infinity.
Should I return std::numeric_limits<double>::infinity()
(How many platforms do NOT support infinity?)
Nearby, it should return std::numeric_limits<double>::max()
Or should I just return std::numeric_limits<double>::max()
instead of infinity? (This is what Cody does).
Or how should I let the user decide?
(Eric may have some ideas here?)
I note that there are no restrictions on arguments
(unless you want to throw an exceptions if infinity would result).
PS Does any one have worries about a boost licence for
re-implementations of code which is ACM copyright?
Dr Paul A Bristow, hetp Chromatography
LA8 8AB UK
+44 1539 561830
Mobile +44 7714 33 02 04
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk