htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/zeta.c

   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 <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_elementary.h>
  26 #include <gsl/gsl_sf_exp.h>
  27 #include <gsl/gsl_sf_gamma.h>
  28 #include <gsl/gsl_sf_pow_int.h>
  29 #include <gsl/gsl_sf_zeta.h>
  30
  31 #include "error.h"
  32
  33 #include "chebyshev.h"
  34 #include "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 loose 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 "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 }