(posted on c++.stl - then realized no one has posted there in a long time)

Hi folks,

I'm having a little trouble here - hope this is the right place to post. I am used to using gcc/g++ under Linux and Mac OS X. This code compiles no problem under gcc:

/* --- Begin code snippet --- */

/* --- erftest.c --- */

#include <stdlib.h>
#include <math.h>

int main()
{
printf("%g\n", erf(1));
return 0;
}

/* --- end --- */

I have checked out a number of compilers for win32 and DM looks like the only one that has the error and complimentary error functions defined in math.h. I need these to complete my win32 port. Using digital mars to compile I get:

erftest.obj(erftest)
Error 42: Symbol Undefined _erf

--- errorlevel 1


Am I doing something wrong? Any help would be greatly appreciated. I'm not using STLport as it does not have erf or any other C99 extented functions. If these functions need to be written I would be happy to contribute.

Brent W.

The erf() function is not implemented yet. Sorry (you're the first to ask
for it!).

"fprimex" <fprimex_member@pathlink.com> wrote in message news:bv7cmq$1ie9$1@digitaldaemon.com...
> (posted on c++.stl - then realized no one has posted there in a long time)
>
> Hi folks,
>
> I'm having a little trouble here - hope this is the right place to post. I
am
> used to using gcc/g++ under Linux and Mac OS X. This code compiles no problem under gcc:
>
> /* --- Begin code snippet --- */
>
> /* --- erftest.c --- */
>
> #include <stdlib.h>
> #include <math.h>
>
> int main()
> {
> printf("%g\n", erf(1));
> return 0;
> }
>
> /* --- end --- */
>
> I have checked out a number of compilers for win32 and DM looks like the
only
> one that has the error and complimentary error functions defined in
math.h. I
> need these to complete my win32 port. Using digital mars to compile I get:
>
> erftest.obj(erftest)
> Error 42: Symbol Undefined _erf
>
> --- errorlevel 1
>
>
> Am I doing something wrong? Any help would be greatly appreciated. I'm not using STLport as it does not have erf or any other C99 extented functions.
If
> these functions need to be written I would be happy to contribute.
>
> Brent W.
>
>

Any chance it's coming soon? As I said, I'd be happy to contribute.

Brent W.

In article <bv7iev$1se8$3@digitaldaemon.com>, Walter says...
>
>The erf() function is not implemented yet. Sorry (you're the first to ask
>for it!).
>

I don't know. There is a lot to do!

"fpirmex" <fpirmex_member@pathlink.com> wrote in message news:bv8t25$11v0$1@digitaldaemon.com...
> Any chance it's coming soon? As I said, I'd be happy to contribute.
>
> Brent W.
>
> In article <bv7iev$1se8$3@digitaldaemon.com>, Walter says...
> >
> >The erf() function is not implemented yet. Sorry (you're the first to ask
> >for it!).
> >
>
>

I have had to implement erf() and erfc() before so I remember exactly how to do it. Here is the code.

P.S. to Walter: feel free to add this code to the DM library

/**************************
*   erf.cpp
*   author:  Steve Strand
*   written: 29-Jan-04
***************************/

#include <iostream.h>
#include <iomanip.h>
#include <strstream.h>
#include <math.h>


static const double rel_error= 1E-12;        //calculate 12 significant figures
//you can adjust rel_error to trade off between accuracy and speed
//but don't ask for > 15 figures (assuming usual 52 bit mantissa in a double)


double erf(double x)
//erf(x) = 2/sqrt(pi)*integral(exp(-t^2),t,0,x)
//         = 2/sqrt(pi)*[x - x^3/3 + x^5/5*2! - x^7/7*3! + ...]
//         = 1-erfc(x)
{
    static const double two_sqrtpi=  1.128379167095512574;        // 2/sqrt(pi)
    if (fabs(x) > 2.2) {
        return 1.0 - erfc(x);        //use continued fraction when fabs(x) > 2.2
    }
    double sum= x, term= x, xsqr= x*x;
    int j= 1;
    do {
        term*= xsqr/j;
        sum-= term/(2*j+1);
        ++j;
        term*= xsqr/j;
        sum+= term/(2*j+1);
        ++j;
    } while (fabs(term)/sum > rel_error);
    return two_sqrtpi*sum;
}


