Re: Math Functions: digamm functions, isnan and infinity

From: Eric Ford (eford_at_[hidden])
Date: 2001-09-13 14:00:08

> 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.

Rational function approximations will likely be common.  I hope this
means you implemented a rational function function/class that can be
easily reused for approximating other functions (e.g. they're used by
the error function routine you posted earlier).  I would think a
rational function function/class based on the boost::array could be
efficient 

> At the other end of the scale, long double precision, the algorithms
> (coefficients)
> are not available.  Most users will find double ample precision.
> All number theoretic workers are likely to use arbitrary precision,
> and quite different algorithms.)

Coefficients aren't magical.  They come from somewhere.  Does your
source explain how they were calculated or a reference to another
source that does?

> For these integral values, it really should return infinity.
> 
> Should I return std::numeric_limits<double>::infinity()
> Or should I just return std::numeric_limits<double>::max()
> instead of infinity?  (This is what Cody does).

I think it should return inf, iff the correct mathematical answer
really is inf.

> Nearby, it should return std::numeric_limits<double>::max()

The important thing is that let the user know it isn't giving the
right answer, or at least give the user the option of testing this. 
Whether that's done through a designated error value, a global
variable, throwing an exception, etc..  As you know, I think it would
be nice if the user could choose which of these facilities is used.
But that's not as important as there being some way for the user to
find out.  Otherwise there's a risk they'll get a wrong answer and not
know it.

> Or how should I let the user decide?

I'd recommend a template argument.  If it's a yes no choice, then it
can be a bool.  If you want to provide multiple choices, you can use
an int with some constants/enums to help make the code legible.  If
you want to make it fully user customizable, then I think a class
providing a couple of functions is the way to go.

> I note that there are no restrictions on arguments
> (unless you want to throw an exceptions if infinity would result).

Well, it definitely doesn't make sense for a nan.  If I remember
correctly, \lim_{x \to -\infnty} \psi(x) isn't defined.  And I'm
guessing your algorithm is only for reals.  As long as your function
is strictly double digamma(double), this isn't an issue.  However, if
I were extending it to take and return a type specified by a template
argument, I'd need to make sure either that the function could deal
with complex arguments or that it would generate a compile-time
error.  

> PS Does any one have worries about a boost licence for
> re-implementations of code which is ACM copyright?

What restrictions do they place on their code?  

Thanks,
E