Return to BSD News archive
Path: sserve!manuel.anu.edu.au!munnari.oz.au!uunet!ukma!cs.widener.edu!dsinc!ub!acsu.buffalo.edu!usenet
From: jones@acsu.buffalo.edu (Terry A. Jones)
Newsgroups: comp.unix.bsd
Subject: Re: Where are the functions erf(), and erfc()?
Message-ID: <Bx9vtF.2Bo@acsu.buffalo.edu>
Date: 6 Nov 92 02:20:02 GMT
References: <Bwz4qG.9p6@acsu.buffalo.edu>
Sender: nntp@acsu.buffalo.edu
Organization: State University of New York at Buffalo/Comp Sci
Lines: 117
Nntp-Posting-Host: hyades.cs.buffalo.edu
Well I did not get any feedback on this one so I implemented
the functions myself. I couldn't wait to get Spice compiled. Cut into
erf.c, compile and add to your libm.a if you like.
Terry Jones jones@cs.buffalo.edu
---------------------------cut here-----------------------------------------
#include "mathimpl.h"
/*
* Predefined terms in expansion series for the error function.
*/
#define T2 0.6666667
#define T3 0.2666667
#define T4 0.07619048
#define T5 0.01693122
#define T6 3.078403e-3
#define T7 4.736005e-4
#define T8 6.314673e-5
#define T9 7.429027e-6
#define T10 7.820028e-7
#define T11 7.447646e-8
#define T12 6.476214e-9
/*
* Internal version of the error function.
*/
static double errorf( x )
double x;
{
double x2;
x2 = x * x;
return( ( 2.0 * exp( -x2 ) / sqrt( M_PI ) ) * ( x * ( 1 + x2 * ( T2 + x2 *
( T3 + x2 *( T4 + x2 *( T5 + x2 *( T6 + x2 *( T7 + x2 *( T8 + x2 *
( T9 + x2 *( T10 + x2 *( T11 + x2 * T12)))))))))))));
}
/*
* Compute the complement to the error function.
*/
static double errorf_complement( x )
double x;
{
double x2;
double v;
x2 = x * x;
v = 1.0 / ( 2.0 * x2 );
return( 1.0 / ( exp( x2 ) * x * sqrt( M_PI ) * ( 1 + v / ( 1 + 2 * v /
( 1 + 3 * v / ( 1 + 4 * v / ( 1 + 5 * v / ( 1 + 6 * v /
( 1 + 7 * v / ( 1 + 8 * v / ( 1 + 9 * v / ( 1 + 10 * v /
( 1 + 11 * v / ( 1 + 12 * v / ( 1 + 13 * v / ( 1 + 14 * v /
( 1 + 15 * v / ( 1 + 16 * v / ( 1 + 17 * v )))))))))))))))))));
}
/*
* Global entry point for the error function.
*/
double xerf( x )
double x;
{
double sign;
if ( 0.0 == x )
return( 0.0 );
/*
* Handle case of negative argument.
*/
if ( x < 0.0 )
{
sign = -1.0;
x = -x;
}
else
sign = 1.0;
if ( x < 1.5 )
return( errorf( x ) * sign );
else
return( ( 1.0 - errorf_complement( x )) * sign );
}
/*
* Global entry point for the error function complement.
*/
double xerfc( x )
double x;
{
if ( 0.0 == x )
return( 1.0 );
if ( x < 0.0 )
{
x = -x;
if ( x < 1.5 )
return( 1.0 + errorf( x ) );
else
return( 2.0 - errorf_complement( x ) );
}
else
{
if ( x < 1.5 )
return( 1.0 - errorf( x ) );
else
return( errorf_complement( x ) );
}
}