diff options
Diffstat (limited to 'libm/double/struve.c')
| -rw-r--r-- | libm/double/struve.c | 312 |
1 files changed, 0 insertions, 312 deletions
diff --git a/libm/double/struve.c b/libm/double/struve.c deleted file mode 100644 index fabf0735e..000000000 --- a/libm/double/struve.c +++ /dev/null @@ -1,312 +0,0 @@ -/* struve.c - * - * Struve function - * - * - * - * SYNOPSIS: - * - * double v, x, y, struve(); - * - * y = struve( 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: - * - * Not accurately characterized, but spot checked against tables. - * - */ - - -/* -Cephes Math Library Release 2.81: June, 2000 -Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier -*/ -#include <math.h> -#define DEBUG 0 -#ifdef ANSIPROT -extern double gamma ( double ); -extern double pow ( double, double ); -extern double sqrt ( double ); -extern double yn ( int, double ); -extern double jv ( double, double ); -extern double fabs ( double ); -extern double floor ( double ); -extern double sin ( double ); -extern double cos ( double ); -double yv ( double, double ); -double onef2 (double, double, double, double, double * ); -double threef0 (double, double, double, double, double * ); -#else -double gamma(), pow(), sqrt(), yn(), yv(), jv(), fabs(), floor(); -double sin(), cos(); -double onef2(), threef0(); -#endif -static double stop = 1.37e-17; -extern double MACHEP; - -double onef2( a, b, c, x, err ) -double a, b, c, x; -double *err; -{ -double n, a0, sum, t; -double an, bn, cn, max, z; - -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 = fabs( a0 ); - if( z > max ) - max = z; - if( sum != 0 ) - t = fabs( a0 / sum ); - else - t = z; - } -while( t > stop ); - -done: - -*err = fabs( MACHEP*max /sum ); - -#if DEBUG - printf(" onef2 cancellation error %.5E\n", *err ); -#endif - -goto xit; - -error: -#if DEBUG -printf("onef2 does not converge\n"); -#endif -*err = 1.0e38; - -xit: - -#if DEBUG -printf("onef2( %.2E %.2E %.2E %.5E ) = %.3E %.6E\n", a, b, c, x, n, sum); -#endif -return(sum); -} - - - - -double threef0( a, b, c, x, err ) -double a, b, c, x; -double *err; -{ -double n, a0, sum, t, conv, conv1; -double an, bn, cn, max, z; - -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 = fabs( 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 = fabs( a0 / sum ); - else - t = z; - } -while( t > stop ); - -done: - -t = fabs( MACHEP*max/sum ); -#if DEBUG - printf(" threef0 cancellation error %.5E\n", t ); -#endif - -max = fabs( conv/sum ); -if( max > t ) - t = max; -#if DEBUG - printf(" threef0 convergence %.5E\n", max ); -#endif - -goto xit; - -error: -#if DEBUG -printf("threef0 does not converge\n"); -#endif -t = 1.0e38; - -xit: - -#if DEBUG -printf("threef0( %.2E %.2E %.2E %.5E ) = %.3E %.6E\n", a, b, c, x, n, sum); -#endif - -*err = t; -return(sum); -} - - - - -extern double PI; - -double struve( v, x ) -double v, x; -{ -double y, ya, f, g, h, t; -double onef2err, threef0err; - -f = floor(v); -if( (v < 0) && ( v-f == 0.5 ) ) - { - y = jv( -v, x ); - f = 1.0 - f; - g = 2.0 * floor(f/2.0); - if( g != f ) - y = -y; - return(y); - } -t = 0.25*x*x; -f = fabs(x); -g = 1.5 * fabs(v); -if( (f > 30.0) && (f > g) ) - { - onef2err = 1.0e38; - y = 0.0; - } -else - { - y = onef2( 1.0, 1.5, 1.5+v, -t, &onef2err ); - } - -if( (f < 18.0) || (x < 0.0) ) - { - threef0err = 1.0e38; - ya = 0.0; - } -else - { - ya = threef0( 1.0, 0.5, 0.5-v, -1.0/t, &threef0err ); - } - -f = sqrt( PI ); -h = pow( 0.5*x, v-1.0 ); - -if( onef2err <= threef0err ) - { - g = gamma( v + 1.5 ); - y = y * h * t / ( 0.5 * f * g ); - return(y); - } -else - { - g = gamma( v + 0.5 ); - ya = ya * h / ( f * g ); - ya = ya + yv( v, x ); - return(ya); - } -} - - - - -/* Bessel function of noninteger order - */ - -double yv( v, x ) -double v, x; -{ -double y, t; -int n; - -y = floor( v ); -if( y == v ) - { - n = v; - y = yn( n, x ); - return( y ); - } -t = PI * v; -y = (cos(t) * jv( v, x ) - jv( -v, x ))/sin(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 - */ |