diff --git a/base/docs/helpdb/Base.jl b/base/docs/helpdb/Base.jl
index 487488f77dd59..83eeb1b5c9fc3 100644
--- a/base/docs/helpdb/Base.jl
+++ b/base/docs/helpdb/Base.jl
@@ -9579,14 +9579,6 @@ Riemann zeta function ``\\zeta(s)``.
 """
 zeta(s)
 
-"""
-    zeta(s, z)
-
-Hurwitz zeta function ``\\zeta(s, z)``.  (This is equivalent to the Riemann zeta function
-``\\zeta(s)`` for the case of `z=1`.)
-"""
-zeta(s,z)
-
 """
     A_mul_Bt(A, B)
 
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 67a8741bf6f2f..6cdf38412c569 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -220,15 +220,21 @@ macro pg_horner(x, m, p...)
     :(($me + 1) * ($(p[1]) + $xe * $ex))
 end
 
-# compute inv(oftype(x, y)) efficiently, choosing the correct branch cut
-inv_oftype(x::Complex, y::Complex) = oftype(x, inv(y))
-function inv_oftype(x::Complex, y::Real)
-    yi = inv(y) # using real arithmetic for efficiency
-    oftype(x, Complex(yi, -zero(yi))) # get correct sign of zero!
+# compute oftype(x, y)^p efficiently, choosing the correct branch cut
+pow_oftype(x, y, p) = oftype(x, y)^p
+pow_oftype(x::Complex, y::Real, p::Complex) = oftype(x, y^p)
+function pow_oftype(x::Complex, y::Real, p::Real)
+    if p >= 0
+        # note: this will never be called for y < 0,
+        # which would throw an error for non-integer p here
+        return oftype(x, y^p)
+    else
+        yp = y^-p # use real power for efficiency
+        return oftype(x, Complex(yp, -zero(yp))) # get correct sign of zero!
+    end
 end
-inv_oftype(x::Real, y::Real) = oftype(x, inv(y))
 
-# Hurwitz zeta function, which is related to polygamma
+# Generalized zeta function, which is related to polygamma
 # (at least for integer m > 0 and real(z) > 0) by:
 #    polygamma(m, z) = (-1)^(m+1) * gamma(m+1) * zeta(m+1, z).
 # Our algorithm for the polygamma is just the m-th derivative
@@ -238,15 +244,27 @@ inv_oftype(x::Real, y::Real) = oftype(x, inv(y))
 # So identifying the (something) with the -zeta function, we get
 # the zeta function for free and might as well export it, especially
 # since this is a common generalization of the Riemann zeta function
-# (which Julia already exports).
+# (which Julia already exports).   Note that this geneneralization
+# is equivalent to Mathematica's Zeta[s,z], and is equivalent to the
+# Hurwitz zeta function for real(z) > 0.
+
+"""
+    zeta(s, z)
+
+Generalized zeta function ``\\zeta(s, z)``, defined
+by the sum ``\\sum_{k=0}^\\infty ((k+z)^2)^{-s/2}``, where
+any term with ``k+z=0`` is excluded.  For ``\\Re z > 0``,
+this definition is equivalent to the Hurwitz zeta function
+``\\sum_{k=0}^\\infty (k+z)^{-s}``.   For ``z=1``, it yields
+the Riemann zeta function ``\\zeta(s)``.
+"""
+zeta(s,z)
+
 function zeta(s::Union{Int,Float64,Complex{Float64}},
               z::Union{Float64,Complex{Float64}})
     ζ = zero(promote_type(typeof(s), typeof(z)))
 
-    # like sqrt, require complex inputs to get complex outputs
-    !isa(s,Integer) && isa(ζ, Real) && z < 0 && throw(DomainError())
-
-    z == 1 && return oftype(ζ, zeta(s))
+    (z == 1 || z == 0) && return oftype(ζ, zeta(s))
     s == 2 && return oftype(ζ, trigamma(z))
 
     x = real(z)
