;;; Science Collection
;;; special-functions/zeta-imp.rkt
;;; Copyright (c) 2004-2010 M. Douglas Williams
;;;
;;; This file is part of the Science Collection.
;;;
;;; The Science Collection is free software: you can redistribute it and/or
;;; modify it under the terms of the GNU Lesser General Public License as
;;; published by the Free Software Foundation, either version 3 of the License
;;; or (at your option) any later version.
;;;
;;; The Science Collection is distributed in the hope that it will be useful,
;;; but WITHOUT WARRANTY; without even the implied warranty of MERCHANTABILITY
;;; or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
;;; License for more details.
;;;
;;; You should have received a copy of the GNU Lesser General Public License
;;; along with the Science Collection.  If not, see
;;; <http://www.gnu.org/licenses/>.
;;;
;;; -----------------------------------------------------------------------------
;;;
;;; This file contains the routines that implement the zeta functions. This code
;;; is based on the Zeta Functions in the GNU Scientific Library (GSL).
;;;
;;; The gamma, psi, and zeta implementations are interdependent.  A
;;; single module, gamma-psi-zeta, includes the implementations of all
;;; all three sets of functions.
;;;
;;; Version  Date      Description
;;; 0.1.0    09/14/04  Initial implementation of the zeta functions. (MDW)
;;; 1.0.0    09/28/02  Marked as ready for Release 1.0. Added contracts for
;;;                    functions. (MDW)
;;; 2.0.0    11/17/07  Added unchecked versions of functions and getting ready
;;;                    for PLT Scheme V4.0. (MDW)
;;; 3.0.0    06/09/08  Changes required for V4.0. (MDW)
;;; 4.0.0    07/03/10  Changed the header and restructured the code. (MDW)

;;; Chebyshev fit for (s(t)-1_Zeta[s(t)]
;;; s(t) = (t+1)/2
;;; -1 <= t <= 1
(define zeta-xlt1-data
  #( 1.48018677156931561235192914649
     0.25012062539889426471999938167
     0.00991137302135360774243761467
    -0.00012084759656676410329833091
    -4.7585866367662556504652535281e-06
     2.2229946694466391855561441361e-07
    -2.2237496498030257121309056582e-09
    -1.0173226513229028319420799028e-10
     4.3756643450424558284466248449e-12
    -6.2229632593100551465504090814e-14
    -6.6116201003272207115277520305e-16
     4.9477279533373912324518463830e-17
    -1.0429819093456189719660003522e-18
     6.9925216166580021051464412040e-21))

(define zeta-xlt1-cs
  (make-chebyshev-series
   zeta-xlt1-data
   13 -1.0 1.0))

;;; Chebyshev fot for (s(t)-1)Zeta[s(t)]
;;; s(t) = (19t+21)/2
;;; -1 <= t <= 1
(define zeta-xgt1-data
  #(19.3918515726724119415911269006
     9.1525329692510756181581271500
     0.2427897658867379985365270155
    -0.1339000688262027338316641329
     0.0577827064065028595578410202
    -0.0187625983754002298566409700
     0.0039403014258320354840823803
    -0.0000581508273158127963598882
    -0.0003756148907214820704594549
     0.0001892530548109214349092999
    -0.0000549032199695513496115090
     8.7086484008939038610413331863e-06
     6.4609477924811889068410083425e-07
    -9.6749773915059089205835337136e-07
     3.6585400766767257736982342461e-07
    -8.4592516427275164351876072573e-08
     9.9956786144497936572288988883e-09
     1.4260036420951118112457144842e-09
    -1.1761968823382879195380320948e-09
     3.7114575899785204664648987295e-10
    -7.4756855194210961661210215325e-11
     7.8536934209183700456512982968e-12
     9.9827182259685539619810406271e-13
    -7.5276687030192221587850302453e-13
     2.1955026393964279988917878654e-13
    -4.1934859852834647427576319246e-14
     4.6341149635933550715779074274e-15
     2.3742488509048340106830309402e-16
    -2.7276516388124786119323824391e-16
     7.8473570134636044722154797225e-17))

(define zeta-xgt1-cs
  (make-chebyshev-series
   zeta-xgt1-data
   29 -1.0 1.0))

