summaryrefslogtreecommitdiff
path: root/libm/float/dawsnf.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/float/dawsnf.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/float/dawsnf.c')
-rw-r--r--libm/float/dawsnf.c168
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 );
+}