double erfc(double x)
//erfc(x) = 2/sqrt(pi)*integral(exp(-t^2),t,x,inf)
//           = exp(-x^2)/sqrt(pi) * [1/x+ (1/2)/x+ (2/2)/x+ (3/2)/x+ (4/2)/x+ ...]
//           = 1-erf(x)
//expression inside [] is a continued fraction so '+' means add to denominator only
{
    static const double one_sqrtpi=  0.564189583547756287;        // 1/sqrt(pi)
    if (fabs(x) < 2.2) {
        return 1.0 - erf(x);        //use series when fabs(x) < 2.2
    }
    if (signbit(x)) {               //continued fraction only valid for x>0
        return 2.0 - erfc(-x);
    }
    double a=1, b=x;                //last two convergent numerators
    double c=x, d=x*x+0.5;          //last two convergent denominators
    double q1,q2;                   //last two convergents (a/c and b/d)
    double n= 1.0, t;
    do {
        t= a*n+b*x;
        a= b;
         b= t;
        t= c*n+d*x;
        c= d;
        d= t;
        n+= 0.5;
        q1= q2;
        q2= b/d;
      } while (fabs(q1-q2)/q2 > rel_error);
    return one_sqrtpi*exp(-x*x)*q2;
}


int main(int argc, char *argv[])
//print table of erf(x) and erfc(x)
{
    double x0= 0.0;                        // starting x
    if (argc>1)
        {istrstream(argv[1]) >> x0;}

    double x1= 1.0;                        // ending x
    if (argc>2)
        {istrstream(argv[2]) >> x1;}

    double dx= 0.1;                         // x increment
    if (argc>3)
        {istrstream(argv[3]) >> dx;}

    cout.precision(10);
    cout << "       x              erf(x)          erfc(x)\n";
    cout << "---------------  ---------------  ---------------\n";

    for (double x= x0; x<x1+dx/2; x+= dx)
        cout << setw(15) << x << setw(17) << erf(x) << setw(17) << erfc(x) << '\n';
}

Thanks!

After reviewing my erf() code I see that I forgot to initialize q2 in erfc(). Find the corrected code below (one line is
different).

The chance of the uninitialized 'q' matching the first computed 'q' to within 10^-12 and causing premature loop exit is about 10^-52 (since 12 exponent bits and 40 mantissa bits have to match) but of course if the code is released uncorrected some user will manage to do it instantly.

/***************************
*   erf.cpp
*   author:  Steve Strand
*   written: 29-Jan-04
***************************/

#include <iostream.h>
#include <iomanip.h>
#include <strstream.h>
#include <math.h>


static const double rel_error= 1E-12;        //calculate 12 significant figures
//you can adjust rel_error to trade off between accuracy and speed
//but don't ask for > 15 figures (assuming usual 52 bit mantissa in a double)


double erf(double x)
//erf(x) = 2/sqrt(pi)*integral(exp(-t^2),t,0,x)
//       = 2/sqrt(pi)*[x - x^3/3 + x^5/5*2! - x^7/7*3! + ...]
//       = 1-erfc(x)
{
    static const double two_sqrtpi=  1.128379167095512574;        // 2/sqrt(pi)
    if (fabs(x) > 2.2) {
        return 1.0 - erfc(x);        //use continued fraction when fabs(x) > 2.2
    }
    double sum= x, term= x, xsqr= x*x;
    int j= 1;
    do {
        term*= xsqr/j;
        sum-= term/(2*j+1);
        ++j;
        term*= xsqr/j;
        sum+= term/(2*j+1);
        ++j;
    } while (fabs(term)/sum > rel_error);
    return two_sqrtpi*sum;
}


double erfc(double x)
//erfc(x) = 2/sqrt(pi)*integral(exp(-t^2),t,x,inf)
//        = exp(-x^2)/sqrt(pi) * [1/x+ (1/2)/x+ (2/2)/x+ (3/2)/x+ (4/2)/x+ ...]
//        = 1-erf(x)
//expression inside [] is a continued fraction so '+' means add to denominator only
{
    static const double one_sqrtpi=  0.564189583547756287;        // 1/sqrt(pi)
    if (fabs(x) < 2.2) {
        return 1.0 - erf(x);        //use series when fabs(x) < 2.2
    }
    if (signbit(x)) {               //continued fraction only valid for x>0
        return 2.0 - erfc(-x);
    }
    double a=1, b=x;                //last two convergent numerators
    double c=x, d=x*x+0.5;          //last two convergent denominators
    double q1, q2= b/d;             //last two convergents (a/c and b/d)
    double n= 1.0, t;
    do {
        t= a*n+b*x;
        a= b;
        b= t;
        t= c*n+d*x;
        c= d;
        d= t;
        n+= 0.5;
        q1= q2;
        q2= b/d;
      } while (fabs(q1-q2)/q2 > rel_error);
    return one_sqrtpi*exp(-x*x)*q2;
}


int main(int argc, char *argv[])
//print table of erf(x) and erfc(x)
{
    double x0= 0.0;                        // starting x
    if (argc>1)
        {istrstream(argv[1]) >> x0;}

    double x1= 1.0;                        // ending x
    if (argc>2)
        {istrstream(argv[2]) >> x1;}

    double dx= 0.1;                         // x increment
    if (argc>3)
        {istrstream(argv[3]) >> dx;}

    cout.precision(10);
    cout << "       x              erf(x)          erfc(x)\n";
    cout << "---------------  ---------------  ---------------\n";

    for (double x= x0; x<x1+dx/2; x+= dx)
        cout << setw(15) << x << setw(17) << erf(x) << setw(17) << erfc(x) << '\n';
}