@@ -263,9 +281,13 @@ function zeta(s::Union{Int,Float64,Complex{Float64}},
         end
         throw(DomainError()) # nothing clever to return
     end
-
-    # We need a different algorithm for the real(s) < 1 domain
-    real(s) < 1 && throw(ArgumentError("order $s < 1 is not implemented (issue #7228)"))
+    if isnan(x)
+        if imag(z)==0 && imag(s)==0
+            return oftype(ζ, x)
+        else
+            return oftype(ζ, Complex(x,x))
+        end
+    end
 
     m = s - 1
 
@@ -275,32 +297,55 @@ function zeta(s::Union{Int,Float64,Complex{Float64}},
     # Note: we multiply by -(-1)^m m! in polygamma below, so this factor is
     #       pulled out of all of our derivatives.
 
-    isnan(x) && return oftype(ζ, imag(z)==0 && isa(s,Int) ? x : Complex(x,x))
-
-    cutoff = 7 + real(m) + imag(m) # TODO: this cutoff is too conservative?
+    cutoff = 7 + real(m) + abs(imag(m)) # TODO: this cutoff is too conservative?
     if x < cutoff
         # shift using recurrence formula
         xf = floor(x)
-        if x <= 0 && xf == z
-            if isa(s, Int)
-                iseven(s) && return oftype(ζ, Inf)
-                x == 0 && return oftype(ζ, inv(x))
-            end
-            throw(DomainError()) # or return NaN?
-        end
         nx = Int(xf)
         n = ceil(Int,cutoff - nx)
-        ζ += inv_oftype(ζ, z)^s
-        for ν = -nx:-1:1
-            ζₒ= ζ
-            ζ += inv_oftype(ζ, z + ν)^s
-            ζ == ζₒ && break # prevent long loop for large -x > 0
-                             # FIXME: still slow for small m, large Im(z)
+        minus_s = -s
+        if nx < 0 # x < 0
+            # need to use (-z)^(-s) recurrence to be correct for real z < 0
+            # [the general form of the recurrence term is (z^2)^(-s/2)]
+            minus_z = -z
+            ζ += pow_oftype(ζ, minus_z, minus_s) # ν = 0 term
+            if xf != z
+                ζ += pow_oftype(ζ, z - nx, minus_s) # real(z - nx) > 0, so use correct branch cut
+                # otherwise, if xf==z, then the definition skips this term
+            end
+            # do loop in different order, depending on the sign of s,
+            # so that we are looping from largest to smallest summands and
+            # can halt the loop early if possible; see issue #15946
+            # FIXME: still slow for small m, large Im(z)
+            if real(s) > 0
+                for ν in -nx-1:-1:1
+                    ζₒ= ζ
+                    ζ += pow_oftype(ζ, minus_z - ν, minus_s)
+                    ζ == ζₒ && break # prevent long loop for large -x > 0
+                end
+            else
+                for ν in 1:-nx-1
+                    ζₒ= ζ
+                    ζ += pow_oftype(ζ, minus_z - ν, minus_s)
+                    ζ == ζₒ && break # prevent long loop for large -x > 0
+                end
+            end
+        else # x ≥ 0 && z != 0
+            ζ += pow_oftype(ζ, z, minus_s)
         end
-        for ν = max(1,1-nx):n-1
-            ζₒ= ζ
-            ζ += inv_oftype(ζ, z + ν)^s
-            ζ == ζₒ && break # prevent long loop for large m
+        # loop order depends on sign of s, as above
+        if real(s) > 0
+            for ν in max(1,1-nx):n-1
+                ζₒ= ζ
+                ζ += pow_oftype(ζ, z + ν, minus_s)
+                ζ == ζₒ && break # prevent long loop for large m
+            end
+        else
+            for ν in n-1:-1:max(1,1-nx)
+                ζₒ= ζ
+                ζ += pow_oftype(ζ, z + ν, minus_s)
+                ζ == ζₒ && break # prevent long loop for large m
+            end
         end
         z += n
     end
diff --git a/doc/stdlib/math.rst b/doc/stdlib/math.rst
index a6353a50d9508..4d3be3e2b70d7 100644
--- a/doc/stdlib/math.rst
+++ b/doc/stdlib/math.rst
@@ -1674,7 +1674,7 @@ Mathematical Functions
 
    .. Docstring generated from Julia source
 
-   Hurwitz zeta function :math:`\zeta(s, z)`\ .  (This is equivalent to the Riemann zeta function :math:`\zeta(s)` for the case of ``z=1``\ .)
+   Generalized zeta function :math:`\zeta(s, z)`\ , defined by the sum :math:`\sum_{k=0}^\infty ((k+z)^2)^{-s/2}`\ , where any term with :math:`k+z=0` is excluded.  For :math:`\Re z > 0`\ , this definition is equivalent to the Hurwitz zeta function :math:`\sum_{k=0}^\infty (k+z)^{-s}`\ .   For :math:`z=1`\ , it yields the Riemann zeta function :math:`\zeta(s)`\ .
 
 .. function:: ndigits(n, b)
 
diff --git a/test/math.jl b/test/math.jl
index 1c2c89e9d45c1..51bcd915edddb 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -615,52 +615,71 @@ end
 @test set_zero_subnormals(false)
 @test !get_zero_subnormals()
 
-# useful test functions for relative error
-err(z, x) = abs(z - x) / abs(x)
+# useful test functions for relative error, which differ from isapprox
+# in that errc separately looks at the real and imaginay parts
+err(z, x) = z == x ? 0.0 : abs(z - x) / abs(x)
 errc(z, x) = max(err(real(z),real(x)), err(imag(z),imag(x)))
+≅(a,b) = errc(a,b) ≤ 1e-13
 
 for x in -10.2:0.3456:50
     @test 1e-12 > err(digamma(x+0im), digamma(x))
 end
 
 # digamma, trigamma, polygamma & zeta test cases (compared to Wolfram Alpha)
-@test 1e-13 > err(digamma(7+0im), 1.872784335098467139393487909917597568957840664060076401194232)
-@test 1e-13 > errc(digamma(7im), 1.94761433458434866917623737015561385331974500663251349960124 + 1.642224898223468048051567761191050945700191089100087841536im)
-@test 1e-13 > errc(digamma(-3.2+0.1im), 4.65022505497781398615943030397508454861261537905047116427511+2.32676364843128349629415011622322040021960602904363963042380im)
-@test 1e-13 > err(trigamma(8+0im), 0.133137014694031425134546685920401606452509991909746283540546)
-@test 1e-13 > errc(trigamma(8im), -0.0078125000000000000029194973110119898029284994355721719150 - 0.12467345030312762782439017882063360876391046513966063947im)
-@test 1e-13 > errc(trigamma(-3.2+0.1im), 15.2073506449733631753218003030676132587307964766963426965699+15.7081038855113567966903832015076316497656334265029416039199im)
-@test 1e-13 > err(polygamma(2, 8.1+0im), -0.01723882695611191078960494454602091934457319791968308929600)
-@test 1e-13 > errc(polygamma(30, 8.1+2im), -2722.8895150799704384107961215752996280795801958784600407589+6935.8508929338093162407666304759101854270641674671634631058im)
-@test 1e-13 > errc(polygamma(3, 2.1+1im), 0.00083328137020421819513475400319288216246978855356531898998-0.27776110819632285785222411186352713789967528250214937861im)
+@test digamma(7+0im) ≅ 1.872784335098467139393487909917597568957840664060076401194232
+@test digamma(7im) ≅ 1.94761433458434866917623737015561385331974500663251349960124 + 1.642224898223468048051567761191050945700191089100087841536im
+@test digamma(-3.2+0.1im) ≅ 4.65022505497781398615943030397508454861261537905047116427511+2.32676364843128349629415011622322040021960602904363963042380im
+@test trigamma(8+0im) ≅ 0.133137014694031425134546685920401606452509991909746283540546
+@test trigamma(8im) ≅ -0.0078125000000000000029194973110119898029284994355721719150 - 0.12467345030312762782439017882063360876391046513966063947im
+@test trigamma(-3.2+0.1im) ≅ 15.2073506449733631753218003030676132587307964766963426965699+15.7081038855113567966903832015076316497656334265029416039199im
+@test polygamma(2, 8.1+0im) ≅ -0.01723882695611191078960494454602091934457319791968308929600
+@test polygamma(30, 8.1+2im) ≅ -2722.8895150799704384107961215752996280795801958784600407589+6935.8508929338093162407666304759101854270641674671634631058im
+@test polygamma(3, 2.1+1im) ≅ 0.00083328137020421819513475400319288216246978855356531898998-0.27776110819632285785222411186352713789967528250214937861im
 @test 1e-11 > err(polygamma(3, -4.2 + 2im),-0.0037752884324358856340054736472407163991189965406070325067-0.018937868838708874282432870292420046797798431078848805822im)
-@test 1e-13 > err(polygamma(13, 5.2 - 2im), 0.08087519202975913804697004241042171828113370070289754772448-0.2300264043021038366901951197725318713469156789541415899307im)
+@test polygamma(13, 5.2 - 2im) ≅ 0.08087519202975913804697004241042171828113370070289754772448-0.2300264043021038366901951197725318713469156789541415899307im
 @test 1e-11 > err(polygamma(123, -47.2 + 0im), 5.7111648667225422758966364116222590509254011308116701029e291)
-@test 1e-13 > errc(zeta(4.1+0.3im, -3.2+0.1im), -461.95403678374488506025596495576748255121001107881278765917+926.02552636148651929560277856510991293536052745360005500774im)
-@test 1e-13 > errc(zeta(4.1+0.3im, 3.2+0.1im), 0.0121197525131633219465301571139288562254218365173899270675-0.00687228692565614267981577154948499247518236888933925740902im)
-@test 1e-13 > errc(zeta(4.1, 3.2+0.1im),0.0137637451187986846516125754047084829556100290057521276517-0.00152194599531628234517456529686769063828217532350810111482im)
-@test 1e-12 > errc(zeta(1.0001, -4.5e2+3.2im), 9993.89099199843392251301993718413132850540848778561412270571-3.13257480938495907945892330398176989805350557816701044268548im)
-@test_throws DomainError zeta(3.1,-4.2)
-@test 1e-13 > errc(zeta(3.1,-4.2+0im), -138.06320182025311080661516120845508778572835942189570145952+45.586579397698817209431034568162819207622092308850063038062im)
+@test zeta(4.1+0.3im, -3.2+0.1im) ≅ -281.34474134962502296077659347175501181994490498591796647 + 286.55601240093672668066037366170168712249413003222992205im
+@test zeta(4.1+0.3im, 3.2+0.1im) ≅ 0.0121197525131633219465301571139288562254218365173899270675-0.00687228692565614267981577154948499247518236888933925740902im
+@test zeta(4.1, 3.2+0.1im) ≅ 0.0137637451187986846516125754047084829556100290057521276517-0.00152194599531628234517456529686769063828217532350810111482im
+@test 1e-12 > errc(zeta(1.0001, -4.5e2+3.2im), 10003.765660925877260544923069342257387254716966134385170 - 0.31956240712464746491659767831985629577542514145649468090im)
+@test zeta(3.1,-4.2) ≅ zeta(3.1,-4.2+0im) ≅ 149.7591329008219102939352965761913107060718455168339040295
 @test 1e-15 > errc(zeta(3.1+0im,-4.2), zeta(3.1,-4.2+0im))
-@test 1e-13 > errc(zeta(3.1,4.2), 0.029938344862645948405021260567725078588893266227472565010234)
-@test 1e-13 > err(zeta(27, 3.1), 5.413318813037879056337862215066960774064332961282599376e-14)
-@test 1e-13 > err(zeta(27, 2), 7.4507117898354294919810041706041194547190318825658299932e-9)
-@test 1e-12 > err(zeta(27, -105.3), -1.311372652244914148556295810515903234635727465138859603e14)
+@test zeta(3.1,4.2) ≅ 0.029938344862645948405021260567725078588893266227472565010234
+@test zeta(27, 3.1) ≅ 5.413318813037879056337862215066960774064332961282599376e-14
+@test zeta(27, 2) ≅ 7.4507117898354294919810041706041194547190318825658299932e-9
+@test 1e-12 > err(zeta(27, -105.3), 1.3113726525492708826840989036205762823329453315093955e14)
 @test polygamma(4, -3.1+Inf*im) == polygamma(4, 3.1+Inf*im) == 0
 @test polygamma(4, -0.0) == Inf == -polygamma(4, +0.0)
-@test zeta(4, +0.0) == Inf == zeta(4, -0.0)
-@test zeta(5, +0.0) == Inf == -zeta(5, -0.0)
+@test zeta(4, +0.0) == zeta(4, -0.0) ≅ pi^4 / 90
+@test zeta(5, +0.0) == zeta(5, -0.0) ≅ 1.036927755143369926331365486457034168057080919501912811974
 @test zeta(Inf, 1.) == 1
 @test zeta(Inf, 2.) == 0
 @test isnan(zeta(NaN, 1.))
 @test isa([digamma(x) for x in [1.0]], Vector{Float64})
 @test isa([trigamma(x) for x in [1.0]], Vector{Float64})
 @test isa([polygamma(3,x) for x in [1.0]], Vector{Float64})
-@test 1e-13 > errc(zeta(2 + 1im, -1.1), zeta(2 + 1im, -1.1+0im))
-@test 1e-13 > errc(zeta(2 + 1im, -1.1), -1525.8095173321060982383023516086563741006869909580583246557 + 1719.4753293650912305811325486980742946107143330321249869576im)
+@test zeta(2 + 1im, -1.1) ≅ zeta(2 + 1im, -1.1+0im) ≅ -64.580137707692178058665068045847533319237536295165484548 + 73.992688148809018073371913557697318846844796582012921247im
 @test_approx_eq polygamma(3,5) polygamma(3,5.)
 
+@test zeta(-3.0, 7.0) ≅ -52919/120
+@test zeta(-3.0, -7.0) ≅ 94081/120
+@test zeta(-3.1, 7.2) ≅ -587.457736596403704429103489502662574345388906240906317350719
+@test zeta(-3.1, -7.2) ≅ 1042.167459863862249173444363794330893294733001752715542569576
+@test zeta(-3.1, 7.0) ≅ -518.431785723446831868686653718848680989961581500352503093748
+@test zeta(-3.1, -7.0) ≅ 935.1284612957581823462429983411337864448020149908884596048161
+@test zeta(-3.1-0.1im, 7.2) ≅ -579.29752287650299181119859268736614017824853865655709516268 - 96.551907752211554484321948972741033127192063648337407683877im
+@test zeta(-3.1-0.1im, -7.2) ≅ 1025.17607931184231774568797674684390615895201417983173984531 + 185.732454778663400767583204948796029540252923367115805842138im
+@test zeta(-3.1-0.1im, 7.2 + 0.1im) ≅ -571.66133526455569807299410569274606007165253039948889085762 - 131.86744836357808604785415199791875369679879576524477540653im
+@test zeta(-3.1-0.1im, -7.2 + 0.1im) ≅ 1035.35760409421020754141207226034191979220047873089445768189 + 130.905870774271320475424492384335798304480814695778053731045im
+@test zeta(-3.1-0.1im, -7.0 + 0.1im) ≅ 929.546530292101383210555114424269079830017210969572819344670 + 113.646687807533854478778193456684618838875194573742062527301im
+@test zeta(-3.1, 7.2 + 0.1im) ≅ -586.61801005507638387063781112254388285799318636946559637115 - 36.148831292706044180986261734913443701649622026758378669700im
+@test zeta(-3.1, -7.2 + 0.1im) ≅ 1041.04241628770682295952302478199641560576378326778432301623 - 55.7154858634145071137760301929537184886497572057171143541058im
+@test zeta(-13.4, 4.1) ≅ -3.860040842156185186414774125656116135638705768861917e6
+@test zeta(3.2, -4) ≅ 2.317164896026427640718298719837102378116771112128525719078
+@test zeta(3.2, 0) ≅ 1.166773370984467020452550350896512026772734054324169010977
+@test zeta(-3.2+0.1im, 0.0) ≅ zeta(-3.2+0.1im, 0.0+0im) ≅ 0.0070547946138977701155565365569392198424378109226519905493 + 0.00076891821792430587745535285452496914239014050562476729610im
+@test zeta(-3.2, 0.0) ≅ zeta(-3.2, 0.0+0im) ≅ 0.007011972077091051091698102914884052994997144629191121056378
+
 @test @evalpoly(2,3,4,5,6) == 3+2*(4+2*(5+2*6)) == @evalpoly(2+0im,3,4,5,6)
 @test let evalcounts=0
           @evalpoly(begin