diff options
Diffstat (limited to 'libm/float/zetacf.c')
| -rw-r--r-- | libm/float/zetacf.c | 266 |
1 files changed, 266 insertions, 0 deletions
diff --git a/libm/float/zetacf.c b/libm/float/zetacf.c new file mode 100644 index 000000000..da2ace6a4 --- /dev/null +++ b/libm/float/zetacf.c @@ -0,0 +1,266 @@ + /* zetacf.c + * + * Riemann zeta function + * + * + * + * SYNOPSIS: + * + * float x, y, zetacf(); + * + * y = zetacf( x ); + * + * + * + * DESCRIPTION: + * + * + * + * inf. + * - -x + * zetac(x) = > k , x > 1, + * - + * k=2 + * + * is related to the Riemann zeta function by + * + * Riemann zeta(x) = zetac(x) + 1. + * + * Extension of the function definition for x < 1 is implemented. + * Zero is returned for x > log2(MAXNUM). + * + * An overflow error may occur for large negative x, due to the + * gamma function in the reflection formula. + * + * ACCURACY: + * + * Tabulated values have full machine accuracy. + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 1,50 30000 5.5e-7 7.5e-8 + * + * + */ + +/* +Cephes Math Library Release 2.2: July, 1992 +Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +#include <math.h> + + +/* Riemann zeta(x) - 1 + * for integer arguments between 0 and 30. + */ +static float azetacf[] = { +-1.50000000000000000000E0, + 1.70141183460469231730E38, /* infinity. */ + 6.44934066848226436472E-1, + 2.02056903159594285400E-1, + 8.23232337111381915160E-2, + 3.69277551433699263314E-2, + 1.73430619844491397145E-2, + 8.34927738192282683980E-3, + 4.07735619794433937869E-3, + 2.00839282608221441785E-3, + 9.94575127818085337146E-4, + 4.94188604119464558702E-4, + 2.46086553308048298638E-4, + 1.22713347578489146752E-4, + 6.12481350587048292585E-5, + 3.05882363070204935517E-5, + 1.52822594086518717326E-5, + 7.63719763789976227360E-6, + 3.81729326499983985646E-6, + 1.90821271655393892566E-6, + 9.53962033872796113152E-7, + 4.76932986787806463117E-7, + 2.38450502727732990004E-7, + 1.19219925965311073068E-7, + 5.96081890512594796124E-8, + 2.98035035146522801861E-8, + 1.49015548283650412347E-8, + 7.45071178983542949198E-9, + 3.72533402478845705482E-9, + 1.86265972351304900640E-9, + 9.31327432419668182872E-10 +}; + + +/* 2**x (1 - 1/x) (zeta(x) - 1) = P(1/x)/Q(1/x), 1 <= x <= 10 */ +static float P[9] = { + 5.85746514569725319540E11, + 2.57534127756102572888E11, + 4.87781159567948256438E10, + 5.15399538023885770696E9, + 3.41646073514754094281E8, + 1.60837006880656492731E7, + 5.92785467342109522998E5, + 1.51129169964938823117E4, + 2.01822444485997955865E2, +}; +static float Q[8] = { +/* 1.00000000000000000000E0,*/ + 3.90497676373371157516E11, + 5.22858235368272161797E10, + 5.64451517271280543351E9, + 3.39006746015350418834E8, + 1.79410371500126453702E7, + 5.66666825131384797029E5, + 1.60382976810944131506E4, + 1.96436237223387314144E2, +}; + +/* log(zeta(x) - 1 - 2**-x), 10 <= x <= 50 */ +static float A[11] = { + 8.70728567484590192539E6, + 1.76506865670346462757E8, + 2.60889506707483264896E10, + 5.29806374009894791647E11, + 2.26888156119238241487E13, + 3.31884402932705083599E14, + 5.13778997975868230192E15, +-1.98123688133907171455E15, +-9.92763810039983572356E16, + 7.82905376180870586444E16, + 9.26786275768927717187E16, +}; +static float B[10] = { +/* 1.00000000000000000000E0,*/ +-7.92625410563741062861E6, +-1.60529969932920229676E8, +-2.37669260975543221788E10, +-4.80319584350455169857E11, +-2.07820961754173320170E13, +-2.96075404507272223680E14, +-4.86299103694609136686E15, + 5.34589509675789930199E15, + 5.71464111092297631292E16, +-1.79915597658676556828E16, +}; + +/* (1-x) (zeta(x) - 1), 0 <= x <= 1 */ + +static float R[6] = { +-3.28717474506562731748E-1, + 1.55162528742623950834E1, +-2.48762831680821954401E2, + 1.01050368053237678329E3, + 1.26726061410235149405E4, +-1.11578094770515181334E5, +}; +static float S[5] = { +/* 1.00000000000000000000E0,*/ + 1.95107674914060531512E1, + 3.17710311750646984099E2, + 3.03835500874445748734E3, + 2.03665876435770579345E4, + 7.43853965136767874343E4, +}; + + +#define MAXL2 127 + +/* + * Riemann zeta function, minus one + */ + +extern float MACHEPF, PIO2F, MAXNUMF, PIF; + +#ifdef ANSIC +extern float sinf ( float xx ); +extern float floorf ( float x ); +extern float gammaf ( float xx ); +extern float powf ( float x, float y ); +extern float expf ( float xx ); +extern float polevlf ( float xx, float *coef, int N ); +extern float p1evlf ( float xx, float *coef, int N ); +#else +float sinf(), floorf(), gammaf(), powf(), expf(); +float polevlf(), p1evlf(); +#endif + +float zetacf(float xx) +{ +int i; +float x, a, b, s, w; + +x = xx; +if( x < 0.0 ) + { + if( x < -30.8148 ) + { + mtherr( "zetacf", OVERFLOW ); + return(0.0); + } + s = 1.0 - x; + w = zetacf( s ); + b = sinf(PIO2F*x) * powf(2.0*PIF, x) * gammaf(s) * (1.0 + w) / PIF; + return(b - 1.0); + } + +if( x >= MAXL2 ) + return(0.0); /* because first term is 2**-x */ + +/* Tabulated values for integer argument */ +w = floorf(x); +if( w == x ) + { + i = x; + if( i < 31 ) + { + return( azetacf[i] ); + } + } + + +if( x < 1.0 ) + { + w = 1.0 - x; + a = polevlf( x, R, 5 ) / ( w * p1evlf( x, S, 5 )); + return( a ); + } + +if( x == 1.0 ) + { + mtherr( "zetacf", SING ); + return( MAXNUMF ); + } + +if( x <= 10.0 ) + { + b = powf( 2.0, x ) * (x - 1.0); + w = 1.0/x; + s = (x * polevlf( w, P, 8 )) / (b * p1evlf( w, Q, 8 )); + return( s ); + } + +if( x <= 50.0 ) + { + b = powf( 2.0, -x ); + w = polevlf( x, A, 10 ) / p1evlf( x, B, 10 ); + w = expf(w) + b; + return(w); + } + + +/* Basic sum of inverse powers */ + + +s = 0.0; +a = 1.0; +do + { + a += 2.0; + b = powf( a, -x ); + s += b; + } +while( b/s > MACHEPF ); + +b = powf( 2.0, -x ); +s = (s + b)/(1.0-b); +return(s); +} |