Math Functions

From: boost_at_[hidden]
Date: 2001-09-04 05:45:37

I am working on some math functions like erf, beta, gamma,
and have some questions on preferred practice  (hoping to
avoid lots of criticisms at a late stage)!

1  My instinct is to only produce a double version(64 bit reals),
rather than template versions.
This is because the polynomial coefficients and algorithms
are only readily available for double.
Most people will use these.
Requirement for float can be met by downsizing the double result.
(I realise that embedded solutions might want to use float,
but the size of a templated float version is likely to drive them
to brew their own  special short-cut versions anyway).

Anyone who wants long double is likely to use an even high precision
arithmetic package. (Perhaps NIST will do something on these lines
in their forthcoming project?)

2  Arguments

Many need checks on validity of arguments.
For example, probability p must be 0 <= p <= 1.

What should I do after the check fails?
throw standard exception out_of_range?

Should I distinguish between invalid_argument for wrong sign,
and out_or_range for too big?

I propose to provide the value of the offending argument,
and perhaps an indication of which argument (many functions have 
several).

Does anyone have any specially good examples of doing this?

3  Handling math errors like overflow/underflow.

MSVC provides _matherr, _isfinite, _isnan to help detect when
things have gone pear-shaped, but are NOT portable.

Some Standard limits things work, but not all:

Checking for denormalisation works OK.

        product = std::numeric_limits<double>::min() /2.; // should be
denormalised!
        cout << "std::numeric_limits<double>::min() /2. " << product
         << ( (product > 0 && product < 
std::numeric_limits<double>::min()) ? " IS
denormalized!" : " is NOT denormalised") << nl;

And so does checking for infinity:

        product = std::numeric_limits<double>::infinity();
        if (product == std::numeric_limits<double>::infinity() )
        { // This does work OK.
                throw My_overflow_error("infinity");
        }

But can't use this for NaNs???

        if (product == std::numeric_limits<double>::quiet_NaN() )
        {	// If compare with standard limits quiet_NaN, then 
result is ALWAYS a
NaN!!!
                throw My_underflow_error("NaN");
        }

Does any one know of PORTABLE ways of doing these things?  Examples 
please.

5  When bad things are detected, should I throw exceptions derived 
from the
standard exceptions like

class My_overflow_error : public std::overflow_error ...
class My_range_error : public std::range_error

5  When checking for an argument value like zero or unity,
should I also allow within an espilon of it?

if (p <= 0)  ...

or

if (p <= 0 + standard::numeric_limits<double>::epsilon() ) ...

(The logic of this is that computation error in previous stages might 
result
in
this instead of a exact zero).

6  If the result is infinity (for limit values like 0 and 1 for 
example)
should I return standard::numeric_limits<double>::infinity()

or standard::numeric_limits<double>::max()

Views and help please.

Paul