;;; chebyshev fit for Ln[Zeta[s(t)] - 1 - 2^(-s(t))]
;;; s(t)= 10 + 5t
;;; -1 <= t <= 1; 5 <= s <= 15
(define zetam1-inter-data
  #(-21.7509435653088483422022339374
     -5.63036877698121782876372020472
      0.0528041358684229425504861579635
     -0.0156381809179670789342700883562
      0.00408218474372355881195080781927
     -0.0010264867349474874045036628282
      0.000260469880409886900143834962387
     -0.0000676175847209968878098566819447
      0.0000179284472587833525426660171124
     -4.83238651318556188834107605116e-6
      1.31913788964999288471371329447e-6
     -3.63760500656329972578222188542e-7
      1.01146847513194744989748396574e-7
     -2.83215225141806501619105289509e-8
      7.97733710252021423361012829496e-9
     -2.25850168553956886676250696891e-9
      6.42269392950164306086395744145e-10
     -1.83363861846127284505060843614e-10
      5.25309763895283179960368072104e-11
     -1.50958687042589821074710575446e-11
      4.34997545516049244697776942981e-12
     -1.25597782748190416118082322061e-12
      3.61280740072222650030134104162e-13
     -9.66437239205745207188920348801e-14))

(define zetam1-inter-cs
  (make-chebyshev-series
   zetam1-inter-data
   23 -1.0 1.0))

;;; riemann-zeta-sgt0: real -> real
;;; Assumes s >- 0 and s != 1.0
(define (riemann-zeta-sgt0 s)
  (cond ((< s 1.0)
         (let ((val (unchecked-chebyshev-eval zeta-xlt1-cs 
                                    (- (* 2.0 s) 1.0))))
           (/ val (- s 1.0))))
        ((<= s 20.0)
         (let* ((x (/ (- (* 2.0 s) 21.0) 19.0))
                (val (unchecked-chebyshev-eval zeta-xgt1-cs x)))
           (/ val (- s 1.0))))
        (else
         (let ((f2 (- 1.0 (expt 2.0 (- s))))
               (f3 (- 1.0 (expt 3.0 (- s))))
               (f5 (- 1.0 (expt 5.0 (- s))))
               (f7 (- 1.0 (expt 7.0 (- s)))))
           (/ 1.0 (* f2 f3 f5 f7))))))

;;; reimann-zeta1m-slt0: real -> real
(define (riemann-zeta1m-slt0 s)
  (cond ((> s -19.0)
         (let* ((x (/ (- -19.0 (* 2.0 s)) 19.0))
                (val (unchecked-chebyshev-eval zeta-xgt1-cs x)))
           (/ val (- s))))
        (else
         (let ((f2 (- 1.0 (expt 2.0 (- (- 1.0 s)))))
               (f3 (- 1.0 (expt 3.0 (- (- 1.0 s)))))
               (f5 (- 1.0 (expt 5.0 (- (- 1.0 s)))))
               (f7 (- 1.0 (expt 7.0 (- (- 1.0 s))))))
           (/ 1.0 (* f2 f3 f5 f7))))))

;;; reimann-zeta-minus-1-intermediate-s: real -> real
;;; works for 5 < s < 15
(define (reimann-zeta-minus-1-intermediate-s s)
  (let* ((t (/ (- s 10.0) 5.0))
         (val (unchecked-chebyshev-eval zetam1-inter-cs t)))
    (+ (exp val) (expt 2.0 (- s)))))

;;; reimann-zeta-minus-1-large-s: real -> real
;;; assumes s is large and positive
;;; write: zeta(s) - 1 = zeta(s) * (1 - 1/zeta(s))
;;; and expand a few terms of the product formula to evaluate
;;; 1 - 1/zeta(s)
;;;
;;; works well for s > 15
(define (reimann-zeta-minus1-large-s s)
  (let* ((a (expt 2.0 (- s)))
         (b (expt 3.0 (- s)))
         (c (expt 5.0 (- s)))
         (d (expt 7.0 (- s)))
         (e (expt 11.0 (- s)))
         (f (expt 13.0 (- s)))
         (t1 (+ a b c d e f))
         (t2 (+ (* a (+ b c d e f))
                (* b (+ c d e f))
                (* c (+ d e f))
                (* d (+ e f))
                (* e f)))
         ;(t3 (+ (* a (+ (* b (+ c d e f))
         ;               (* c (+ d e f))
         ;               (* d (+ e f))
         ;               (* e f)))
         ;       (* b (+ (* c (+ d e f))
         ;               (* d (+ e f))
         ;               (* e f)))
         ;       (* c (+ (* d (+ e f))
         ;               (* e f)))
         ;       (* d e f)))
         ;(t4 (+ (* a (+ (* b (+ (* c (+ d e f))
         ;                       (* d (+ e f))
         ;                       (* e f)))
         ;               (* c (+ (* d (+ e f))
         ;                       (* e f)))
         ;               (* d e f)))
         ;       (* b (+ (* c (+ (* d (+ e f))
         ;                       (* e f)))
         ;               (* d e f)))
         ;       (* c d e f)))
         ;(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)))
         ;(t6 (* a b c d e f))
         (numt (- t1 t2)) ; + t3 - t4 + t5 - t6
         (zeta 
          (/ 1.0 
             (* (- 1.0 a)
                (- 1.0 b)
                (- 1.0 c)
                (- 1.0 d)
                (- 1.0 e)
                (- 1.0 f)))))
    (* numt zeta)))

