/*************************************************** * Program to demonstrate ASYMERF * * ------------------------------------------------ * * Reference: BASIC Scientific Subroutines, Vol. II * * By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].* * * * C++ Version by J.-P. Moreau, Paris. * * (www.jpmoreau.fr) * * ------------------------------------------------ * * SAMPLE RUN: * * * * Find the value of ERF(X)=2*Exp(-X*X)/SQRT(PI) * * * * Input X ? 3 * * ERF(X)= .9999779 with error estimate= -.00000000 * * Number of terms evaluated was 10 * * * * Input X ? 4 * * ERF(X)= 1.0000000 with error estimate= 0.0000000 * * Number of terms evaluated was 17 * * * ***************************************************/ #include #include #define PI 4*atan(1) double e, x, y; int n; /********************************************************** * Asymptotic series expansion of the integral of * * 2 exp(-X*X)/(X*sqrt(PI)), the normalized error function * * (ASYMERF). This program determines the values of the * * above integrand using an asymptotic series which is * * evaluated to the level of maximum accuracy. * * The integral is from 0 to X. The input parameter, X * * must be > 0. The results are returned in Y and Y1, * * with the error measure in E. The number of terms used * * is returned in N. The error is roughly equal to first * * term neglected in the series summation. * * ------------------------------------------------------- * * Reference: A short table of integrals by B.O. Peirce, * * Ginn and Company, 1957. * **********************************************************/ void Asymerf() { // Labels: e100, e200 double c1, c2, y1; n = 1; y = 1; c2 = 1 / (2 * x * x); e100: y = y - c2; n = n + 2; c1 = c2; c2 = -c1 * n / (2 * x * x); // Test for divergence - The break point is roughly N=X*X if ( fabs(c2) > fabs(c1)) goto e200; // Continue summation goto e100; e200: n = (n + 1) / 2; e = exp(-x * x) / (x * sqrt(PI)); y1 = y * e; y = 1.0 - y1; e = e * c2; } void main() { printf("\n Input X: "); scanf("%lf",&x); Asymerf(); printf("\n ERF(%5.2f)= %9.7f with error estimate= %9.7f\n\n",x,y,e); printf(" Number of terms evaluated was %d.\n\n\n\n", n); } // End of file asymerf.cpp