MODULE Riemann_zeta ! From CERNLIB, with permission. ! Code converted using TO_F90 by Alan Miller ! Date: 2002-02-27 Time: 16:11:48 ! Author: K.S. Kolbig, 7 June 1992 ! Revision 1.1.1.1 1996/04/01 15:02:00 mclareni ! Mathlib gen/C (C315) ! The Riemann Zeta function ! Reference: ! Cody, W.J., Hillstrom, K.E. & Thather, H.C., `Chebyshev approximations ! for the Riemann zeta function', Math. Comp., vol.25 (1971), 537-547. ! N.B. The value returned by the CERNLIB routine was (zeta - 1) when the ! argument was > 1. This version returns the value of zeta. IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) CONTAINS FUNCTION r_zeta(x) RESULT(fn_val) REAL (dp), INTENT(IN) :: x REAL (dp) :: fn_val REAL (dp), PARAMETER :: delta = 1.0D-13, z1 = 1.0_dp, hf = z1/2, th = z1/3, & pi = 3.14159265358979324_dp, pih = pi/2, pi2 = 2*pi REAL (dp), PARAMETER :: p1(0:8) = (/ 1.28716812148244639D+10, & 1.37539693203702511D+10, 5.10665591836440610D+09, 8.56147100243331486D+08, & 7.48361812438023298D+07, 4.86010658546188251D+06, 2.73957499022140609D+05, & 4.63171084318342712D+03, 5.78758100409666066D+01 /) REAL (dp), PARAMETER :: q1(0:8) = (/ 2.57433624296484624D+10, & 5.93816564867959016D+09, 9.00633037326123344D+08, 8.04253663428328989D+07, & 5.60971175954192006D+06, 2.24743120289913752D+05, 7.57457890934153756D+03, & -2.37383578137377262D+01, 1.0_dp /) REAL (dp), PARAMETER :: p2(0:8) = (/ -6.88197293216348954D+06, & 7.48218916305315972D+06, -2.07584504810211014D+06, 3.55302557096214295D+05, & -4.06706449551854889D+04, 3.19804864027146911D+03, -1.69820937033722853D+02, & 5.61485842394289048D+00, -8.93888705926154944D-02 /) REAL (dp), PARAMETER :: q2(0:8) = (/ -1.29725624934891554D+09, & -9.48715407579907817D+08, -1.05496193474005203D+08, 4.67774488211993048D+06, & 3.12936040573813534D+06, 4.59581803839305070D+05, 3.88176109610396834D+04, & 1.92561544834491423D+03, 5.12578125000000000D+01 /) REAL (dp), PARAMETER :: p3(0:9) = (/ 1.66156480515774676D-11, & -4.68068827660654529D-09, 5.83519727319147047D-07, -4.17644012643145602D-05, & 1.85468422843597959D-03, -5.11288800220490241D-02, 8.10450231751100353D-01, & -5.69951948768478923D+00, 0.0_dp, 0.0_dp /) REAL (dp), PARAMETER :: q3(0:9) = (/ -6.99562633519191655D-10, & -1.77757961895149257D-08, -9.82231825734078036D-07, -2.84927282759096488D-05, & -5.81727909388048094D-04, -1.15848749169766586D-02, -1.28149124051978196D-01, & -1.11913057349097709D+00, -7.67928761604628813D-01, 1.0_dp /) REAL (dp), PARAMETER :: p4(0:8) = (/ 1.03144877188859712D-15, & -5.12584613964688241D-13, 1.12948794194873548D-10, -1.44234665373130952D-08, & 1.16824676984458098D-06, -6.14975167990314806D-05, 2.05594677988830328D-03, & -3.99339429394668869D-02, 3.45234976736178457D-01 /) REAL (dp), PARAMETER :: q4(0:8) = (/ 5.93959417288419050D-11, & -6.04755359079991806D-09, 3.64680208668388563D-07, -1.29456905568011812D-05, & 3.20189498470229250D-04, -5.07801557099994077D-03, 5.49628907881587266D-02, & -3.24517611155972419D-01, 1.0_dp /) REAL (dp) :: alfa, ap, aq, b0, b1, b2, f, h, t, v INTEGER :: ix, j v = x f = 1.0_dp IF (x /= 0 .AND. x < hf) THEN ix = x - delta IF (ABS(ix-x) <= delta) THEN IF (MOD(-ix,2) == 0) THEN h = 0.0_dp GO TO 70 ELSE v = 1.0_dp - x f = 2 * (-z1) ** ((1-ix)/2) * dgamma(v) / pi2 ** v END IF ELSE v = 1.0_dp - x f = 2 * SIN(pih*x) * dgamma(v) / pi2 ** v END IF END IF IF (x == 0) THEN h = -3 * hf ELSE IF (x == 1) THEN fn_val = 0.0_dp WRITE(*, *) 'Riemanns ZETA(X) HAS POLE AT X = 1' RETURN ELSE IF (v <= 5) THEN ap = p1(8) aq = q1(8) DO j = 7, 0, -1 ap = p1(j) + v * ap aq = q1(j) + v * aq END DO h = ap / (aq*(v-1)) - 1 ELSE IF (v <= 11) THEN t = th * (v-8) alfa = t + t b1 = 0 b2 = 0 DO j = 8, 0, -1 b0 = p2(j) + alfa * b1 - b2 b2 = b1 b1 = b0 END DO h = b0 - t * b2 b1 = 0 b2 = 0 DO j = 8, 0, -1 b0 = q2(j) + alfa * b1 - b2 b2 = b1 b1 = b0 END DO h = h / (b0 - t*b2) ELSE IF (v <= 25) THEN t = 1 / v ap = p3(7) DO j = 6, 0, -1 ap = p3(j) + t * ap END DO aq = q3(9) DO j = 8, 0, -1 aq = q3(j) + t * aq END DO h = hf ** (v - t*ap/aq) ELSE IF (v <= 55) THEN t = 1 / v ap = p4(8) aq = q4(8) DO j = 7, 0, -1 ap = p4(j) + t * ap aq = q4(j) + t * aq END DO h = hf ** (v-t*ap/aq) ELSE IF (v <= 90) THEN h = hf ** v + th ** v ELSE h = hf ** v END IF IF (x < 1) h = f * (1 + h) 70 IF (x > 1.0_dp) THEN fn_val = 1.0_dp + h ELSE fn_val = h END IF RETURN END FUNCTION r_zeta ! $Log: gamma64.F,v $ ! Revision 1.1.1.1 1996/04/01 15:01:54 mclareni ! Mathlib gen FUNCTION dgamma(x) RESULT(fn_val) REAL (dp), INTENT(IN) :: x REAL (dp) :: fn_val REAL (dp), PARAMETER :: c(0:15) = (/ 3.65738772508338244_dp, & 1.95754345666126827_dp, 0.33829711382616039_dp, 0.04208951276557549_dp, & 0.00428765048212909_dp, 0.00036521216929462_dp, 0.00002740064222642_dp, & 0.00000181240233365_dp, 0.00000010965775866_dp, 0.00000000598718405_dp, & 0.00000000030769081_dp, 0.00000000001431793_dp, 0.00000000000065109_dp, & 0.00000000000002596_dp, 0.00000000000000111_dp, 0.00000000000000004_dp /) REAL (dp) :: alfa, b0, b1, b2, f, h, u INTEGER :: i u = x IF (u <= 0) THEN WRITE(*, '(a, g13.5, a)') ' Argument value: ', u, ' <= 0 for routine DGAMMA' h = 0 RETURN END IF f = 1 IF (u < 3) THEN DO i = 1, INT(4-u) f = f / u u = u + 1 END DO ELSE DO i = 1, INT(u-3) u = u - 1 f = f * u END DO END IF h = u + u - 7 alfa = h + h b1 = 0 b2 = 0 DO i = 15, 0, -1 b0 = c(i) + alfa * b1 - b2 b2 = b1 b1 = b0 END DO fn_val = f * (b0 - h*b2) RETURN END FUNCTION dgamma END MODULE Riemann_zeta PROGRAM Test_Riemann_zeta USE Riemann_zeta IMPLICIT NONE ! Output the LHS of Table 23.3 (page 811) of Abramowitz & Stegun INTEGER :: n REAL (dp) :: zeta CHARACTER (LEN=21) :: text DO n = 2, 42 zeta = r_zeta( DBLE(n) ) WRITE(text, '(f18.16)') zeta text = text(:7) // ' ' // text(8:12) // ' ' // text(13:17) // ' ' // text(18:18) WRITE(*, '(i3, " ", a)') n, text END DO STOP END PROGRAM Test_Riemann_zeta