Digital Mars - c++ - erf, erfc, and math.h

↑ ↓ ← → fprimex <fprimex_member pathlink.com> writes:

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

Jan 27 2004

↑ ↓ ← → "Walter" <walter digitalmars.com> writes:

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


 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


 one that has the error and complimentary error functions defined in


 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.


 these functions need to be written I would be happy to contribute.

 Brent W.

Jan 27 2004

↑ ↓ ← → "Walter" <walter digitalmars.com> writes:

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

Jan 28 2004

↑ ↓ ← → "Steve Strand" <snstrand comcast.net> writes:

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';
}

Jan 29 2004

↑ ↓ ← → "Steve Strand" <snstrand comcast.net> writes:

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';
}

Jan 29 2004

↑ ↓ ← → "Steve Strand" <snstrand comcast.net> writes:

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

Feb 01 2004

↑ ↓ ← → pbaraduc ion.ucl.ac.uk writes:

In article <bvbnom$2pc6$1 digitaldaemon.com>, Steve Strand says...
After reviewing my erf() code I see that I forgot to initialize q2 in erfc().
Find the corrected code below (one line is
different).


Actually there's far bigger bug in erf: the absolute value (fabs) is not applied
to the right term so erf gives wrong results for values between -2.2 and 0...

Here's the corrected code:

#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);   // CORRECTED LINE
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;
}



--
Pierre

Apr 26 2004