1 /* specfunc/zeta.c 2 * 3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman 4 * 5 * This program is free software; you can redistribute it and/or modify 6 * it under the terms of the GNU General Public License as published by 7 * the Free Software Foundation; either version 3 of the License, or (at 8 * your option) any later version. 9 * 10 * This program is distributed in the hope that it will be useful, but 11 * WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 13 * General Public License for more details. 14 * 15 * You should have received a copy of the GNU General Public License 16 * along with this program; if not, write to the Free Software 17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 18 */ 19 20 /* Author: G. Jungman */ 21 22 #include "gsl__config.h" 23 #include "gsl_math.h" 24 #include "gsl_errno.h" 25 #include "gsl_sf_elementary.h" 26 #include "gsl_sf_exp.h" 27 #include "gsl_sf_gamma.h" 28 #include "gsl_sf_pow_int.h" 29 #include "gsl_sf_zeta.h" 30 31 #include "gsl_specfunc__error.h" 32 33 #include "gsl_specfunc__chebyshev.h" 34 #include "gsl_specfunc__cheb_eval.c" 35 36 #define LogTwoPi_ 1.8378770664093454835606594728111235279723 37 38 39 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/ 40 41 /* chebyshev fit for (s(t)-1)Zeta[s(t)] 42 * s(t)= (t+1)/2 43 * -1 <= t <= 1 44 */ 45 static double zeta_xlt1_data[14] = { 46 1.48018677156931561235192914649, 47 0.25012062539889426471999938167, 48 0.00991137502135360774243761467, 49 -0.00012084759656676410329833091, 50 -4.7585866367662556504652535281e-06, 51 2.2229946694466391855561441361e-07, 52 -2.2237496498030257121309056582e-09, 53 -1.0173226513229028319420799028e-10, 54 4.3756643450424558284466248449e-12, 55 -6.2229632593100551465504090814e-14, 56 -6.6116201003272207115277520305e-16, 57 4.9477279533373912324518463830e-17, 58 -1.0429819093456189719660003522e-18, 59 6.9925216166580021051464412040e-21, 60 }; 61 static cheb_series zeta_xlt1_cs = { 62 zeta_xlt1_data, 63 13, 64 -1, 1, 65 8 66 }; 67 68 /* chebyshev fit for (s(t)-1)Zeta[s(t)] 69 * s(t)= (19t+21)/2 70 * -1 <= t <= 1 71 */ 72 static double zeta_xgt1_data[30] = { 73 19.3918515726724119415911269006, 74 9.1525329692510756181581271500, 75 0.2427897658867379985365270155, 76 -0.1339000688262027338316641329, 77 0.0577827064065028595578410202, 78 -0.0187625983754002298566409700, 79 0.0039403014258320354840823803, 80 -0.0000581508273158127963598882, 81 -0.0003756148907214820704594549, 82 0.0001892530548109214349092999, 83 -0.0000549032199695513496115090, 84 8.7086484008939038610413331863e-6, 85 6.4609477924811889068410083425e-7, 86 -9.6749773915059089205835337136e-7, 87 3.6585400766767257736982342461e-7, 88 -8.4592516427275164351876072573e-8, 89 9.9956786144497936572288988883e-9, 90 1.4260036420951118112457144842e-9, 91 -1.1761968823382879195380320948e-9, 92 3.7114575899785204664648987295e-10, 93 -7.4756855194210961661210215325e-11, 94 7.8536934209183700456512982968e-12, 95 9.9827182259685539619810406271e-13, 96 -7.5276687030192221587850302453e-13, 97 2.1955026393964279988917878654e-13, 98 -4.1934859852834647427576319246e-14, 99 4.6341149635933550715779074274e-15, 100 2.3742488509048340106830309402e-16, 101 -2.7276516388124786119323824391e-16, 102 7.8473570134636044722154797225e-17 103 }; 104 static cheb_series zeta_xgt1_cs = { 105 zeta_xgt1_data, 106 29, 107 -1, 1, 108 17 109 }; 110 111 112 /* chebyshev fit for Ln[Zeta[s(t)] - 1 - 2^(-s(t))] 113 * s(t)= 10 + 5t 114 * -1 <= t <= 1; 5 <= s <= 15 115 */ 116 static double zetam1_inter_data[24] = { 117 -21.7509435653088483422022339374, 118 -5.63036877698121782876372020472, 119 0.0528041358684229425504861579635, 120 -0.0156381809179670789342700883562, 121 0.00408218474372355881195080781927, 122 -0.0010264867349474874045036628282, 123 0.000260469880409886900143834962387, 124 -0.0000676175847209968878098566819447, 125 0.0000179284472587833525426660171124, 126 -4.83238651318556188834107605116e-6, 127 1.31913788964999288471371329447e-6, 128 -3.63760500656329972578222188542e-7, 129 1.01146847513194744989748396574e-7, 130 -2.83215225141806501619105289509e-8, 131 7.97733710252021423361012829496e-9, 132 -2.25850168553956886676250696891e-9, 133 6.42269392950164306086395744145e-10, 134 -1.83363861846127284505060843614e-10, 135 5.25309763895283179960368072104e-11, 136 -1.50958687042589821074710575446e-11, 137 4.34997545516049244697776942981e-12, 138 -1.25597782748190416118082322061e-12, 139 3.61280740072222650030134104162e-13, 140 -9.66437239205745207188920348801e-14 141 }; 142 static cheb_series zetam1_inter_cs = { 143 zetam1_inter_data, 144 22, 145 -1, 1, 146 12 147 }; 148 149 150 151 /* assumes s >= 0 and s != 1.0 */ 152 inline 153 static int 154 riemann_zeta_sgt0(double s, gsl_sf_result * result) 155 { 156 if(s < 1.0) { 157 gsl_sf_result c; 158 cheb_eval_e(&zeta_xlt1_cs, 2.0*s - 1.0, &c); 159 result->val = c.val / (s - 1.0); 160 result->err = c.err / fabs(s-1.0) + GSL_DBL_EPSILON * fabs(result->val); 161 return GSL_SUCCESS; 162 } 163 else if(s <= 20.0) { 164 double x = (2.0*s - 21.0)/19.0; 165 gsl_sf_result c; 166 cheb_eval_e(&zeta_xgt1_cs, x, &c); 167 result->val = c.val / (s - 1.0); 168 result->err = c.err / (s - 1.0) + GSL_DBL_EPSILON * fabs(result->val); 169 return GSL_SUCCESS; 170 } 171 else { 172 double f2 = 1.0 - pow(2.0,-s); 173 double f3 = 1.0 - pow(3.0,-s); 174 double f5 = 1.0 - pow(5.0,-s); 175 double f7 = 1.0 - pow(7.0,-s); 176 result->val = 1.0/(f2*f3*f5*f7); 177 result->err = 3.0 * GSL_DBL_EPSILON * fabs(result->val); 178 return GSL_SUCCESS; 179 } 180 } 181 182 183 inline 184 static int 185 riemann_zeta1ms_slt0(double s, gsl_sf_result * result) 186 { 187 if(s > -19.0) { 188 double x = (-19 - 2.0*s)/19.0; 189 gsl_sf_result c; 190 cheb_eval_e(&zeta_xgt1_cs, x, &c); 191 result->val = c.val / (-s); 192 result->err = c.err / (-s) + GSL_DBL_EPSILON * fabs(result->val); 193 return GSL_SUCCESS; 194 } 195 else { 196 double f2 = 1.0 - pow(2.0,-(1.0-s)); 197 double f3 = 1.0 - pow(3.0,-(1.0-s)); 198 double f5 = 1.0 - pow(5.0,-(1.0-s)); 199 double f7 = 1.0 - pow(7.0,-(1.0-s)); 200 result->val = 1.0/(f2*f3*f5*f7); 201 result->err = 3.0 * GSL_DBL_EPSILON * fabs(result->val); 202 return GSL_SUCCESS; 203 } 204 } 205 206 207 /* works for 5 < s < 15*/ 208 static int 209 riemann_zeta_minus_1_intermediate_s(double s, gsl_sf_result * result) 210 { 211 double t = (s - 10.0)/5.0; 212 gsl_sf_result c; 213 cheb_eval_e(&zetam1_inter_cs, t, &c); 214 result->val = exp(c.val) + pow(2.0, -s); 215 result->err = (c.err + 2.0*GSL_DBL_EPSILON)*result->val; 216 return GSL_SUCCESS; 217 } 218 219 220 /* assumes s is large and positive 221 * write: zeta(s) - 1 = zeta(s) * (1 - 1/zeta(s)) 222 * and expand a few terms of the product formula to evaluate 1 - 1/zeta(s) 223 * 224 * works well for s > 15 225 */ 226 static int 227 riemann_zeta_minus1_large_s(double s, gsl_sf_result * result) 228 { 229 double a = pow( 2.0,-s); 230 double b = pow( 3.0,-s); 231 double c = pow( 5.0,-s); 232 double d = pow( 7.0,-s); 233 double e = pow(11.0,-s); 234 double f = pow(13.0,-s); 235 double t1 = a + b + c + d + e + f; 236 double t2 = a*(b+c+d+e+f) + b*(c+d+e+f) + c*(d+e+f) + d*(e+f) + e*f; 237 /* 238 double t3 = a*(b*(c+d+e+f) + c*(d+e+f) + d*(e+f) + e*f) + 239 b*(c*(d+e+f) + d*(e+f) + e*f) + 240 c*(d*(e+f) + e*f) + 241 d*e*f; 242 double t4 = a*(b*(c*(d + e + f) + d*(e + f) + e*f) + c*(d*(e+f) + e*f) + d*e*f) + 243 b*(c*(d*(e+f) + e*f) + d*e*f) + 244 c*d*e*f; 245 double t5 = b*c*d*e*f + a*c*d*e*f+ a*b*d*e*f+ a*b*c*e*f+ a*b*c*d*f+ a*b*c*d*e; 246 double t6 = a*b*c*d*e*f; 247 */ 248 double numt = t1 - t2 /* + t3 - t4 + t5 - t6 */; 249 double zeta = 1.0/((1.0-a)*(1.0-b)*(1.0-c)*(1.0-d)*(1.0-e)*(1.0-f)); 250 result->val = numt*zeta; 251 result->err = (15.0/s + 1.0) * 6.0*GSL_DBL_EPSILON*result->val; 252 return GSL_SUCCESS; 253 } 254 255 256 #if 0 257 /* zeta(n) */ 258 #define ZETA_POS_TABLE_NMAX 100 259 static double zeta_pos_int_table_OLD[ZETA_POS_TABLE_NMAX+1] = { 260 -0.50000000000000000000000000000, /* zeta(0) */ 261 0.0 /* FIXME: DirectedInfinity() */, /* zeta(1) */ 262 1.64493406684822643647241516665, /* ... */ 263 1.20205690315959428539973816151, 264 1.08232323371113819151600369654, 265 1.03692775514336992633136548646, 266 1.01734306198444913971451792979, 267 1.00834927738192282683979754985, 268 1.00407735619794433937868523851, 269 1.00200839282608221441785276923, 270 1.00099457512781808533714595890, 271 1.00049418860411946455870228253, 272 1.00024608655330804829863799805, 273 1.00012271334757848914675183653, 274 1.00006124813505870482925854511, 275 1.00003058823630702049355172851, 276 1.00001528225940865187173257149, 277 1.00000763719763789976227360029, 278 1.00000381729326499983985646164, 279 1.00000190821271655393892565696, 280 1.00000095396203387279611315204, 281 1.00000047693298678780646311672, 282 1.00000023845050272773299000365, 283 1.00000011921992596531107306779, 284 1.00000005960818905125947961244, 285 1.00000002980350351465228018606, 286 1.00000001490155482836504123466, 287 1.00000000745071178983542949198, 288 1.00000000372533402478845705482, 289 1.00000000186265972351304900640, 290 1.00000000093132743241966818287, 291 1.00000000046566290650337840730, 292 1.00000000023283118336765054920, 293 1.00000000011641550172700519776, 294 1.00000000005820772087902700889, 295 1.00000000002910385044497099687, 296 1.00000000001455192189104198424, 297 1.00000000000727595983505748101, 298 1.00000000000363797954737865119, 299 1.00000000000181898965030706595, 300 1.00000000000090949478402638893, 301 1.00000000000045474737830421540, 302 1.00000000000022737368458246525, 303 1.00000000000011368684076802278, 304 1.00000000000005684341987627586, 305 1.00000000000002842170976889302, 306 1.00000000000001421085482803161, 307 1.00000000000000710542739521085, 308 1.00000000000000355271369133711, 309 1.00000000000000177635684357912, 310 1.00000000000000088817842109308, 311 1.00000000000000044408921031438, 312 1.00000000000000022204460507980, 313 1.00000000000000011102230251411, 314 1.00000000000000005551115124845, 315 1.00000000000000002775557562136, 316 1.00000000000000001387778780973, 317 1.00000000000000000693889390454, 318 1.00000000000000000346944695217, 319 1.00000000000000000173472347605, 320 1.00000000000000000086736173801, 321 1.00000000000000000043368086900, 322 1.00000000000000000021684043450, 323 1.00000000000000000010842021725, 324 1.00000000000000000005421010862, 325 1.00000000000000000002710505431, 326 1.00000000000000000001355252716, 327 1.00000000000000000000677626358, 328 1.00000000000000000000338813179, 329 1.00000000000000000000169406589, 330 1.00000000000000000000084703295, 331 1.00000000000000000000042351647, 332 1.00000000000000000000021175824, 333 1.00000000000000000000010587912, 334 1.00000000000000000000005293956, 335 1.00000000000000000000002646978, 336 1.00000000000000000000001323489, 337 1.00000000000000000000000661744, 338 1.00000000000000000000000330872, 339 1.00000000000000000000000165436, 340 1.00000000000000000000000082718, 341 1.00000000000000000000000041359, 342 1.00000000000000000000000020680, 343 1.00000000000000000000000010340, 344 1.00000000000000000000000005170, 345 1.00000000000000000000000002585, 346 1.00000000000000000000000001292, 347 1.00000000000000000000000000646, 348 1.00000000000000000000000000323, 349 1.00000000000000000000000000162, 350 1.00000000000000000000000000081, 351 1.00000000000000000000000000040, 352 1.00000000000000000000000000020, 353 1.00000000000000000000000000010, 354 1.00000000000000000000000000005, 355 1.00000000000000000000000000003, 356 1.00000000000000000000000000001, 357 1.00000000000000000000000000001, 358 1.00000000000000000000000000000, 359 1.00000000000000000000000000000, 360 1.00000000000000000000000000000 361 }; 362 #endif /* 0 */ 363 364 365 /* zeta(n) - 1 */ 366 #define ZETA_POS_TABLE_NMAX 100 367 static double zetam1_pos_int_table[ZETA_POS_TABLE_NMAX+1] = { 368 -1.5, /* zeta(0) */ 369 0.0, /* FIXME: Infinity */ /* zeta(1) - 1 */ 370 0.644934066848226436472415166646, /* zeta(2) - 1 */ 371 0.202056903159594285399738161511, 372 0.082323233711138191516003696541, 373 0.036927755143369926331365486457, 374 0.017343061984449139714517929790, 375 0.008349277381922826839797549849, 376 0.004077356197944339378685238508, 377 0.002008392826082214417852769232, 378 0.000994575127818085337145958900, 379 0.000494188604119464558702282526, 380 0.000246086553308048298637998047, 381 0.000122713347578489146751836526, 382 0.000061248135058704829258545105, 383 0.000030588236307020493551728510, 384 0.000015282259408651871732571487, 385 7.6371976378997622736002935630e-6, 386 3.8172932649998398564616446219e-6, 387 1.9082127165539389256569577951e-6, 388 9.5396203387279611315203868344e-7, 389 4.7693298678780646311671960437e-7, 390 2.3845050272773299000364818675e-7, 391 1.1921992596531107306778871888e-7, 392 5.9608189051259479612440207935e-8, 393 2.9803503514652280186063705069e-8, 394 1.4901554828365041234658506630e-8, 395 7.4507117898354294919810041706e-9, 396 3.7253340247884570548192040184e-9, 397 1.8626597235130490064039099454e-9, 398 9.3132743241966818287176473502e-10, 399 4.6566290650337840729892332512e-10, 400 2.3283118336765054920014559759e-10, 401 1.1641550172700519775929738354e-10, 402 5.8207720879027008892436859891e-11, 403 2.9103850444970996869294252278e-11, 404 1.4551921891041984235929632245e-11, 405 7.2759598350574810145208690123e-12, 406 3.6379795473786511902372363558e-12, 407 1.8189896503070659475848321007e-12, 408 9.0949478402638892825331183869e-13, 409 4.5474737830421540267991120294e-13, 410 2.2737368458246525152268215779e-13, 411 1.1368684076802278493491048380e-13, 412 5.6843419876275856092771829675e-14, 413 2.8421709768893018554550737049e-14, 414 1.4210854828031606769834307141e-14, 415 7.1054273952108527128773544799e-15, 416 3.5527136913371136732984695340e-15, 417 1.7763568435791203274733490144e-15, 418 8.8817842109308159030960913863e-16, 419 4.4408921031438133641977709402e-16, 420 2.2204460507980419839993200942e-16, 421 1.1102230251410661337205445699e-16, 422 5.5511151248454812437237365905e-17, 423 2.7755575621361241725816324538e-17, 424 1.3877787809725232762839094906e-17, 425 6.9388939045441536974460853262e-18, 426 3.4694469521659226247442714961e-18, 427 1.7347234760475765720489729699e-18, 428 8.6736173801199337283420550673e-19, 429 4.3368086900206504874970235659e-19, 430 2.1684043449972197850139101683e-19, 431 1.0842021724942414063012711165e-19, 432 5.4210108624566454109187004043e-20, 433 2.7105054312234688319546213119e-20, 434 1.3552527156101164581485233996e-20, 435 6.7762635780451890979952987415e-21, 436 3.3881317890207968180857031004e-21, 437 1.6940658945097991654064927471e-21, 438 8.4703294725469983482469926091e-22, 439 4.2351647362728333478622704833e-22, 440 2.1175823681361947318442094398e-22, 441 1.0587911840680233852265001539e-22, 442 5.2939559203398703238139123029e-23, 443 2.6469779601698529611341166842e-23, 444 1.3234889800848990803094510250e-23, 445 6.6174449004244040673552453323e-24, 446 3.3087224502121715889469563843e-24, 447 1.6543612251060756462299236771e-24, 448 8.2718061255303444036711056167e-25, 449 4.1359030627651609260093824555e-25, 450 2.0679515313825767043959679193e-25, 451 1.0339757656912870993284095591e-25, 452 5.1698788284564313204101332166e-26, 453 2.5849394142282142681277617708e-26, 454 1.2924697071141066700381126118e-26, 455 6.4623485355705318034380021611e-27, 456 3.2311742677852653861348141180e-27, 457 1.6155871338926325212060114057e-27, 458 8.0779356694631620331587381863e-28, 459 4.0389678347315808256222628129e-28, 460 2.0194839173657903491587626465e-28, 461 1.0097419586828951533619250700e-28, 462 5.0487097934144756960847711725e-29, 463 2.5243548967072378244674341938e-29, 464 1.2621774483536189043753999660e-29, 465 6.3108872417680944956826093943e-30, 466 3.1554436208840472391098412184e-30, 467 1.5777218104420236166444327830e-30, 468 7.8886090522101180735205378276e-31 469 }; 470 471 472 #define ZETA_NEG_TABLE_NMAX 99 473 #define ZETA_NEG_TABLE_SIZE 50 474 static double zeta_neg_int_table[ZETA_NEG_TABLE_SIZE] = { 475 -0.083333333333333333333333333333, /* zeta(-1) */ 476 0.008333333333333333333333333333, /* zeta(-3) */ 477 -0.003968253968253968253968253968, /* ... */ 478 0.004166666666666666666666666667, 479 -0.007575757575757575757575757576, 480 0.021092796092796092796092796093, 481 -0.083333333333333333333333333333, 482 0.44325980392156862745098039216, 483 -3.05395433027011974380395433027, 484 26.4562121212121212121212121212, 485 -281.460144927536231884057971014, 486 3607.5105463980463980463980464, 487 -54827.583333333333333333333333, 488 974936.82385057471264367816092, 489 -2.0052695796688078946143462272e+07, 490 4.7238486772162990196078431373e+08, 491 -1.2635724795916666666666666667e+10, 492 3.8087931125245368811553022079e+11, 493 -1.2850850499305083333333333333e+13, 494 4.8241448354850170371581670362e+14, 495 -2.0040310656516252738108421663e+16, 496 9.1677436031953307756992753623e+17, 497 -4.5979888343656503490437943262e+19, 498 2.5180471921451095697089023320e+21, 499 -1.5001733492153928733711440151e+23, 500 9.6899578874635940656497942895e+24, 501 -6.7645882379292820990945242302e+26, 502 5.0890659468662289689766332916e+28, 503 -4.1147288792557978697665486068e+30, 504 3.5666582095375556109684574609e+32, 505 -3.3066089876577576725680214670e+34, 506 3.2715634236478716264211227016e+36, 507 -3.4473782558278053878256455080e+38, 508 3.8614279832705258893092720200e+40, 509 -4.5892974432454332168863989006e+42, 510 5.7775386342770431824884825688e+44, 511 -7.6919858759507135167410075972e+46, 512 1.0813635449971654696354033351e+49, 513 -1.6029364522008965406067102346e+51, 514 2.5019479041560462843656661499e+53, 515 -4.1067052335810212479752045004e+55, 516 7.0798774408494580617452972433e+57, 517 -1.2804546887939508790190849756e+60, 518 2.4267340392333524078020892067e+62, 519 -4.8143218874045769355129570066e+64, 520 9.9875574175727530680652777408e+66, 521 -2.1645634868435185631335136160e+69, 522 4.8962327039620553206849224516e+71, /* ... */ 523 -1.1549023923963519663954271692e+74, /* zeta(-97) */ 524 2.8382249570693706959264156336e+76 /* zeta(-99) */ 525 }; 526 527 528 /* coefficients for Maclaurin summation in hzeta() 529 * B_{2j}/(2j)! 530 */ 531 static double hzeta_c[15] = { 532 1.00000000000000000000000000000, 533 0.083333333333333333333333333333, 534 -0.00138888888888888888888888888889, 535 0.000033068783068783068783068783069, 536 -8.2671957671957671957671957672e-07, 537 2.0876756987868098979210090321e-08, 538 -5.2841901386874931848476822022e-10, 539 1.3382536530684678832826980975e-11, 540 -3.3896802963225828668301953912e-13, 541 8.5860620562778445641359054504e-15, 542 -2.1748686985580618730415164239e-16, 543 5.5090028283602295152026526089e-18, 544 -1.3954464685812523340707686264e-19, 545 3.5347070396294674716932299778e-21, 546 -8.9535174270375468504026113181e-23 547 }; 548 549 #define ETA_POS_TABLE_NMAX 100 550 static double eta_pos_int_table[ETA_POS_TABLE_NMAX+1] = { 551 0.50000000000000000000000000000, /* eta(0) */ 552 M_LN2, /* eta(1) */ 553 0.82246703342411321823620758332, /* ... */ 554 0.90154267736969571404980362113, 555 0.94703282949724591757650323447, 556 0.97211977044690930593565514355, 557 0.98555109129743510409843924448, 558 0.99259381992283028267042571313, 559 0.99623300185264789922728926008, 560 0.99809429754160533076778303185, 561 0.99903950759827156563922184570, 562 0.99951714349806075414409417483, 563 0.99975768514385819085317967871, 564 0.99987854276326511549217499282, 565 0.99993917034597971817095419226, 566 0.99996955121309923808263293263, 567 0.99998476421490610644168277496, 568 0.99999237829204101197693787224, 569 0.99999618786961011347968922641, 570 0.99999809350817167510685649297, 571 0.99999904661158152211505084256, 572 0.99999952325821554281631666433, 573 0.99999976161323082254789720494, 574 0.99999988080131843950322382485, 575 0.99999994039889239462836140314, 576 0.99999997019885696283441513311, 577 0.99999998509923199656878766181, 578 0.99999999254955048496351585274, 579 0.99999999627475340010872752767, 580 0.99999999813736941811218674656, 581 0.99999999906868228145397862728, 582 0.99999999953434033145421751469, 583 0.99999999976716989595149082282, 584 0.99999999988358485804603047265, 585 0.99999999994179239904531592388, 586 0.99999999997089618952980952258, 587 0.99999999998544809143388476396, 588 0.99999999999272404460658475006, 589 0.99999999999636202193316875550, 590 0.99999999999818101084320873555, 591 0.99999999999909050538047887809, 592 0.99999999999954525267653087357, 593 0.99999999999977262633369589773, 594 0.99999999999988631316532476488, 595 0.99999999999994315658215465336, 596 0.99999999999997157829090808339, 597 0.99999999999998578914539762720, 598 0.99999999999999289457268000875, 599 0.99999999999999644728633373609, 600 0.99999999999999822364316477861, 601 0.99999999999999911182158169283, 602 0.99999999999999955591079061426, 603 0.99999999999999977795539522974, 604 0.99999999999999988897769758908, 605 0.99999999999999994448884878594, 606 0.99999999999999997224442439010, 607 0.99999999999999998612221219410, 608 0.99999999999999999306110609673, 609 0.99999999999999999653055304826, 610 0.99999999999999999826527652409, 611 0.99999999999999999913263826204, 612 0.99999999999999999956631913101, 613 0.99999999999999999978315956551, 614 0.99999999999999999989157978275, 615 0.99999999999999999994578989138, 616 0.99999999999999999997289494569, 617 0.99999999999999999998644747284, 618 0.99999999999999999999322373642, 619 0.99999999999999999999661186821, 620 0.99999999999999999999830593411, 621 0.99999999999999999999915296705, 622 0.99999999999999999999957648353, 623 0.99999999999999999999978824176, 624 0.99999999999999999999989412088, 625 0.99999999999999999999994706044, 626 0.99999999999999999999997353022, 627 0.99999999999999999999998676511, 628 0.99999999999999999999999338256, 629 0.99999999999999999999999669128, 630 0.99999999999999999999999834564, 631 0.99999999999999999999999917282, 632 0.99999999999999999999999958641, 633 0.99999999999999999999999979320, 634 0.99999999999999999999999989660, 635 0.99999999999999999999999994830, 636 0.99999999999999999999999997415, 637 0.99999999999999999999999998708, 638 0.99999999999999999999999999354, 639 0.99999999999999999999999999677, 640 0.99999999999999999999999999838, 641 0.99999999999999999999999999919, 642 0.99999999999999999999999999960, 643 0.99999999999999999999999999980, 644 0.99999999999999999999999999990, 645 0.99999999999999999999999999995, 646 0.99999999999999999999999999997, 647 0.99999999999999999999999999999, 648 0.99999999999999999999999999999, 649 1.00000000000000000000000000000, 650 1.00000000000000000000000000000, 651 1.00000000000000000000000000000, 652 }; 653 654 655 #define ETA_NEG_TABLE_NMAX 99 656 #define ETA_NEG_TABLE_SIZE 50 657 static double eta_neg_int_table[ETA_NEG_TABLE_SIZE] = { 658 0.25000000000000000000000000000, /* eta(-1) */ 659 -0.12500000000000000000000000000, /* eta(-3) */ 660 0.25000000000000000000000000000, /* ... */ 661 -1.06250000000000000000000000000, 662 7.75000000000000000000000000000, 663 -86.3750000000000000000000000000, 664 1365.25000000000000000000000000, 665 -29049.0312500000000000000000000, 666 800572.750000000000000000000000, 667 -2.7741322625000000000000000000e+7, 668 1.1805291302500000000000000000e+9, 669 -6.0523980051687500000000000000e+10, 670 3.6794167785377500000000000000e+12, 671 -2.6170760990658387500000000000e+14, 672 2.1531418140800295250000000000e+16, 673 -2.0288775575173015930156250000e+18, 674 2.1708009902623770590275000000e+20, 675 -2.6173826968455814932120125000e+22, 676 3.5324148876863877826668602500e+24, 677 -5.3042033406864906641493838981e+26, 678 8.8138218364311576767253114668e+28, 679 -1.6128065107490778547354654864e+31, 680 3.2355470001722734208527794569e+33, 681 -7.0876727476537493198506645215e+35, 682 1.6890450341293965779175629389e+38, 683 -4.3639690731216831157655651358e+40, 684 1.2185998827061261322605065672e+43, 685 -3.6670584803153006180101262324e+45, 686 1.1859898526302099104271449748e+48, 687 -4.1120769493584015047981746438e+50, 688 1.5249042436787620309090168687e+53, 689 -6.0349693196941307074572991901e+55, 690 2.5437161764210695823197691519e+58, 691 -1.1396923802632287851130360170e+61, 692 5.4180861064753979196802726455e+63, 693 -2.7283654799994373847287197104e+66, 694 1.4529750514918543238511171663e+69, 695 -8.1705519371067450079777183386e+71, 696 4.8445781606678367790247757259e+74, 697 -3.0246694206649519336179448018e+77, 698 1.9858807961690493054169047970e+80, 699 -1.3694474620720086994386818232e+83, 700 9.9070382984295807826303785989e+85, 701 -7.5103780796592645925968460677e+88, 702 5.9598418264260880840077992227e+91, 703 -4.9455988887500020399263196307e+94, 704 4.2873596927020241277675775935e+97, 705 -3.8791952037716162900707994047e+100, 706 3.6600317773156342245401829308e+103, 707 -3.5978775704117283875784869570e+106 /* eta(-99) */ 708 }; 709 710 711 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ 712 713 714 int gsl_sf_hzeta_e(const double s, const double q, gsl_sf_result * result) 715 { 716 /* CHECK_POINTER(result) */ 717 718 if(s <= 1.0 || q <= 0.0) { 719 DOMAIN_ERROR(result); 720 } 721 else { 722 const double max_bits = 54.0; 723 const double ln_term0 = -s * log(q); 724 725 if(ln_term0 < GSL_LOG_DBL_MIN + 1.0) { 726 UNDERFLOW_ERROR(result); 727 } 728 else if(ln_term0 > GSL_LOG_DBL_MAX - 1.0) { 729 OVERFLOW_ERROR (result); 730 } 731 else if((s > max_bits && q < 1.0) || (s > 0.5*max_bits && q < 0.25)) { 732 result->val = pow(q, -s); 733 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 734 return GSL_SUCCESS; 735 } 736 else if(s > 0.5*max_bits && q < 1.0) { 737 const double p1 = pow(q, -s); 738 const double p2 = pow(q/(1.0+q), s); 739 const double p3 = pow(q/(2.0+q), s); 740 result->val = p1 * (1.0 + p2 + p3); 741 result->err = GSL_DBL_EPSILON * (0.5*s + 2.0) * fabs(result->val); 742 return GSL_SUCCESS; 743 } 744 else { 745 /* Euler-Maclaurin summation formula 746 * [Moshier, p. 400, with several typo corrections] 747 */ 748 const int jmax = 12; 749 const int kmax = 10; 750 int j, k; 751 const double pmax = pow(kmax + q, -s); 752 double scp = s; 753 double pcp = pmax / (kmax + q); 754 double ans = pmax*((kmax+q)/(s-1.0) + 0.5); 755 756 for(k=0; k<kmax; k++) { 757 ans += pow(k + q, -s); 758 } 759 760 for(j=0; j<=jmax; j++) { 761 double delta = hzeta_c[j+1] * scp * pcp; 762 ans += delta; 763 if(fabs(delta/ans) < 0.5*GSL_DBL_EPSILON) break; 764 scp *= (s+2*j+1)*(s+2*j+2); 765 pcp /= (kmax + q)*(kmax + q); 766 } 767 768 result->val = ans; 769 result->err = 2.0 * (jmax + 1.0) * GSL_DBL_EPSILON * fabs(ans); 770 return GSL_SUCCESS; 771 } 772 } 773 } 774 775 776 int gsl_sf_zeta_e(const double s, gsl_sf_result * result) 777 { 778 /* CHECK_POINTER(result) */ 779 780 if(s == 1.0) { 781 DOMAIN_ERROR(result); 782 } 783 else if(s >= 0.0) { 784 return riemann_zeta_sgt0(s, result); 785 } 786 else { 787 /* reflection formula, [Abramowitz+Stegun, 23.2.5] */ 788 789 gsl_sf_result zeta_one_minus_s; 790 const int stat_zoms = riemann_zeta1ms_slt0(s, &zeta_one_minus_s); 791 const double sin_term = (fmod(s,2.0) == 0.0) ? 0.0 : sin(0.5*M_PI*fmod(s,4.0))/M_PI; 792 793 if(sin_term == 0.0) { 794 result->val = 0.0; 795 result->err = 0.0; 796 return GSL_SUCCESS; 797 } 798 else if(s > -170) { 799 /* We have to be careful about losing digits 800 * in calculating pow(2 Pi, s). The gamma 801 * function is fine because we were careful 802 * with that implementation. 803 * We keep an array of (2 Pi)^(10 n). 804 */ 805 const double twopi_pow[18] = { 1.0, 806 9.589560061550901348e+007, 807 9.195966217409212684e+015, 808 8.818527036583869903e+023, 809 8.456579467173150313e+031, 810 8.109487671573504384e+039, 811 7.776641909496069036e+047, 812 7.457457466828644277e+055, 813 7.151373628461452286e+063, 814 6.857852693272229709e+071, 815 6.576379029540265771e+079, 816 6.306458169130020789e+087, 817 6.047615938853066678e+095, 818 5.799397627482402614e+103, 819 5.561367186955830005e+111, 820 5.333106466365131227e+119, 821 5.114214477385391780e+127, 822 4.904306689854036836e+135 823 }; 824 const int n = floor((-s)/10.0); 825 const double fs = s + 10.0*n; 826 const double p = pow(2.0*M_PI, fs) / twopi_pow[n]; 827 828 gsl_sf_result g; 829 const int stat_g = gsl_sf_gamma_e(1.0-s, &g); 830 result->val = p * g.val * sin_term * zeta_one_minus_s.val; 831 result->err = fabs(p * g.val * sin_term) * zeta_one_minus_s.err; 832 result->err += fabs(p * sin_term * zeta_one_minus_s.val) * g.err; 833 result->err += GSL_DBL_EPSILON * (fabs(s)+2.0) * fabs(result->val); 834 return GSL_ERROR_SELECT_2(stat_g, stat_zoms); 835 } 836 else { 837 /* The actual zeta function may or may not 838 * overflow here. But we have no easy way 839 * to calculate it when the prefactor(s) 840 * overflow. Trying to use log's and exp 841 * is no good because we lose a couple 842 * digits to the exp error amplification. 843 * When we gather a little more patience, 844 * we can implement something here. Until 845 * then just give up. 846 */ 847 OVERFLOW_ERROR(result); 848 } 849 } 850 } 851 852 853 int gsl_sf_zeta_int_e(const int n, gsl_sf_result * result) 854 { 855 /* CHECK_POINTER(result) */ 856 857 if(n < 0) { 858 if(!GSL_IS_ODD(n)) { 859 result->val = 0.0; /* exactly zero at even negative integers */ 860 result->err = 0.0; 861 return GSL_SUCCESS; 862 } 863 else if(n > -ZETA_NEG_TABLE_NMAX) { 864 result->val = zeta_neg_int_table[-(n+1)/2]; 865 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 866 return GSL_SUCCESS; 867 } 868 else { 869 return gsl_sf_zeta_e((double)n, result); 870 } 871 } 872 else if(n == 1){ 873 DOMAIN_ERROR(result); 874 } 875 else if(n <= ZETA_POS_TABLE_NMAX){ 876 result->val = 1.0 + zetam1_pos_int_table[n]; 877 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 878 return GSL_SUCCESS; 879 } 880 else { 881 result->val = 1.0; 882 result->err = GSL_DBL_EPSILON; 883 return GSL_SUCCESS; 884 } 885 } 886 887 888 int gsl_sf_zetam1_e(const double s, gsl_sf_result * result) 889 { 890 if(s <= 5.0) 891 { 892 int stat = gsl_sf_zeta_e(s, result); 893 result->val = result->val - 1.0; 894 return stat; 895 } 896 else if(s < 15.0) 897 { 898 return riemann_zeta_minus_1_intermediate_s(s, result); 899 } 900 else 901 { 902 return riemann_zeta_minus1_large_s(s, result); 903 } 904 } 905 906 907 int gsl_sf_zetam1_int_e(const int n, gsl_sf_result * result) 908 { 909 if(n < 0) { 910 if(!GSL_IS_ODD(n)) { 911 result->val = -1.0; /* at even negative integers zetam1 == -1 since zeta is exactly zero */ 912 result->err = 0.0; 913 return GSL_SUCCESS; 914 } 915 else if(n > -ZETA_NEG_TABLE_NMAX) { 916 result->val = zeta_neg_int_table[-(n+1)/2] - 1.0; 917 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 918 return GSL_SUCCESS; 919 } 920 else { 921 /* could use gsl_sf_zetam1_e here but subtracting 1 makes no difference 922 for such large values, so go straight to the result */ 923 return gsl_sf_zeta_e((double)n, result); 924 } 925 } 926 else if(n == 1){ 927 DOMAIN_ERROR(result); 928 } 929 else if(n <= ZETA_POS_TABLE_NMAX){ 930 result->val = zetam1_pos_int_table[n]; 931 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 932 return GSL_SUCCESS; 933 } 934 else { 935 return gsl_sf_zetam1_e(n, result); 936 } 937 } 938 939 940 int gsl_sf_eta_int_e(int n, gsl_sf_result * result) 941 { 942 if(n > ETA_POS_TABLE_NMAX) { 943 result->val = 1.0; 944 result->err = GSL_DBL_EPSILON; 945 return GSL_SUCCESS; 946 } 947 else if(n >= 0) { 948 result->val = eta_pos_int_table[n]; 949 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 950 return GSL_SUCCESS; 951 } 952 else { 953 /* n < 0 */ 954 955 if(!GSL_IS_ODD(n)) { 956 /* exactly zero at even negative integers */ 957 result->val = 0.0; 958 result->err = 0.0; 959 return GSL_SUCCESS; 960 } 961 else if(n > -ETA_NEG_TABLE_NMAX) { 962 result->val = eta_neg_int_table[-(n+1)/2]; 963 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 964 return GSL_SUCCESS; 965 } 966 else { 967 gsl_sf_result z; 968 gsl_sf_result p; 969 int stat_z = gsl_sf_zeta_int_e(n, &z); 970 int stat_p = gsl_sf_exp_e((1.0-n)*M_LN2, &p); 971 int stat_m = gsl_sf_multiply_e(-p.val, z.val, result); 972 result->err = fabs(p.err * (M_LN2*(1.0-n)) * z.val) + z.err * fabs(p.val); 973 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); 974 return GSL_ERROR_SELECT_3(stat_m, stat_p, stat_z); 975 } 976 } 977 } 978 979 980 int gsl_sf_eta_e(const double s, gsl_sf_result * result) 981 { 982 /* CHECK_POINTER(result) */ 983 984 if(s > 100.0) { 985 result->val = 1.0; 986 result->err = GSL_DBL_EPSILON; 987 return GSL_SUCCESS; 988 } 989 else if(fabs(s-1.0) < 10.0*GSL_ROOT5_DBL_EPSILON) { 990 double del = s-1.0; 991 double c0 = M_LN2; 992 double c1 = M_LN2 * (M_EULER - 0.5*M_LN2); 993 double c2 = -0.0326862962794492996; 994 double c3 = 0.0015689917054155150; 995 double c4 = 0.00074987242112047532; 996 result->val = c0 + del * (c1 + del * (c2 + del * (c3 + del * c4))); 997 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 998 return GSL_SUCCESS; 999 } 1000 else { 1001 gsl_sf_result z; 1002 gsl_sf_result p; 1003 int stat_z = gsl_sf_zeta_e(s, &z); 1004 int stat_p = gsl_sf_exp_e((1.0-s)*M_LN2, &p); 1005 int stat_m = gsl_sf_multiply_e(1.0-p.val, z.val, result); 1006 result->err = fabs(p.err * (M_LN2*(1.0-s)) * z.val) + z.err * fabs(p.val); 1007 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); 1008 return GSL_ERROR_SELECT_3(stat_m, stat_p, stat_z); 1009 } 1010 } 1011 1012 1013 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ 1014 1015 #include "gsl_specfunc__eval.h" 1016 1017 double gsl_sf_zeta(const double s) 1018 { 1019 EVAL_RESULT(gsl_sf_zeta_e(s, &result)); 1020 } 1021 1022 double gsl_sf_hzeta(const double s, const double a) 1023 { 1024 EVAL_RESULT(gsl_sf_hzeta_e(s, a, &result)); 1025 } 1026 1027 double gsl_sf_zeta_int(const int s) 1028 { 1029 EVAL_RESULT(gsl_sf_zeta_int_e(s, &result)); 1030 } 1031 1032 double gsl_sf_zetam1(const double s) 1033 { 1034 EVAL_RESULT(gsl_sf_zetam1_e(s, &result)); 1035 } 1036 1037 double gsl_sf_zetam1_int(const int s) 1038 { 1039 EVAL_RESULT(gsl_sf_zetam1_int_e(s, &result)); 1040 } 1041 1042 double gsl_sf_eta_int(const int s) 1043 { 1044 EVAL_RESULT(gsl_sf_eta_int_e(s, &result)); 1045 } 1046 1047 double gsl_sf_eta(const double s) 1048 { 1049 EVAL_RESULT(gsl_sf_eta_e(s, &result)); 1050 } 1051