summaryrefslogtreecommitdiff
path: root/libm/float/struvef.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/struvef.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/struvef.c')
-rw-r--r--libm/float/struvef.c315
1 files changed, 0 insertions, 315 deletions
diff --git a/libm/float/struvef.c b/libm/float/struvef.c
deleted file mode 100644
index 4cf8854ed..000000000
--- a/libm/float/struvef.c
+++ /dev/null
@@ -1,315 +0,0 @@
-/* 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
- */