;;; zeta(n) - 1
(define zeta-pos-table-nmax 100)
(define zetam1-pos-int-table
  #(-1.5                               ; zeta(0)
     0.0          ;; FIXME: Infinity   ; zeta(1) - 1
     0.644934066848226436472415166646  ; zeta(2) - 1
     0.202056903159594285399738161511
     0.082323233711138191516003696541
     0.036927755143369926331365486457
     0.017343061984449139714517929790
     0.008349277381922826839797549849
     0.004077356197944339378685238508
     0.002008392826082214417852769232
     0.000994575127818085337145958900
     0.000494188604119464558702282526
     0.000246086553308048298637998047
     0.000122713347578489146751836526
     0.000061248135058704829258545105
     0.000030588236307020493551728510
     0.000015282259408651871732571487
     7.6371976378997622736002935630e-6
     3.8172932649998398564616446219e-6
     1.9082127165539389256569577951e-6
     9.5396203387279611315203868344e-7
     4.7693298678780646311671960437e-7
     2.3845050272773299000364818675e-7
     1.1921992596531107306778871888e-7
     5.9608189051259479612440207935e-8
     2.9803503514652280186063705069e-8
     1.4901554828365041234658506630e-8
     7.4507117898354294919810041706e-9
     3.7253340247884570548192040184e-9
     1.8626597235130490064039099454e-9
     9.3132743241966818287176473502e-10
     4.6566290650337840729892332512e-10
     2.3283118336765054920014559759e-10
     1.1641550172700519775929738354e-10
     5.8207720879027008892436859891e-11
     2.9103850444970996869294252278e-11
     1.4551921891041984235929632245e-11
     7.2759598350574810145208690123e-12
     3.6379795473786511902372363558e-12
     1.8189896503070659475848321007e-12
     9.0949478402638892825331183869e-13
     4.5474737830421540267991120294e-13
     2.2737368458246525152268215779e-13
     1.1368684076802278493491048380e-13
     5.6843419876275856092771829675e-14
     2.8421709768893018554550737049e-14
     1.4210854828031606769834307141e-14
     7.1054273952108527128773544799e-15
     3.5527136913371136732984695340e-15
     1.7763568435791203274733490144e-15
     8.8817842109308159030960913863e-16
     4.4408921031438133641977709402e-16
     2.2204460507980419839993200942e-16
     1.1102230251410661337205445699e-16
     5.5511151248454812437237365905e-17
     2.7755575621361241725816324538e-17
     1.3877787809725232762839094906e-17
     6.9388939045441536974460853262e-18
     3.4694469521659226247442714961e-18
     1.7347234760475765720489729699e-18
     8.6736173801199337283420550673e-19
     4.3368086900206504874970235659e-19
     2.1684043449972197850139101683e-19
     1.0842021724942414063012711165e-19
     5.4210108624566454109187004043e-20
     2.7105054312234688319546213119e-20
     1.3552527156101164581485233996e-20
     6.7762635780451890979952987415e-21
     3.3881317890207968180857031004e-21
     1.6940658945097991654064927471e-21
     8.4703294725469983482469926091e-22
     4.2351647362728333478622704833e-22
     2.1175823681361947318442094398e-22
     1.0587911840680233852265001539e-22
     5.2939559203398703238139123029e-23
     2.6469779601698529611341166842e-23
     1.3234889800848990803094510250e-23
     6.6174449004244040673552453323e-24
     3.3087224502121715889469563843e-24
     1.6543612251060756462299236771e-24
     8.2718061255303444036711056167e-25
     4.1359030627651609260093824555e-25
     2.0679515313825767043959679193e-25
     1.0339757656912870993284095591e-25
     5.1698788284564313204101332166e-26
     2.5849394142282142681277617708e-26
     1.2924697071141066700381126118e-26
     6.4623485355705318034380021611e-27
     3.2311742677852653861348141180e-27
     1.6155871338926325212060114057e-27
     8.0779356694631620331587381863e-28
     4.0389678347315808256222628129e-28
     2.0194839173657903491587626465e-28
     1.0097419586828951533619250700e-28
     5.0487097934144756960847711725e-29
     2.5243548967072378244674341938e-29
     1.2621774483536189043753999660e-29
     6.3108872417680944956826093943e-30
     3.1554436208840472391098412184e-30
     1.5777218104420236166444327830e-30
     7.8886090522101180735205378276e-31))

