porting some code for complex lgamma, beta, eta, and zeta functions

Author: JeffBezanson

examples/specfun.j

@@ -0,0 +1,165 @@
+function angle_restrict_symm(theta)
+    P1 = 4 * 7.8539812564849853515625e-01
+    P2 = 4 * 3.7748947079307981766760e-08
+    P3 = 4 * 2.6951514290790594840552e-15
+
+    y = 2*floor(theta/(2*pi))
+    r = ((theta - y*P1) - y*P2) - y*P3
+    if (r > pi)
+        r -= (2*pi)
+    end
+    return r
+end
+
+const clg_coeff = [76.18009172947146,
+                   -86.50532032941677,
+                   24.01409824083091,
+                   -1.231739572450155,
+                   0.1208650973866179e-2,
+                   -0.5395239384953e-5]
+
+function clgamma_lanczos(z)
+    sqrt2pi = 2.5066282746310005
+    
+    y = x = z
+    temp = x + 5.5
+    zz = log(temp)
+    zz = zz * (x+0.5)
+    temp -= zz
+    ser = complex(1.000000000190015, 0)
+    for j=1:6
+        y += 1.0
+        zz = clg_coeff[j]/y
+        ser += zz
+    end
+    zz = sqrt2pi*ser / x
+    return log(zz) - temp
+end
+
+function lgamma(z::Complex)
+    if real(z) <= 0.5
+        a = clgamma_lanczos(1-z)
+        b = log(sin(pi * z))
+        logpi = 1.14472988584940017
+        z = logpi - b - a
+    else
+        z = clgamma_lanczos(z)
+    end
+    complex(real(z), angle_restrict_symm(imag(z)))
+end
+
+beta(x::Number, w::Number) = exp(lgamma(x)+lgamma(w)-lgamma(x+w))
+
+const eta_coeffs =
+    [.99999999999999999997,
+     -.99999999999999999821,
+     .99999999999999994183,
+     -.99999999999999875788,
+     .99999999999998040668,
+     -.99999999999975652196,
+     .99999999999751767484,
+     -.99999999997864739190,
+     .99999999984183784058,
+     -.99999999897537734890,
+     .99999999412319859549,
+     -.99999996986230482845,
+     .99999986068828287678,
+     -.99999941559419338151,
+     .99999776238757525623,
+     -.99999214148507363026,
+     .99997457616475604912,
+     -.99992394671207596228,
+     .99978893483826239739,
+     -.99945495809777621055,
+     .99868681159465798081,
+     -.99704078337369034566,
+     .99374872693175507536,
+     -.98759401271422391785,
+     .97682326283354439220,
+     -.95915923302922997013,
+     .93198380256105393618,
+     -.89273040299591077603,
+     .83945793215750220154,
+     -.77148960729470505477,
+     .68992761745934847866,
+     -.59784149990330073143,
+     .50000000000000000000,
+     -.40215850009669926857,
+     .31007238254065152134,
+     -.22851039270529494523,
+     .16054206784249779846,
+     -.10726959700408922397,
+     .68016197438946063823e-1,
+     -.40840766970770029873e-1,
+     .23176737166455607805e-1,
+     -.12405987285776082154e-1,
+     .62512730682449246388e-2,
+     -.29592166263096543401e-2,
+     .13131884053420191908e-2,
+     -.54504190222378945440e-3,
+     .21106516173760261250e-3,
+     -.76053287924037718971e-4,
+     .25423835243950883896e-4,
+     -.78585149263697370338e-5,
+     .22376124247437700378e-5,
+     -.58440580661848562719e-6,
+     .13931171712321674741e-6,
+     -.30137695171547022183e-7,
+     .58768014045093054654e-8,
+     -.10246226511017621219e-8,
+     .15816215942184366772e-9,
+     -.21352608103961806529e-10,
+     .24823251635643084345e-11,
+     -.24347803504257137241e-12,
+     .19593322190397666205e-13,
+     -.12421162189080181548e-14,
+     .58167446553847312884e-16,
+     -.17889335846010823161e-17,
+     .27105054312137610850e-19]
+
+function eta(z)
+    if z == 0
+        return complex(0.5)
+    end
+    re, im = reim(z)
+    if im==0 && re < 0 && integer_valued(re) && re==round(re/2)*2
+        return complex(0.0)
+    end
+    reflect = false
+    if re < 0.5
+        re = 1-re
+        im = -im
+        reflect = true
+    end
+    dn = float64(length(eta_coeffs))
+    sr = 0.0
+    si = 0.0
+    for n = length(eta_coeffs):-1:1
+        p = (dn^-re) * eta_coeffs[n]
+        lnn = -im * log(dn)
+        sr += p * cos(lnn)
+        si += p * sin(lnn)
+        dn -= 1
+    end
+    if reflect
+        z = complex(re, im)
+        b = 2.0 - 2.0^complex(re+1,im)
+        
+        f = 2.0^z - 2
+        piz = pi^z
+        
+        b = b/f/piz
+        
+        return complex(sr,si) * exp(lgamma(z)) * b * cos(pi/2*z)
+    end
+    return complex(sr, si)
+end
+
+eta(x::Real) = real(eta(complex(float64(x))))
+
+function zeta(z::Complex)
+    zz = 2.0^z
+    eta(z) * zz/(zz-2)
+end
+
+zeta(x::Real) = real(zeta(complex(float64(x))))