#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

double Erf( double x )
{
	if( x == 0.0 )
		return 0.0;

	double		f, fZm, fN2, fn1, fOld, fT;
	bool	bAdd = true;
	long	nMaxIter = 1000;

	fZm = x;
	fZm *= x;	// x^2

	fn1 = 2.0;
	fN2 = 3.0;

	fT = x * fZm;

	f = x - fT / fN2;

	fOld = f * 0.9;

	while( f != fOld && nMaxIter )
	{
		fOld = f;

		fN2 += 2.0;

		fT /= fn1;
		fT *= fZm;

		fn1++;

		if( bAdd )
			f += fT / fN2;
		else
			f -= fT / fN2;

		bAdd = !bAdd;

		nMaxIter--;
	}

	f *= 1.128379167095512573896;
	return f;
}

//# define M_2_SQRTPIl	1.1283791670955125738961589031215452L  /* 2/sqrt(pi) */

void ErfAsymptotic( double x, double& fErf, double& fErfc )
{
	if( x == 0.0 )
	{
		fErf  = 0.0;
		fErfc = 1.0;
		return;
	}

	bool bNegative = false;
	if ( x < 0.0 )
	{
		x = fabs( x );
		bNegative = true;
	}

	long nMaxIter = 1000;

	double fFct2 = 1.0, fDenom = 1.0, fX = x, fn = 0.0;
	double fXCnt = 1.0;
	double f = 1.0 / fX;
	double fOld = f * 0.9;

	bool bAdd = false;
	double fDiff = 0.0, fDiffOld = 1000.0;
	while( f != fOld && nMaxIter-- )
	{
		fOld = f;

		fFct2 *= 2.0*++fn - 1.0;
		fDenom *= 2.0;
		fX *= x*x;
		fXCnt += 2.0;

		if ( bAdd )
			f += fFct2 / ( fDenom * fX );
		else
			f -= fFct2 / ( fDenom * fX );

		bAdd = !bAdd;
		
		fDiff = fabs( f - fOld );
		if ( fDiffOld < fDiff )
		{
			f = fOld;
			break;
		}
		fDiffOld = fDiff;
	}

	fErf = 1.0 - f*exp( -1.0*x*x )*0.5641895835477562869480794515607726;
	if ( bNegative )
		fErf *= -1.0;
	fErfc = 1.0 - fErf;
}

int main( int argc, char** argv )
{
	
	if ( argc < 2 )
	{
		cout << "x\terf(x) <maclaurin>\terfc(x) <maclaurin>\terf(x) <asymptotic>\terfc(x) <asymptotic>" << endl;

		for ( double i = -10.0; i <= 10.0; i += 0.1 )
		{
			cout << setprecision( 2 ) << i;
			double fErf = Erf( i );
			double fErfc = 1.0 - fErf;
			
			cout << setprecision( 16 ) << "\t" << fErf << "\t" << fErfc << "\t";
			
			ErfAsymptotic( i, fErf, fErfc );
			cout << fErf << "\t" << fErfc << endl;
		}	
	}
	else
	{
		cout << setprecision( 16 );
		double fErf, fErfc;
		double x = atof( argv[1] );
		fErf = Erf( x );
		fErfc = 1.0 - fErf;
		cout << "<maclaurin>" << endl;
		cout << "erf:  " << fErf << endl;
		cout << "erfc: " << fErfc << endl;

		ErfAsymptotic( x, fErf, fErfc );
		cout << "<asymptotic>" << endl;
		cout << "erf:  " << fErf << endl;
		cout << "erfc: " << fErfc << endl;
	}

	return EXIT_SUCCESS;
}