(define zeta-neg-table-nmax 99)
(define zeta-neg-table-size 50)
(define zeta-neg-int-table
  #(-0.083333333333333333333333333333     ; zeta(-1)
     0.008333333333333333333333333333     ; zeta(-3)
    -0.003968253968253968253968253968     ; ...
     0.004166666666666666666666666667
    -0.007575757575757575757575757576
     0.021092796092796092796092796093
    -0.083333333333333333333333333333
     0.44325980392156862745098039216
    -3.05395433027011974380395433027
     26.4562121212121212121212121212
    -281.460144927536231884057971014
     3607.5105463980463980463980464
    -54827.583333333333333333333333
     974936.82385057471264367816092
    -2.0052695796688078946143462272e+07
     4.7238486772162990196078431373e+08
    -1.2635724795916666666666666667e+10
     3.8087931125245368811553022079e+11
    -1.2850850499305083333333333333e+13
     4.8241448354850170371581670362e+14
    -2.0040310656516252738108421663e+16
     9.1677436031953307756992753623e+17
    -4.5979888343656503490437943262e+19
     2.5180471921451095697089023320e+21
    -1.5001733492153928733711440151e+23
     9.6899578874635940656497942895e+24
    -6.7645882379292820990945242302e+26
     5.0890659468662289689766332916e+28
    -4.1147288792557978697665486068e+30
     3.5666582095375556109684574609e+32
    -3.3066089876577576725680214670e+34
     3.2715634236478716264211227016e+36
    -3.4473782558278053878256455080e+38
     3.8614279832705258893092720200e+40
    -4.5892974432454332168863989006e+42
     5.7775386342770431824884825688e+44
    -7.6919858759507135167410075972e+46
     1.0813635449971654696354033351e+49
    -1.6029364522008965406067102346e+51
     2.5019479041560462843656661499e+53
    -4.1067052335810212479752045004e+55
     7.0798774408494580617452972433e+57
    -1.2804546887939508790190849756e+60
     2.4267340392333524078020892067e+62
    -4.8143218874045769355129570066e+64
     9.9875574175727530680652777408e+66
    -2.1645634868435185631335136160e+69
     4.8962327039620553206849224516e+71    ; ...
    -1.1549023923963519663954271692e+74    ; zeta(-97)
     2.8382249570693706959264156336e+76))  ; zeta(-99)

;;; coefficients for Maclaurin summation in hzeta()
;;; B_{2j}/(2j)!
(define hzeta-c
  #( 1.00000000000000000000000000000
     0.083333333333333333333333333333
    -0.00138888888888888888888888888889
     0.000033068783068783068783068783069
    -8.2671957671957671957671957672e-07
     2.0876756987868098979210090321e-08
    -5.2841901386874931848476822022e-10
     1.3382536530684678832826980975e-11
    -3.3896802963225828668301953912e-13
     8.5860620562778445641359054504e-15
    -2.1748686985580618730415164239e-16
     5.5090028283602295152026526089e-18
    -1.3954464685812523340707686264e-19
     3.5347070396294674716932299778e-21
    -8.9535174270375468504026113181e-23))