Great! Do you also have some test vectors for it?

Here is the data I used to test my code for erf() and erfc().
This test data was calculated by the symbolic algebra program
MuPAD with DIGITS set to 20. I compared against my code
run with rel_error set to 1E-12 and found no problems.

For small x test erf(x) which is near 0.  erfc(x)= 1-erf(x) is
near 1 and cannot carry as many significant figures

For x > 2 or so test erfc(x) which is near 0. erf(x)= 1-erfc(x)
is near 1 and cannot carry as many significant figures

erf(x) for small x
1.0e-1,  1.1246291601828489220e-1
1.0e-2,  1.1283415555849616916e-2
1.0e-3,  1.1283787909692363799e-3
1.0e-4,  1.1283791633342486949e-4
1.0e-5,  1.1283791670578999350e-5
1.0e-6,  1.1283791670951364475e-6
1.0e-7,  1.1283791670955088126e-7
1.0e-8,  1.1283791670955125363e-8
1.0e-9,  1.1283791670955125735e-9
1.0e-10, 1.1283791670955125739e-10

erf(x) for 0 to 2
0.0, 0.00000000000000000000
0.1, 0.11246291601828489220
0.2, 0.22270258921047845414
0.3, 0.32862675945912742764
0.4, 0.42839235504666845510
0.5, 0.52049987781304653768
0.6, 0.60385609084792592256
0.7, 0.67780119383741847298
0.8, 0.74210096470766048617
0.9, 0.79690821242283212852
1.0, 0.84270079294971486934
1.1, 0.88020506957408169977
1.2, 0.91031397822963538024
1.3, 0.93400794494065243660
1.4, 0.95228511976264881052
1.5, 0.96610514647531072707
1.6, 0.97634838334464400777
1.7, 0.98379045859077456363
1.8, 0.98909050163573071418
1.9, 0.99279042923525746995
2.0, 0.99532226501895273416

erfc(x) for 1 to 5
1.0, 1.5729920705028513066e-1
1.1, 1.1979493042591830023e-1
1.2, 8.9686021770364619762e-2
1.3, 6.5992055059347563396e-2
1.4, 4.7714880237351189484e-2
1.5, 3.3894853524689272933e-2
1.6, 2.3651616655355992226e-2
1.7, 1.6209541409225436374e-2
1.8, 1.0909498364269285816e-2
1.9, 7.2095707647425300516e-3
2.0, 4.6777349810472658379e-3
2.1, 2.9794666563329855039e-3
2.2, 1.8628462979818914435e-3
2.3, 1.1431765973566514653e-3
2.4, 6.8851389664507856974e-4
2.5, 4.0695201744495893956e-4
2.6, 2.3603441652934920399e-4
2.7, 1.3433273994052432914e-4
2.8, 7.5013194665459024223e-5
2.9, 4.1097878099458835684e-5
3.0, 2.2090496998585441373e-5
3.1, 1.1648657367199596034e-5
3.2, 6.0257611517620949717e-6
3.3, 3.0577097964381614618e-6
3.4, 1.5219933628622853618e-6
3.5, 7.4309837234141274552e-7
3.6, 3.5586299300768529882e-7
3.7, 1.6715105790914620237e-7
3.8, 7.7003927456964128699e-8
3.9, 3.4792248597231742279e-8
4.0, 1.5417257900280018853e-8
4.1, 6.7000276540848983736e-9
4.2, 2.8554941795921886166e-9
4.3, 1.1934717937220413048e-9
4.4, 4.8917102706058884270e-10
4.5, 1.9661604415428874870e-10
4.6, 7.7495995974418319916e-11
4.7, 2.9952597863796604121e-11
4.8, 1.1352143584921962156e-11
4.9, 4.2189365240057826046e-12
5.0, 1.5374597944280357689e-12

erfc(x) for 5 to 20
 5.0, 1.5374597944280356932e-12
 6.0, 2.1519736712499091684e-17
 7.0, 4.1838256077794143987e-23
 8.0, 1.1224297172982927080e-29
 9.0, 4.1370317465138102381e-37
10.0, 2.0884875837625447570e-45
11.0, 1.4408661379436946803e-54
12.0, 1.3562611692059042128e-64
13.0, 1.7395573154667245218e-75
14.0, 3.0372298477503116651e-87
15.0, 7.2129941724512066666e-100
16.0, 2.3284857515715306934e-113
17.0, 1.0212280150942608811e-127
18.0, 6.0823692318163993077e-143
19.0, 4.9177228392564754464e-159
20.0, 5.3958656116079009289e-176

That'll do. Thanks!