lib/libm/zetacf.c

   1  /*                                                     zetacf.c
   2  *
   3  *      Riemann zeta function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, zetacf();
  10  *
  11  * y = zetacf( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  *
  18  *
  19  *                inf.
  20  *                 -    -x
  21  *   zetac(x)  =   >   k   ,   x > 1,
  22  *                 -
  23  *                k=2
  24  *
  25  * is related to the Riemann zeta function by
  26  *
  27  *      Riemann zeta(x) = zetac(x) + 1.
  28  *
  29  * Extension of the function definition for x < 1 is implemented.
  30  * Zero is returned for x > log2(MAXNUM).
  31  *
  32  * An overflow error may occur for large negative x, due to the
  33  * gamma function in the reflection formula.
  34  *
  35  * ACCURACY:
  36  *
  37  * Tabulated values have full machine accuracy.
  38  *
  39  *                      Relative error:
  40  * arithmetic   domain     # trials      peak         rms
  41  *    IEEE      1,50        30000       5.5e-7      7.5e-8
  42  *
  43  *
  44  */
  45 \f
  46 /*
  47 Cephes Math Library Release 2.2:  July, 1992
  48 Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
  49 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  50 */
  51
  52 #include "mconf.h"
  53
  54
  55 /* Riemann zeta(x) - 1
  56  * for integer arguments between 0 and 30.
  57  */
  58 static float azetacf[] = {
  59 -1.50000000000000000000E0,
  60  1.70141183460469231730E38, /* infinity. */
  61  6.44934066848226436472E-1,
  62  2.02056903159594285400E-1,
  63  8.23232337111381915160E-2,
  64  3.69277551433699263314E-2,
  65  1.73430619844491397145E-2,
  66  8.34927738192282683980E-3,
  67  4.07735619794433937869E-3,
  68  2.00839282608221441785E-3,
  69  9.94575127818085337146E-4,
  70  4.94188604119464558702E-4,
  71  2.46086553308048298638E-4,
  72  1.22713347578489146752E-4,
  73  6.12481350587048292585E-5,
  74  3.05882363070204935517E-5,
  75  1.52822594086518717326E-5,
  76  7.63719763789976227360E-6,
  77  3.81729326499983985646E-6,
  78  1.90821271655393892566E-6,
  79  9.53962033872796113152E-7,
  80  4.76932986787806463117E-7,
  81  2.38450502727732990004E-7,
  82  1.19219925965311073068E-7,
  83  5.96081890512594796124E-8,
  84  2.98035035146522801861E-8,
  85  1.49015548283650412347E-8,
  86  7.45071178983542949198E-9,
  87  3.72533402478845705482E-9,
  88  1.86265972351304900640E-9,
  89  9.31327432419668182872E-10
  90 };
  91
  92
  93 /* 2**x (1 - 1/x) (zeta(x) - 1) = P(1/x)/Q(1/x), 1 <= x <= 10 */
  94 static float P[9] = {
  95   5.85746514569725319540E11,
  96   2.57534127756102572888E11,
  97   4.87781159567948256438E10,
  98   5.15399538023885770696E9,
  99   3.41646073514754094281E8,
 100   1.60837006880656492731E7,
 101   5.92785467342109522998E5,
 102   1.51129169964938823117E4,
 103   2.01822444485997955865E2,
 104 };
 105 static float Q[8] = {
 106 /*  1.00000000000000000000E0,*/
 107   3.90497676373371157516E11,
 108   5.22858235368272161797E10,
 109   5.64451517271280543351E9,
 110   3.39006746015350418834E8,
 111   1.79410371500126453702E7,
 112   5.66666825131384797029E5,
 113   1.60382976810944131506E4,
 114   1.96436237223387314144E2,
 115 };
 116
 117 /* log(zeta(x) - 1 - 2**-x), 10 <= x <= 50 */
 118 static float A[11] = {
 119  8.70728567484590192539E6,
 120  1.76506865670346462757E8,
 121  2.60889506707483264896E10,
 122  5.29806374009894791647E11,
 123  2.26888156119238241487E13,
 124  3.31884402932705083599E14,
 125  5.13778997975868230192E15,
 126 -1.98123688133907171455E15,
 127 -9.92763810039983572356E16,
 128  7.82905376180870586444E16,
 129  9.26786275768927717187E16,
 130 };
 131 static float B[10] = {
 132 /* 1.00000000000000000000E0,*/
 133 -7.92625410563741062861E6,
 134 -1.60529969932920229676E8,
 135 -2.37669260975543221788E10,
 136 -4.80319584350455169857E11,
 137 -2.07820961754173320170E13,
 138 -2.96075404507272223680E14,
 139 -4.86299103694609136686E15,
 140  5.34589509675789930199E15,
 141  5.71464111092297631292E16,
 142 -1.79915597658676556828E16,
 143 };
 144
 145 /* (1-x) (zeta(x) - 1), 0 <= x <= 1 */
 146
 147 static float R[6] = {
 148 -3.28717474506562731748E-1,
 149  1.55162528742623950834E1,
 150 -2.48762831680821954401E2,
 151  1.01050368053237678329E3,
 152  1.26726061410235149405E4,
 153 -1.11578094770515181334E5,
 154 };
 155 static float S[5] = {
 156 /* 1.00000000000000000000E0,*/
 157  1.95107674914060531512E1,
 158  3.17710311750646984099E2,
 159  3.03835500874445748734E3,
 160  2.03665876435770579345E4,
 161  7.43853965136767874343E4,
 162 };
 163
 164
 165 #define MAXL2 127
 166
 167 /*
 168  * Riemann zeta function, minus one
 169  */
 170
 171 extern float MACHEPF, PIO2F, MAXNUMF, PIF;
 172
 173 #ifndef ANSIC
 174 float sinf(), floorf(), gammaf(), powf(), expf();
 175 float polevlf(), p1evlf();
 176 #endif
 177
 178 #ifdef ANSIC
 179 float zetacf(float xx)
 180 #else
 181 float zetacf(xx)
 182 double xx;
 183 #endif
 184 {
 185 int i;
 186 float x, a, b, s, w;
 187
 188 x = xx;
 189 if( x < 0.0 )
 190         {
 191         if( x < -30.8148 )
 192                 {
 193                 mtherr( "zetacf", OVERFLOW );
 194                 return(0.0);
 195                 }
 196         s = 1.0 - x;
 197         w = zetacf( s );
 198         b = sinf(PIO2F*x) * powf(2.0*PIF, x) * gammaf(s) * (1.0 + w) / PIF;
 199         return(b - 1.0);
 200         }
 201
 202 if( x >= MAXL2 )
 203         return(0.0);    /* because first term is 2**-x */
 204
 205 /* Tabulated values for integer argument */
 206 w = floorf(x);
 207 if( w == x )
 208         {
 209         i = x;
 210         if( i < 31 )
 211                 {
 212                 return( azetacf[i] );
 213                 }
 214         }
 215
 216
 217 if( x < 1.0 )
 218         {
 219         w = 1.0 - x;
 220         a = polevlf( x, R, 5 ) / ( w * p1evlf( x, S, 5 ));
 221         return( a );
 222         }
 223
 224 if( x == 1.0 )
 225         {
 226         mtherr( "zetacf", SING );
 227         return( MAXNUMF );
 228         }
 229
 230 if( x <= 10.0 )
 231         {
 232         b = powf( 2.0, x ) * (x - 1.0);
 233         w = 1.0/x;
 234         s = (x * polevlf( w, P, 8 )) / (b * p1evlf( w, Q, 8 ));
 235         return( s );
 236         }
 237
 238 if( x <= 50.0 )
 239         {
 240         b = powf( 2.0, -x );
 241         w = polevlf( x, A, 10 ) / p1evlf( x, B, 10 );
 242         w = expf(w) + b;
 243         return(w);
 244         }
 245
 246
 247 /* Basic sum of inverse powers */
 248
 249
 250 s = 0.0;
 251 a = 1.0;
 252 do
 253         {
 254         a += 2.0;
 255         b = powf( a, -x );
 256         s += b;
 257         }
 258 while( b/s > MACHEPF );
 259
 260 b = powf( 2.0, -x );
 261 s = (s + b)/(1.0-b);
 262 return(s);
 263 }