(define eta-pos-table-nmax 100)
(define eta-pos-int-table
  `#( 0.50000000000000000000000000000  ; eta(0)
      ,ln2                             ; eta(1)
      0.82246703342411321823620758332  ; ...
      0.90154267736969571404980362113
      0.94703282949724591757650323447
      0.97211977044690930593565514355
      0.98555109129743510409843924448
      0.99259381992283028267042571313
      0.99623300185264789922728926008
      0.99809429754160533076778303185
      0.99903950759827156563922184570
      0.99951714349806075414409417483
      0.99975768514385819085317967871
      0.99987854276326511549217499282
      0.99993917034597971817095419226
      0.99996955121309923808263293263
      0.99998476421490610644168277496
      0.99999237829204101197693787224
      0.99999618786961011347968922641
      0.99999809350817167510685649297
      0.99999904661158152211505084256
      0.99999952325821554281631666433
      0.99999976161323082254789720494
      0.99999988080131843950322382485
      0.99999994039889239462836140314
      0.99999997019885696283441513311
      0.99999998509923199656878766181
      0.99999999254955048496351585274
      0.99999999627475340010872752767
      0.99999999813736941811218674656
      0.99999999906868228145397862728
      0.99999999953434033145421751469
      0.99999999976716989595149082282
      0.99999999988358485804603047265
      0.99999999994179239904531592388
      0.99999999997089618952980952258
      0.99999999998544809143388476396
      0.99999999999272404460658475006
      0.99999999999636202193316875550
      0.99999999999818101084320873555
      0.99999999999909050538047887809
      0.99999999999954525267653087357
      0.99999999999977262633369589773
      0.99999999999988631316532476488
      0.99999999999994315658215465336
      0.99999999999997157829090808339
      0.99999999999998578914539762720
      0.99999999999999289457268000875
      0.99999999999999644728633373609
      0.99999999999999822364316477861
      0.99999999999999911182158169283
      0.99999999999999955591079061426
      0.99999999999999977795539522974
      0.99999999999999988897769758908
      0.99999999999999994448884878594
      0.99999999999999997224442439010
      0.99999999999999998612221219410
      0.99999999999999999306110609673
      0.99999999999999999653055304826
      0.99999999999999999826527652409
      0.99999999999999999913263826204
      0.99999999999999999956631913101
      0.99999999999999999978315956551
      0.99999999999999999989157978275
      0.99999999999999999994578989138
      0.99999999999999999997289494569
      0.99999999999999999998644747284
      0.99999999999999999999322373642
      0.99999999999999999999661186821
      0.99999999999999999999830593411
      0.99999999999999999999915296705
      0.99999999999999999999957648353
      0.99999999999999999999978824176
      0.99999999999999999999989412088
      0.99999999999999999999994706044
      0.99999999999999999999997353022
      0.99999999999999999999998676511
      0.99999999999999999999999338256
      0.99999999999999999999999669128
      0.99999999999999999999999834564
      0.99999999999999999999999917282
      0.99999999999999999999999958641
      0.99999999999999999999999979320
      0.99999999999999999999999989660
      0.99999999999999999999999994830
      0.99999999999999999999999997415
      0.99999999999999999999999998708
      0.99999999999999999999999999354
      0.99999999999999999999999999677
      0.99999999999999999999999999838
      0.99999999999999999999999999919
      0.99999999999999999999999999960
      0.99999999999999999999999999980
      0.99999999999999999999999999990
      0.99999999999999999999999999995
      0.99999999999999999999999999997
      0.99999999999999999999999999999
      0.99999999999999999999999999999
      1.00000000000000000000000000000
      1.00000000000000000000000000000
      1.00000000000000000000000000000))

(define eta-neg-table-nmax 99)
(define eta-neg-table-size 50)
(define eta-neg-int-table
  #( 0.25000000000000000000000000000   ; eta(-1)
    -0.12500000000000000000000000000   ; eta(-3)
     0.25000000000000000000000000000   ; ...
    -1.06250000000000000000000000000
     7.75000000000000000000000000000
    -86.3750000000000000000000000000
     1365.25000000000000000000000000
    -29049.0312500000000000000000000
     800572.750000000000000000000000
    -2.7741322625000000000000000000e+7
     1.1805291302500000000000000000e+9
    -6.0523980051687500000000000000e+10
     3.6794167785377500000000000000e+12
    -2.6170760990658387500000000000e+14
     2.1531418140800295250000000000e+16
    -2.0288775575173015930156250000e+18
     2.1708009902623770590275000000e+20
    -2.6173826968455814932120125000e+22
     3.5324148876863877826668602500e+24
    -5.3042033406864906641493838981e+26
     8.8138218364311576767253114668e+28
    -1.6128065107490778547354654864e+31
     3.2355470001722734208527794569e+33
    -7.0876727476537493198506645215e+35
     1.6890450341293965779175629389e+38
    -4.3639690731216831157655651358e+40
     1.2185998827061261322605065672e+43
    -3.6670584803153006180101262324e+45
     1.1859898526302099104271449748e+48
    -4.1120769493584015047981746438e+50
     1.5249042436787620309090168687e+53
    -6.0349693196941307074572991901e+55
     2.5437161764210695823197691519e+58
    -1.1396923802632287851130360170e+61
     5.4180861064753979196802726455e+63
    -2.7283654799994373847287197104e+66
     1.4529750514918543238511171663e+69
    -8.1705519371067450079777183386e+71
     4.8445781606678367790247757259e+74
    -3.0246694206649519336179448018e+77
     1.9858807961690493054169047970e+80
    -1.3694474620720086994386818232e+83
     9.9070382984295807826303785989e+85
    -7.5103780796592645925968460677e+88
     5.9598418264260880840077992227e+91
    -4.9455988887500020399263196307e+94
     4.2873596927020241277675775935e+97
    -3.8791952037716162900707994047e+100
     3.6600317773156342245401829308e+103
    -3.5978775704117283875784869570e+106)) ; eta(-99)

;;; hzeta: real x real -> real
;;; This routine computes the Hurwitz zeta function zeta(s,q)
;;; for s > 1, q > 0.
(define (hzeta s q)
  (if (or (<= s 1.0)
          (<= q 0.0))
      +nan.0
      (let ((max-bits 54)
            (ln-term0 (* (- s) (log q))))
        (cond ((< ln-term0 (+ log-double-min 1))
               -inf.0)
              ((> ln-term0 (- log-double-max 1))
               +inf.0)
              ((or (and (> s max-bits)
                        (< q 1.0))
                   (and (> s (* 0.5 max-bits))
                        (< q 0.25)))
               (expt q (- s)))
              ((and (> s (* 0.5 max-bits))
                    (< q 1.0))
               (let ((p1 (expt q (- s)))
                     (p2 (expt (/ q (+ 1.0 q)) s))
                     (p3 (expt (/ q (+ 2.0 q)) s)))
                 (* p1 (+ 1.0 p2 p2))))
              (else
               ;; Euler-Maclaurin summation formula
               ;; [Moshier, p. 400, with several typo corrections
               (let* ((jmax 12)
                      (kmax 10)
                      (pmax (expt (+ kmax q) (- s)))
                      (scp s)
                      (pcp (/ pmax (+ kmax q)))
                      (ans (* pmax (+ (/ (+ kmax q) (- s 1.0)) 0.5))))
                 (do ((k 0 (+ k 1)))
                   ((= k kmax) (void))
                   (set! ans (+ ans (expt (+ k q) (- s)))))
                 (let/ec break
                   (do ((j 0 (+ j 1)))
                     ((> j jmax) (void))
                     (let ((delta (* (vector-ref hzeta-c (+ j 1))
                                     scp pcp)))
                       (set! ans (+ ans delta))
                       (when (< (abs (/ delta ans))
                                (* 0.5 double-epsilon))
                         (break))
                       (set! scp (* scp 
                                    (+ s (* 2 j) 1.0) (+ s (* 2 j) 2.0)))
                       (set! pcp (/ pcp
                                    (* (+ kmax q) (+ kmax q)))))))
                 ans))))))

;;; An array of (2 pi)^(10 n)
(define twopi-pow
  #(1.0
    9.589560061550901348e+007
    9.195966217409212684e+015
    8.818527036583869903e+023
    8.456579467173150313e+031
    8.109487671573504384e+039
    7.776641909496069036e+047
    7.457457466828644277e+055
    7.151373628461452286e+063
    6.857852693272229709e+071
    6.576379029540265771e+079
    6.306458169130020789e+087
    6.047615938853066678e+095
    5.799397627482402614e+103
    5.561367186955830005e+111
    5.333106466365131227e+119
    5.114214477385391780e+127
    4.904306689854036836e+135))

;;; zeta: real -> real
;;; This routine computes the Riemann zeta function zeta(s)
;;; for arbitrary s, s /= 0.
(define (zeta s)
  (cond ((= s 1.0)
         +nan.0)
        ((>= s 0.0)
         (riemann-zeta-sgt0 s))
        (else
         ;; reflection formula, [Abramowitz & Stegun, 23.2.5]
         (let ((zeta-one-minus-s (riemann-zeta1m-slt0 s))
               (sin-term (/ (sin (* 0.5 pi s)) pi)))
           (cond ((= sin-term 0.0)
                  0.0)
                 ((> s -170.0)
                  ;; We have to be careful about losing digits
                  ;; in calculating pow(2 pi, s).  The gamma
                  ;; function is fine because we were careful
                  ;; with the implementation.
                  (let* ((n (inexact->exact
                             (floor (/ (- s) 10.0))))
                         (fs (+ s (* 10.0 n)))
                         (p (/ (expt (* 2.0 pi) fs)
                               (vector-ref twopi-pow n)))
                         (val (gamma (- 1.0 s))))
                    (* p val sin-term zeta-one-minus-s)))
                 (else
                  ;; The actual zeta function may or may not
                  ;; overflow here.  But, we have no easy way
                  ;; to calculate it when the prefactor(s)
                  ;; overflow.  Trying to use log's and exp
                  ;; is no good because we loose a couple
                  ;; digits to the exp error amplification.
                  ;; When we gather a little more patience
                  ;; we can implement something here.  Until
                  ;; then, just give up.
                  +inf.0))))))

;;; zeta-int: integer -> real
;;; This routine computes the Riemann zeta function zeta(n)
;;; for integer n, n /= 1.
(define (zeta-int n)
  (cond ((< n 0)
         (cond ((even? n) 0)
               ((> n (- zeta-neg-table-nmax))
                (vector-ref zeta-neg-int-table (quotient (- (+ n 1)) 2)))
               (else
                (zeta n))))
        ((= n 1)
         +nan.0)
        ((<= n zeta-pos-table-nmax)
         (+ 1 (vector-ref zetam1-pos-int-table (inexact->exact n))))
        (else
         1.0)))

;;; zetam1: real -> real
;;; This routine computes zeta(s) - 1 for arbitrary s, s /= 1.
(define (zetam1 s)
  (cond ((<= s 5.0)
         (- (zeta s) 1.0))
        ((< s 15.0)
         (reimann-zeta-minus-1-intermediate-s s))
        (else
         (reimann-zeta-minus1-large-s s))))

;;; zetam1-int: integer -> real
;;; This routine computes zeta(n) - 1 for integer n, n /= 1.
(define (zetam1-int n)
  (cond ((< n 0)
         (cond ((even? n)
                0.0)
               ((> n (- zeta-neg-table-nmax))
                (- (vector-ref zeta-neg-int-table (quotient (- (+ n 1)) 2))
                   1.0))
               (else
                (zeta n))))
        ((= n 1)
         +nan.0)
        ((<= n zeta-pos-table-nmax)
         (vector-ref zetam1-pos-int-table (inexact->exact n)))
        (else
         (zetam1 n))))

;;; eta-int: integer -> real
;;; This routine computes the eta function eta(n) for integer n.
(define (eta-int n)
  (cond ((> n eta-pos-table-nmax)
         1.0)
        ((>= n 0)
         (vector-ref eta-pos-int-table (inexact->exact n)))
        (else
         ;; n < 0
         (cond ((even? n)
                0.0)
               ((> n (- eta-neg-table-nmax))
                (vector-ref eta-neg-int-table
                            (quotient (- (+ n 1)) 2)))
               (else
                (let* ((z-val (zeta-int n))
                       (p-val (exp (* (- 1.0 n) ln2)))
                       (m-val (* (- p-val) z-val)))
                  m-val))))))

;;; eta: real -> real
;;; This routine computes the eta function eta(s) for arbitrary s.
(define (eta s)
  (cond ((> s 100.0)
         1.0)
        ((< (abs (- s 1.0)) (* 10.0 root5-double-epsilon))
         (let ((del (- s 1.0))
               (c0 ln2)
               (c1 (* ln2 (- euler (* 0.5 ln2))))
               (c2 -0.0326862962794492996)
               (c3  0.0015689917054155150)
               (c4  0.00074987242112047532))
           (+ c0 (* del 
                    (+ c1 (* del 
                             (+ c2 (* del 
                                      (+ c3 (* del c4))))))))))
        (else
         (let* ((z-val (zeta s))
                (p-val (exp (* (- 1.0 s) ln2)))
                (m-val (* (- 1.0 p-val) z-val)))
           m-val))))