diff options
Diffstat (limited to 'libm/float/struvef.c')
-rw-r--r-- | libm/float/struvef.c | 315 |
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 - */ |