/***************************************************
*          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 <stdio.h>
#include <math.h>

#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