summaryrefslogtreecommitdiff
path: root/libm/float/dawsnf.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/float/dawsnf.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/float/dawsnf.c')
-rw-r--r--libm/float/dawsnf.c168
1 files changed, 0 insertions, 168 deletions
diff --git a/libm/float/dawsnf.c b/libm/float/dawsnf.c
deleted file mode 100644
index d00607719..000000000
--- a/libm/float/dawsnf.c
+++ /dev/null
@@ -1,168 +0,0 @@
-/* dawsnf.c
- *
- * Dawson's Integral
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, dawsnf();
- *
- * y = dawsnf( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Approximates the integral
- *
- * x
- * -
- * 2 | | 2
- * dawsn(x) = exp( -x ) | exp( t ) dt
- * | |
- * -
- * 0
- *
- * Three different rational approximations are employed, for
- * the intervals 0 to 3.25; 3.25 to 6.25; and 6.25 up.
- *
- *
- * ACCURACY:
- *
- * Relative error:
- * arithmetic domain # trials peak rms
- * IEEE 0,10 50000 4.4e-7 6.3e-8
- *
- *
- */
-
-/* dawsn.c */
-
-
-/*
-Cephes Math Library Release 2.1: January, 1989
-Copyright 1984, 1987, 1989 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-*/
-
-#include <math.h>
-/* Dawson's integral, interval 0 to 3.25 */
-static float AN[10] = {
- 1.13681498971755972054E-11,
- 8.49262267667473811108E-10,
- 1.94434204175553054283E-8,
- 9.53151741254484363489E-7,
- 3.07828309874913200438E-6,
- 3.52513368520288738649E-4,
--8.50149846724410912031E-4,
- 4.22618223005546594270E-2,
--9.17480371773452345351E-2,
- 9.99999999999999994612E-1,
-};
-static float AD[11] = {
- 2.40372073066762605484E-11,
- 1.48864681368493396752E-9,
- 5.21265281010541664570E-8,
- 1.27258478273186970203E-6,
- 2.32490249820789513991E-5,
- 3.25524741826057911661E-4,
- 3.48805814657162590916E-3,
- 2.79448531198828973716E-2,
- 1.58874241960120565368E-1,
- 5.74918629489320327824E-1,
- 1.00000000000000000539E0,
-};
-
-/* interval 3.25 to 6.25 */
-static float BN[11] = {
- 5.08955156417900903354E-1,
--2.44754418142697847934E-1,
- 9.41512335303534411857E-2,
--2.18711255142039025206E-2,
- 3.66207612329569181322E-3,
--4.23209114460388756528E-4,
- 3.59641304793896631888E-5,
--2.14640351719968974225E-6,
- 9.10010780076391431042E-8,
--2.40274520828250956942E-9,
- 3.59233385440928410398E-11,
-};
-static float BD[10] = {
-/* 1.00000000000000000000E0,*/
--6.31839869873368190192E-1,
- 2.36706788228248691528E-1,
--5.31806367003223277662E-2,
- 8.48041718586295374409E-3,
--9.47996768486665330168E-4,
- 7.81025592944552338085E-5,
--4.55875153252442634831E-6,
- 1.89100358111421846170E-7,
--4.91324691331920606875E-9,
- 7.18466403235734541950E-11,
-};
-
-/* 6.25 to infinity */
-static float CN[5] = {
--5.90592860534773254987E-1,
- 6.29235242724368800674E-1,
--1.72858975380388136411E-1,
- 1.64837047825189632310E-2,
--4.86827613020462700845E-4,
-};
-static float CD[5] = {
-/* 1.00000000000000000000E0,*/
--2.69820057197544900361E0,
- 1.73270799045947845857E0,
--3.93708582281939493482E-1,
- 3.44278924041233391079E-2,
--9.73655226040941223894E-4,
-};
-
-
-extern float PIF, MACHEPF;
-#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
-#ifdef ANSIC
-float polevlf(float, float *, int);
-float p1evlf(float, float *, int);
-#else
-float polevlf(), p1evlf();
-#endif
-
-float dawsnf( float xxx )
-{
-float xx, x, y;
-int sign;
-
-xx = xxx;
-sign = 1;
-if( xx < 0.0 )
- {
- sign = -1;
- xx = -xx;
- }
-
-if( xx < 3.25 )
- {
- x = xx*xx;
- y = xx * polevlf( x, AN, 9 )/polevlf( x, AD, 10 );
- return( sign * y );
- }
-
-
-x = 1.0/(xx*xx);
-
-if( xx < 6.25 )
- {
- y = 1.0/xx + x * polevlf( x, BN, 10) / (p1evlf( x, BD, 10) * xx);
- return( sign * 0.5 * y );
- }
-
-
-if( xx > 1.0e9 )
- return( (sign * 0.5)/xx );
-
-/* 6.25 to infinity */
-y = 1.0/xx + x * polevlf( x, CN, 4) / (p1evlf( x, CD, 5) * xx);
-return( sign * 0.5 * y );
-}