summaryrefslogtreecommitdiff
path: root/libm/float/struvef.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/struvef.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/float/struvef.c')
-rw-r--r--libm/float/struvef.c315
1 files changed, 315 insertions, 0 deletions
diff --git a/libm/float/struvef.c b/libm/float/struvef.c
new file mode 100644
index 000000000..4cf8854ed
--- /dev/null
+++ b/libm/float/struvef.c
@@ -0,0 +1,315 @@
+/* struvef.c
+ *
+ * Struve function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float v, x, y, struvef();
+ *
+ * y = struvef( v, x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Computes the Struve function Hv(x) of order v, argument x.
+ * Negative x is rejected unless v is an integer.
+ *
+ * This module also contains the hypergeometric functions 1F2
+ * and 3F0 and a routine for the Bessel function Yv(x) with
+ * noninteger v.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * v varies from 0 to 10.
+ * Absolute error (relative error when |Hv(x)| > 1):
+ * arithmetic domain # trials peak rms
+ * IEEE -10,10 100000 9.0e-5 4.0e-6
+ *
+ */
+
+
+/*
+Cephes Math Library Release 2.2: July, 1992
+Copyright 1984, 1987, 1989 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+#include <math.h>
+#define DEBUG 0
+
+extern float MACHEPF, MAXNUMF, PIF;
+
+#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
+
+#ifdef ANSIC
+float gammaf(float), powf(float, float), sqrtf(float);
+float yvf(float, float);
+float floorf(float), ynf(int, float);
+float jvf(float, float);
+float sinf(float), cosf(float);
+#else
+float gammaf(), powf(), sqrtf(), yvf();
+float floorf(), ynf(), jvf(), sinf(), cosf();
+#endif
+
+float onef2f( float aa, float bb, float cc, float xx, float *err )
+{
+float a, b, c, x, n, a0, sum, t;
+float an, bn, cn, max, z;
+
+a = aa;
+b = bb;
+c = cc;
+x = xx;
+an = a;
+bn = b;
+cn = c;
+a0 = 1.0;
+sum = 1.0;
+n = 1.0;
+t = 1.0;
+max = 0.0;
+
+do
+ {
+ if( an == 0 )
+ goto done;
+ if( bn == 0 )
+ goto error;
+ if( cn == 0 )
+ goto error;
+ if( (a0 > 1.0e34) || (n > 200) )
+ goto error;
+ a0 *= (an * x) / (bn * cn * n);
+ sum += a0;
+ an += 1.0;
+ bn += 1.0;
+ cn += 1.0;
+ n += 1.0;
+ z = fabsf( a0 );
+ if( z > max )
+ max = z;
+ if( sum != 0 )
+ t = fabsf( a0 / sum );
+ else
+ t = z;
+ }
+while( t > MACHEPF );
+
+done:
+
+*err = fabsf( MACHEPF*max /sum );
+
+#if DEBUG
+ printf(" onef2f cancellation error %.5E\n", *err );
+#endif
+
+goto xit;
+
+error:
+#if DEBUG
+printf("onef2f does not converge\n");
+#endif
+*err = MAXNUMF;
+
+xit:
+
+#if DEBUG
+printf("onef2( %.2E %.2E %.2E %.5E ) = %.3E %.6E\n", a, b, c, x, n, sum);
+#endif
+return(sum);
+}
+
+
+
+float threef0f( float aa, float bb, float cc, float xx, float *err )
+{
+float a, b, c, x, n, a0, sum, t, conv, conv1;
+float an, bn, cn, max, z;
+
+a = aa;
+b = bb;
+c = cc;
+x = xx;
+an = a;
+bn = b;
+cn = c;
+a0 = 1.0;
+sum = 1.0;
+n = 1.0;
+t = 1.0;
+max = 0.0;
+conv = 1.0e38;
+conv1 = conv;
+
+do
+ {
+ if( an == 0.0 )
+ goto done;
+ if( bn == 0.0 )
+ goto done;
+ if( cn == 0.0 )
+ goto done;
+ if( (a0 > 1.0e34) || (n > 200) )
+ goto error;
+ a0 *= (an * bn * cn * x) / n;
+ an += 1.0;
+ bn += 1.0;
+ cn += 1.0;
+ n += 1.0;
+ z = fabsf( a0 );
+ if( z > max )
+ max = z;
+ if( z >= conv )
+ {
+ if( (z < max) && (z > conv1) )
+ goto done;
+ }
+ conv1 = conv;
+ conv = z;
+ sum += a0;
+ if( sum != 0 )
+ t = fabsf( a0 / sum );
+ else
+ t = z;
+ }
+while( t > MACHEPF );
+
+done:
+
+t = fabsf( MACHEPF*max/sum );
+#if DEBUG
+ printf(" threef0f cancellation error %.5E\n", t );
+#endif
+
+max = fabsf( conv/sum );
+if( max > t )
+ t = max;
+#if DEBUG
+ printf(" threef0f convergence %.5E\n", max );
+#endif
+
+goto xit;
+
+error:
+#if DEBUG
+printf("threef0f does not converge\n");
+#endif
+t = MAXNUMF;
+
+xit:
+
+#if DEBUG
+printf("threef0f( %.2E %.2E %.2E %.5E ) = %.3E %.6E\n", a, b, c, x, n, sum);
+#endif
+
+*err = t;
+return(sum);
+}
+
+
+
+
+float struvef( float vv, float xx )
+{
+float v, x, y, ya, f, g, h, t;
+float onef2err, threef0err;
+
+v = vv;
+x = xx;
+f = floorf(v);
+if( (v < 0) && ( v-f == 0.5 ) )
+ {
+ y = jvf( -v, x );
+ f = 1.0 - f;
+ g = 2.0 * floorf(0.5*f);
+ if( g != f )
+ y = -y;
+ return(y);
+ }
+t = 0.25*x*x;
+f = fabsf(x);
+g = 1.5 * fabsf(v);
+if( (f > 30.0) && (f > g) )
+ {
+ onef2err = MAXNUMF;
+ y = 0.0;
+ }
+else
+ {
+ y = onef2f( 1.0, 1.5, 1.5+v, -t, &onef2err );
+ }
+
+if( (f < 18.0) || (x < 0.0) )
+ {
+ threef0err = MAXNUMF;
+ ya = 0.0;
+ }
+else
+ {
+ ya = threef0f( 1.0, 0.5, 0.5-v, -1.0/t, &threef0err );
+ }
+
+f = sqrtf( PIF );
+h = powf( 0.5*x, v-1.0 );
+
+if( onef2err <= threef0err )
+ {
+ g = gammaf( v + 1.5 );
+ y = y * h * t / ( 0.5 * f * g );
+ return(y);
+ }
+else
+ {
+ g = gammaf( v + 0.5 );
+ ya = ya * h / ( f * g );
+ ya = ya + yvf( v, x );
+ return(ya);
+ }
+}
+
+
+
+
+/* Bessel function of noninteger order
+ */
+
+float yvf( float vv, float xx )
+{
+float v, x, y, t;
+int n;
+
+v = vv;
+x = xx;
+y = floorf( v );
+if( y == v )
+ {
+ n = v;
+ y = ynf( n, x );
+ return( y );
+ }
+t = PIF * v;
+y = (cosf(t) * jvf( v, x ) - jvf( -v, x ))/sinf(t);
+return( y );
+}
+
+/* Crossover points between ascending series and asymptotic series
+ * for Struve function
+ *
+ * v x
+ *
+ * 0 19.2
+ * 1 18.95
+ * 2 19.15
+ * 3 19.3
+ * 5 19.7
+ * 10 21.35
+ * 20 26.35
+ * 30 32.31
+ * 40 40.0
+ */