1 files changed, 168 insertions, 0 deletions
diff --git a/libm/float/dawsnf.c b/libm/float/dawsnf.c
new file mode 100644
index 000000000..d00607719
--- /dev/null
+++ b/libm/float/dawsnf.c
@@ -0,0 +1,168 @@
+/*							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 );
+}