RE: [boost] Re: C++ Standard Library proposal -MathfunctionsforStatistics

From: Paul A. Bristow (boost_at_[hidden])
Date: 2003-04-24 13:43:45

| -----Original Message-----
| From: boost-bounces_at_[hidden]
| [mailto:boost-bounces_at_[hidden]]On Behalf Of Jason House
| Sent: Thursday, April 24, 2003 3:57 PM
| To: boost_at_[hidden]
| Subject: [boost] Re: C++ Standard Library proposal
| -MathfunctionsforStatistics
|
| Well, take erf using a single parameter for precision
|
template <precision = 0.00001>
double erf(double x)
{
        // ...
  // return erfx;
};

For many of the proposed functions, like this one erf, IMHO this is making an
unnecessary complexity.

The implementation is probably an evaluation of a polynomial and there is little
point in not returning the full precision for the floating point type. Possible
saving from not evaluating all coefficients will be small compared to the
overheads.

It probably uses some (purely imaginary) coefficients like:

static double P[] = {
        2.46196981473530512524E-10,
        5.64189564831068821977E-1,
        7.46321056442269912687E0,
        4.86371970985681366614E1,
        1.96520832956077098242E2,
};

static float P[] = { // note fewer terms needed for 32-bit precison?
        2.461969814E-10,
        5.6418956E-1,
        7.46321056E0,
        4.8637197E1,
};

static long double P[] = { // More coefficients for 80 or 128 bit precision?
        2.4619698147353051252446196981473530512524E-10,
        5.6418956483106882197746196981473530512524E-1,
        7.4632105644226991268746196981473530512524E0,
        4.8637197098568136461969814735305125246614E1,
        1.9652083295646196981473530512524077098242E2,
        5.2644519499547735846196981473530512524631E2,
        9.3446196981473530512524528527171957607540E2,
        1.0275518864619698147353051252489515710272E3,
        5.5753533461969814735305125245369399327526E2
};

Of course, there may be different sets of coefficients for float, double and
long double.
So three templated versions are all I believe are needed for this function.

FPType erf<FPType>(FPType)
{
// ... evaluate polynomial using the set of coefficient for FPType.
        return erfFPType;
}

float erff = erf<float>(1.234F);
double erfd = erf<double>(1.234); // erf(1.234) will chose this of course.
long double erfl = erf<long double>(1.234L);

One may want to calculate with a higher precision like double and simply return
a float( erf<double>(1.234)) which is probably accurate to the full 32-bit.

C versions: float erff(float x){ return erf<float>(x); }  ...

Any rounding (for output) should take place as late as possible.

Of course, most of the inverse functions will be implemented using some
iterative Newton-style root finding, and here there IS a compromise between
accuracy and speed.  Even here I am unconvinced that we need the complexity of
precision control (though Eric Ford provided a different view with a complex
system).

For most purposes, one wants the natural precision the floating point type
provides (though sadly for some functions this is several decimal places short
of the full precision for that type - but much better than the A&S tables).

Paul

Paul A Bristow, Prizet Farmhouse, Kendal, Cumbria, LA8 8AB  UK
[phone]   Mobile [phone]
Mobile mailto:pabristow_at_[hidden]
mailto:pbristow_at_[hidden]