erf(x) and math.h

Question

According to this site the error function erf(x) comes from math.h. But actually looking in math.h, it isn't there, and gcc cannot compile the following test program while g++ can:

#include <math.h>
#include <stdio.h>

int main(int argc, char* argv[]) {
  double x;
  double erfX;
  x = 1.0;
  erfX = erf(x);

  printf("erf(%f) = %f", x, erfX);
}

$ gcc mathHTest.c
/tmp/ccWfNox5.o: In function `main':
mathHTest.c:(.text+0x28): undefined reference to `erf'
collect2: ld returned 1 exit status
$ g++ mathHTest.c

What does g++ pull in that gcc doesn't? Looking in /usr/include, the only place I could find erf(x) was in tgmath.h, which I don't include. So g++ must be grabbing different headers than gcc, but which ones?

EDIT: I wasn't linking in libm with gcc, hence the link error. However, I still don't understand why erf() is not in math.h. Where is it coming from?

Answer 1

I had a similar problem and needed to find the exact definition of erf so let me expand on this. As said by Chris Dodd the function is declared in bits/mathcalls.h which is included by maths.h.

bits/mathcalls.h:

...
#if defined __USE_MISC || defined __USE_XOPEN || defined __USE_ISOC99
__BEGIN_NAMESPACE_C99
/* Error and gamma functions.  */
__MATHCALL (erf,, (_Mdouble_));
__MATHCALL (erfc,, (_Mdouble_));
__MATHCALL (lgamma,, (_Mdouble_));
__END_NAMESPACE_C99
#endif
...

Macro magic expands __MATHCALL (erf,, (_Mdouble_)); to

extern double erf (double) throw (); extern double __erf (double) throw ();

The actual code is in libm.a or libm.so (gcc -lm):

$ nm /usr/lib/libm.a
...
s_erf.o:
00000400 T __erf
00000000 T __erfc
         U __ieee754_exp
00000400 W erf
00000000 W erfc
...

The source can be obtained from the gnu libc webpage. For a rough idea on the actual implementation here a few lines of the source:

sysdeps/ieee754/dbl-64/s_erf.c:

/* double erf(double x)
 * double erfc(double x)
 *                           x
 *                    2      |\
 *     erf(x)  =  ---------  | exp(-t*t)dt
 *                 sqrt(pi) \|
 *                           0
 *
 *     erfc(x) =  1-erf(x)
 *  Note that
 *              erf(-x) = -erf(x)
 *              erfc(-x) = 2 - erfc(x)
 *
 * Method:
 *      1. For |x| in [0, 0.84375]
 *          erf(x)  = x + x*R(x^2)
 *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
 *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
 *         where R = P/Q where P is an odd poly of degree 8 and
 *         Q is an odd poly of degree 10.
 *                                               -57.90
 *                      | R - (erf(x)-x)/x | <= 2
 *
 *
 *         Remark. The formula is derived by noting
 *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
 *         and that
 *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
 *         is close to one. The interval is chosen because the fix
 *         point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
 *         near 0.6174), and by some experiment, 0.84375 is chosen to
 *         guarantee the error is less than one ulp for erf.
 *
 *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and     
 ...

Answer 2

'erf' is actually declared in bits/mathcalls.h, which is #included by math.h. The actual declaration is heavily obscured by macro magic to make it do the right thing for both C and C++