Skip to content

Instantly share code, notes, and snippets.

/-

Created February 9, 2017 19:47
Show Gist options
  • Save anonymous/84de6e5d71d989c549cb20671929302c to your computer and use it in GitHub Desktop.
Save anonymous/84de6e5d71d989c549cb20671929302c to your computer and use it in GitHub Desktop.
Download ZIP
Raw
-
From b0b58a80fdb3670d4ee9b39bfc83c63cadaefd9f Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Mon, 19 Sep 2016 23:55:30 +0800
Subject: [PATCH 01/14] =Work on #17474
---
base/special/gamma.jl | 107 ++++++++++++++++++++++++++++++++++----------------
test/math.jl | 9 +++++
2 files changed, 83 insertions(+), 33 deletions(-)
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 24e76ba..3bf5958 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -485,6 +485,8 @@ f16(z::Complex) = Complex32(z)
# If we really cared about single precision, we could make a faster
# Float32 version by truncating the Stirling series at a smaller cutoff.
+# and if we really cared about half precision, we could make a faster
+# Float16 version, by using a precomputed table look-up.
for (f,T) in ((:f32,Float32),(:f16,Float16))
@eval begin
zeta(s::Integer, z::Union{$T,Complex{$T}}) = $f(zeta(Int(s), f64(z)))
@@ -496,20 +498,71 @@ for (f,T) in ((:f32,Float32),(:f16,Float16))
end
end
-zeta(s::Integer, z::Number) = zeta(Int(s), f64(z))
-zeta(s::Number, z::Number) = zeta(f64(s), f64(z))
-for f in (:digamma, :trigamma)
+
+function zeta(s::Integer, z::Number)
+ x = float(z)
+ t = Int(s)
+ if typeof(x) === typeof(z) && typeof(t) === typeof(s)
+ throw MethodError(zeta,(s,t))
+ end
+ oftype(x, zeta(t, x))
+end
+
+function zeta(s::Number, z::Number)
+ x = float(z)
+ t = float(s)
+ if typeof(x) === typeof(z) && typeof(t) === typeof(s)
+ throw MethodError(zeta,(s,t))
+ end
+ oftype(x, zeta(t, x))
+end
+
+
+function polygamma(m::Integer, z::Number)
+ x = float(z)
+ typeof(x) == typeof(z) && throw(MethodError(polygamma, (m,z))
+ oftype(x,polygamma(m, x))
+end
+
+
+for f in (:digamma, :trigamma, :zeta, :eta, invdigamma)
@eval begin
- $f(z::Number) = $f(f64(z))
+ $f(z::Base.BitInteger) = $f(Float64(z))
+ $f(z::Float32) = Float32($f(Float64(z)))
+ $f(z::Float16) = Float16($f(Float64(z)))
+
+ function $f(z::Number)
+ x = float(z)
+ typeof(x) == typeof(z) && throw(MethodError($f, (z,)))
+ oftype(x, $f(x))
+ end
end
end
-polygamma(m::Integer, z::Number) = polygamma(m, f64(z))
-# Inverse digamma function:
-# Implementation of fixed point algorithm described in
-# "Estimating a Dirichlet distribution" by Thomas P. Minka, 2000
+for f in (:zeta, :eta)
+ @eval begin
+ $f{T<:Union{Base.BitInteger,Float32,Float16}}(x::t) = oftype(float(x), $f(Complex128(z)))
+
+ function $f(z::Complex)
+ x = float(z)
+ typeof(x) == typeof(z) && throw(MethodError($f, (z,)))
+ oftype(x, $f(Complex(x))
+ end
+ end
+end
+
+
+
+"""
+ invdigamma(x)
+
+Compute the inverse [`digamma`](:func:`digamma`) function of `x`.
+"""
function invdigamma(y::Float64)
- # Closed form initial estimates
+ # Implementation of fixed point algorithm described in
+ # "Estimating a Dirichlet distribution" by Thomas P. Minka, 2000
+
+ # Closed form initial estimates
if y >= -2.22
x_old = exp(y) + 0.5
x_new = x_old
@@ -530,14 +583,6 @@ function invdigamma(y::Float64)
return x_new
end
-invdigamma(x::Float32) = Float32(invdigamma(Float64(x)))
-
-"""
- invdigamma(x)
-
-Compute the inverse [`digamma`](:func:`digamma`) function of `x`.
-"""
-invdigamma(x::Real) = invdigamma(Float64(x))
"""
beta(x, y)
@@ -559,9 +604,16 @@ function ``\\log(|\\operatorname{B}(x,y)|)``.
"""
lbeta(x::Number, w::Number) = lgamma(x)+lgamma(w)-lgamma(x+w)
-# Riemann zeta function; algorithm is based on specializing the Hurwitz
-# zeta function above for z==1.
+
+"""
+ zeta(s)
+
+Riemann zeta function ``\\zeta(s)``.
+"""
function zeta(s::Union{Float64,Complex{Float64}})
+ # Riemann zeta function; algorithm is based on specializing the Hurwitz
+ # zeta function above for z==1.
+
# blows up to ±Inf, but get correct sign of imaginary zero
s == 1 && return NaN + zero(s) * imag(s)
@@ -606,16 +658,13 @@ function zeta(s::Union{Float64,Complex{Float64}})
return ζ
end
-zeta(x::Integer) = zeta(Float64(x))
-zeta(x::Real) = oftype(float(x),zeta(Float64(x)))
-"""
- zeta(s)
-Riemann zeta function ``\\zeta(s)``.
"""
-zeta(z::Complex) = oftype(float(z),zeta(Complex128(z)))
+ eta(x)
+Dirichlet eta function ``\\eta(s) = \\sum^\\infty_{n=1}(-1)^{n-1}/n^{s}``.
+"""
function eta(z::Union{Float64,Complex{Float64}})
δz = 1 - z
if abs(real(δz)) + abs(imag(δz)) < 7e-3 # Taylor expand around z==1
@@ -630,12 +679,4 @@ function eta(z::Union{Float64,Complex{Float64}})
return -zeta(z) * expm1(0.6931471805599453094172321214581765*δz)
end
end
-eta(x::Integer) = eta(Float64(x))
-eta(x::Real) = oftype(float(x),eta(Float64(x)))
-
-"""
- eta(x)
-Dirichlet eta function ``\\eta(s) = \\sum^\\infty_{n=1}(-1)^{n-1}/n^{s}``.
-"""
-eta(z::Complex) = oftype(float(z),eta(Complex128(z)))
diff --git a/test/math.jl b/test/math.jl
index 132dd53..66ac7f8 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -941,3 +941,12 @@ end
# no domain error is thrown for negative values
@test invoke(cbrt, Tuple{AbstractFloat}, -1.0) == -1.0
+
+# issue #17474
+@test typeof(eta(big"2")) == BigFloat
+@test typeof(zeta(big"2")) == BigFloat
+@test typeof(digamma(big"2")) == BigFloat
+@test_throws MethodError trigamma(big"2")
+@test_throws MethodError trigamma(big"2.0")
+@test_throws MethodError invdigamma(big"2")
+@test_throws MethodError invdiamma(big"2.0")
From ba62160a6d83941e7f0fa9cb1e1cb5bd0ba0a809 Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Tue, 4 Oct 2016 11:10:09 +0800
Subject: [PATCH 02/14] =Work on #17474 squashme
---
base/special/gamma.jl | 181 ++++++++++++++++++++++++++------------------------
base/sysimg.jl | 24 ++-----
test/math.jl | 6 +-
3 files changed, 105 insertions(+), 106 deletions(-)
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 3bf5958..1552c33 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -1,4 +1,5 @@
# This file is a part of Julia. License is MIT: http://julialang.org/license
+typealias ComplexOrReal{T} Union{T,Complex{T}}
gamma(x::Float64) = nan_dom_err(ccall((:tgamma,libm), Float64, (Float64,), x), x)
gamma(x::Float32) = nan_dom_err(ccall((:tgammaf,libm), Float32, (Float32,), x), x)
@@ -147,7 +148,7 @@ gamma(z::Complex) = exp(lgamma(z))
Compute the digamma function of `x` (the logarithmic derivative of `gamma(x)`)
"""
-function digamma(z::Union{Float64,Complex{Float64}})
+function digamma(z::ComplexOrReal{Float64})
# Based on eq. (12), without looking at the accompanying source
# code, of: K. S. Kölbig, "Programs for computing the logarithm of
# the gamma function, and the digamma function, for complex
@@ -181,7 +182,7 @@ end
Compute the trigamma function of `x` (the logarithmic second derivative of `gamma(x)`).
"""
-function trigamma(z::Union{Float64,Complex{Float64}})
+function trigamma(z::ComplexOrReal{Float64})
# via the derivative of the Kölbig digamma formulation
x = real(z)
if x <= 0 # reflection formula
@@ -343,8 +344,8 @@ the Riemann zeta function ``\\zeta(s)``.
"""
zeta(s,z)
-function zeta(s::Union{Int,Float64,Complex{Float64}},
- z::Union{Float64,Complex{Float64}})
+function zeta(s::Union{Int, ComplexOrReal{Float64}},
+ z::ComplexOrReal{Float64})
ζ = zero(promote_type(typeof(s), typeof(z)))
(z == 1 || z == 0) && return oftype(ζ, zeta(s))
@@ -449,7 +450,7 @@ end
Compute the polygamma function of order `m` of argument `x` (the `(m+1)th` derivative of the
logarithm of `gamma(x)`)
"""
-function polygamma(m::Integer, z::Union{Float64,Complex{Float64}})
+function polygamma(m::Integer, z::ComplexOrReal{Float64})
m == 0 && return digamma(z)
m == 1 && return trigamma(z)
@@ -475,83 +476,6 @@ function polygamma(m::Integer, z::Union{Float64,Complex{Float64}})
end
end
-# TODO: better way to do this
-f64(x::Real) = Float64(x)
-f64(z::Complex) = Complex128(z)
-f32(x::Real) = Float32(x)
-f32(z::Complex) = Complex64(z)
-f16(x::Real) = Float16(x)
-f16(z::Complex) = Complex32(z)
-
-# If we really cared about single precision, we could make a faster
-# Float32 version by truncating the Stirling series at a smaller cutoff.
-# and if we really cared about half precision, we could make a faster
-# Float16 version, by using a precomputed table look-up.
-for (f,T) in ((:f32,Float32),(:f16,Float16))
- @eval begin
- zeta(s::Integer, z::Union{$T,Complex{$T}}) = $f(zeta(Int(s), f64(z)))
- zeta(s::Union{Float64,Complex128}, z::Union{$T,Complex{$T}}) = zeta(s, f64(z))
- zeta(s::Number, z::Union{$T,Complex{$T}}) = $f(zeta(f64(s), f64(z)))
- polygamma(m::Integer, z::Union{$T,Complex{$T}}) = $f(polygamma(Int(m), f64(z)))
- digamma(z::Union{$T,Complex{$T}}) = $f(digamma(f64(z)))
- trigamma(z::Union{$T,Complex{$T}}) = $f(trigamma(f64(z)))
- end
-end
-
-
-function zeta(s::Integer, z::Number)
- x = float(z)
- t = Int(s)
- if typeof(x) === typeof(z) && typeof(t) === typeof(s)
- throw MethodError(zeta,(s,t))
- end
- oftype(x, zeta(t, x))
-end
-
-function zeta(s::Number, z::Number)
- x = float(z)
- t = float(s)
- if typeof(x) === typeof(z) && typeof(t) === typeof(s)
- throw MethodError(zeta,(s,t))
- end
- oftype(x, zeta(t, x))
-end
-
-
-function polygamma(m::Integer, z::Number)
- x = float(z)
- typeof(x) == typeof(z) && throw(MethodError(polygamma, (m,z))
- oftype(x,polygamma(m, x))
-end
-
-
-for f in (:digamma, :trigamma, :zeta, :eta, invdigamma)
- @eval begin
- $f(z::Base.BitInteger) = $f(Float64(z))
- $f(z::Float32) = Float32($f(Float64(z)))
- $f(z::Float16) = Float16($f(Float64(z)))
-
- function $f(z::Number)
- x = float(z)
- typeof(x) == typeof(z) && throw(MethodError($f, (z,)))
- oftype(x, $f(x))
- end
- end
-end
-
-for f in (:zeta, :eta)
- @eval begin
- $f{T<:Union{Base.BitInteger,Float32,Float16}}(x::t) = oftype(float(x), $f(Complex128(z)))
-
- function $f(z::Complex)
- x = float(z)
- typeof(x) == typeof(z) && throw(MethodError($f, (z,)))
- oftype(x, $f(Complex(x))
- end
- end
-end
-
-
"""
invdigamma(x)
@@ -559,10 +483,10 @@ end
Compute the inverse [`digamma`](:func:`digamma`) function of `x`.
"""
function invdigamma(y::Float64)
- # Implementation of fixed point algorithm described in
- # "Estimating a Dirichlet distribution" by Thomas P. Minka, 2000
+ # Implementation of fixed point algorithm described in
+ # "Estimating a Dirichlet distribution" by Thomas P. Minka, 2000
- # Closed form initial estimates
+ # Closed form initial estimates
if y >= -2.22
x_old = exp(y) + 0.5
x_new = x_old
@@ -610,7 +534,7 @@ lbeta(x::Number, w::Number) = lgamma(x)+lgamma(w)-lgamma(x+w)
Riemann zeta function ``\\zeta(s)``.
"""
-function zeta(s::Union{Float64,Complex{Float64}})
+function zeta(s::ComplexOrReal{Float64})
# Riemann zeta function; algorithm is based on specializing the Hurwitz
# zeta function above for z==1.
@@ -665,7 +589,7 @@ end
Dirichlet eta function ``\\eta(s) = \\sum^\\infty_{n=1}(-1)^{n-1}/n^{s}``.
"""
-function eta(z::Union{Float64,Complex{Float64}})
+function eta(z::ComplexOrReal{Float64})
δz = 1 - z
if abs(real(δz)) + abs(imag(δz)) < 7e-3 # Taylor expand around z==1
return 0.6931471805599453094172321214581765 *
@@ -680,3 +604,86 @@ function eta(z::Union{Float64,Complex{Float64}})
end
end
+
+# Converting types that we can convert, and not ones we can not
+# If we really cared about single precision, we could make a faster
+# Float32 version by truncating the Stirling series at a smaller cutoff.
+# and if we really cared about half precision, we could make a faster
+# Float16 version, by using a precomputed table look-up.
+
+f64(x::Union{Base.BitInteger, Float16, Float32, Float64}) = Float64(x)
+f64(x::Complex{Union{Base.BitInteger, Float16, Float32, Float64}}) = Complex128(x)
+ofpromotedtype(as, c) = oftype(promote(as...), c)
+
+for T in (Float32, Float16)
+ @eval begin
+ polygamma(m::Integer, z::ComplexOrReal{T}) = ofpromotedtype((m,z), polygamma(Int(m), f64(z)))
+ digamma(z::ComplexOrReal{$T}) = oftype(z, digamma(f64(z)))
+ trigamma(z::ComplexOrReal{$T}) = oftype(z, trigamma(f64(z)))
+ zeta(s::Integer, z::ComplexOrReal{$T}}) = ofpromotedtype((s,z), zeta(Int(s), f64(z)))
+ end
+end
+
+function zeta(s::ComplexOrReal{Union{Float16, Float32}},
+ z::ComplexOrReal{Union{Float16, Float32, Float64, Base.BitInteger})
+ ofpromotedtype((s, z), zeta(f64(s), f64(z)))
+end
+
+
+function zeta(s::Integer, z::Number)
+ x = float(z)
+ t = Int(s)
+ if typeof(x) === typeof(z) && typeof(t) === typeof(s)
+ # There is nothing to fallback to, since this didn't work
+ throw(MethodError(zeta,(s,t)))
+ end
+ ofpromotedtype((x,y), zeta(t, x))
+end
+
+function zeta(s::Number, z::Number)
+ x = float(z)
+ t = float(s)
+ if typeof(x) === typeof(z) && typeof(t) === typeof(s)
+ # There is nothing to fallback to, since this didn't work
+ throw(MethodError(zeta,(s,t)))
+ end
+ ofpromotedtype((x,t), zeta(t, x))
+end
+
+
+function polygamma(m::Integer, z::Number)
+ x = float(z)
+ typeof(x) == typeof(z) && throw(MethodError(polygamma, (m,z)))
+ # There is nothing to fallback to, since this didn't work
+ oftype(x, polygamma(m, x))
+end
+
+
+for f in (:digamma, :trigamma, :zeta, :eta, :invdigamma)
+ @eval begin
+ $f(z::Base.BitInteger) = $f(Float64(z))
+ $f(z::Float32) = Float32($f(Float64(z)))
+ $f(z::Float16) = Float16($f(Float64(z)))
+
+ function $f(z::Number)
+ x = float(z)
+ typeof(x) == typeof(z) && throw(MethodError($f, (z,)))
+ # There is nothing to fallback to, since this didn't work
+ oftype(x, $f(x))
+ end
+ end
+end
+
+for f in (:zeta, :eta)
+ @eval begin
+ $f{T<:Union{Base.BitInteger,Float32,Float16}}(z::Complex{T}) = oftype(float(z), $f(Complex128(z)))
+
+ function $f(z::Complex)
+ x = float(z)
+ typeof(x) == typeof(z) && throw(MethodError($f, (z,)))
+ # There is nothing to fallback to, since this didn't work
+ oftype(x, $f(x))
+ end
+ end
+end
+
diff --git a/base/sysimg.jl b/base/sysimg.jl
index 8e5dd92..6821e37 100644
--- a/base/sysimg.jl
+++ b/base/sysimg.jl
@@ -4,19 +4,7 @@ baremodule Base
using Core.Intrinsics
ccall(:jl_set_istopmod, Void, (Bool,), true)
-function include(path::AbstractString)
- local result
- if INCLUDE_STATE === 1
- result = Core.include(path)
- elseif INCLUDE_STATE === 2
- result = _include(path)
- elseif INCLUDE_STATE === 3
- result = include_from_node1(path)
- end
- result
-end
-INCLUDE_STATE = 1 # include = Core.include
-
+include = Core.include
include("coreio.jl")
eval(x) = Core.eval(Base,x)
@@ -80,9 +68,9 @@ end
Symbol(x...) = Symbol(string(x...))
# array structures
-include("array.jl")
include("abstractarray.jl")
include("subarray.jl")
+include("array.jl")
# Array convenience converting constructors
(::Type{Array{T}}){T}(m::Integer) = Array{T,1}(Int(m))
@@ -114,6 +102,7 @@ include("multinverses.jl")
using .MultiplicativeInverses
include("abstractarraymath.jl")
include("arraymath.jl")
+include("float16.jl")
# SIMD loops
include("simdloop.jl")
@@ -216,7 +205,7 @@ include("permuteddimsarray.jl")
using .PermutedDimsArrays
let SOURCE_PATH = ""
- global function _include(path)
+ global include = function(path)
prev = SOURCE_PATH
path = joinpath(dirname(prev),path)
SOURCE_PATH = path
@@ -224,7 +213,6 @@ let SOURCE_PATH = ""
SOURCE_PATH = prev
end
end
-INCLUDE_STATE = 2 # include = _include (from lines above)
# reduction along dims
include("reducedim.jl") # macros in this file relies on string.jl
@@ -387,8 +375,8 @@ function __init__()
init_threadcall()
end
-INCLUDE_STATE = 3 # include = include_from_node1
-include("precompile.jl")
+include = include_from_node1
+#include("precompile.jl") #Don't commit me. Speed up testing l
end # baremodule Base
diff --git a/test/math.jl b/test/math.jl
index 66ac7f8..e1f85f8 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -949,4 +949,8 @@ end
@test_throws MethodError trigamma(big"2")
@test_throws MethodError trigamma(big"2.0")
@test_throws MethodError invdigamma(big"2")
-@test_throws MethodError invdiamma(big"2.0")
+@test_throws MethodError invdigamma(big"2.0")
+@test_throws MethodError eta(Complex(big"2"))
+@test_throws MethodError eta(Complex(big"2.0"))
+@test_throws MethodError zeta(Complex(big"2"))
+@test_throws MethodError zeta(Complex(big"2.0"))
From bf783c97ab87f8e820efbdff6fae2158255f2bde Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Wed, 5 Oct 2016 18:56:17 +0800
Subject: [PATCH 03/14] =Tests Passing
---
base/special/gamma.jl | 100 +++++++++++++++++++++++++++-----------------------
test/math.jl | 6 ++-
2 files changed, 59 insertions(+), 47 deletions(-)
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 1552c33..7605c5c 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -611,35 +611,75 @@ end
# and if we really cared about half precision, we could make a faster
# Float16 version, by using a precomputed table look-up.
-f64(x::Union{Base.BitInteger, Float16, Float32, Float64}) = Float64(x)
-f64(x::Complex{Union{Base.BitInteger, Float16, Float32, Float64}}) = Complex128(x)
-ofpromotedtype(as, c) = oftype(promote(as...), c)
+# Float16, and Float32 and their Complex equivalents can be cast to Float64
+# and results cast back. Similar for BitIntegers (eg Int32), but no casting back
+# Otherwise, we need to make things use their own `float` converting methods
+# and in those cases, we do not cast back either.
-for T in (Float32, Float16)
+
+ofpromotedtype(as::Tuple, c) = convert(promote_type(typeof.(as)...), c)
+ComplexOrRealUnion(TS...) = Union{(ComplexOrReal{T} for T in TS)...}
+
+const types_le_Float64 = (Float16, Float32, Float64, Base.BitInteger.types...)
+
+for T in types_le_Float64
+ @eval f64(x::Complex{$T}) = Complex128(x)
+ @eval f64(x::$T) = Float64(x)
+end
+
+
+for f in (:digamma, :trigamma, :zeta, :eta, :invdigamma)
@eval begin
- polygamma(m::Integer, z::ComplexOrReal{T}) = ofpromotedtype((m,z), polygamma(Int(m), f64(z)))
- digamma(z::ComplexOrReal{$T}) = oftype(z, digamma(f64(z)))
- trigamma(z::ComplexOrReal{$T}) = oftype(z, trigamma(f64(z)))
- zeta(s::Integer, z::ComplexOrReal{$T}}) = ofpromotedtype((s,z), zeta(Int(s), f64(z)))
+ $f(z::ComplexOrRealUnion(Base.BitInteger.types...)) = $f(f64(z))
+ $f(z::ComplexOrRealUnion(Float16,Float32)) = oftype(z, $f(f64(z)))
+
+ function $f(z::Number)
+ x = float(z)
+ typeof(x) == typeof(z) && throw(MethodError($f, (z,)))
+ # There is nothing to fallback to, since this didn't work
+ $f(x)
+ end
end
end
-function zeta(s::ComplexOrReal{Union{Float16, Float32}},
- z::ComplexOrReal{Union{Float16, Float32, Float64, Base.BitInteger})
- ofpromotedtype((s, z), zeta(f64(s), f64(z)))
+
+polygamma(m::Integer, z::ComplexOrRealUnion(Float16,Float32)) = oftype(z, polygamma(m, f64(z)))
+
+
+for T1 in types_le_Float64, T2 in types_le_Float64
+ if (T1 == T2 == Float64) || (T1 == Int && T2 == Float64)
+ continue # Avoid redefining base definition
+ # However this skips `zeta(::Complex{Int}, ::ComplexOrReal{Float64})`
+ # so that will need to be added back after
+ end
+
+ if T1<:Integer && T2<:Integer
+ @eval function zeta(s::ComplexOrReal{$T1}, z::ComplexOrReal{$T2})
+ zeta(f64(s), f64(z)) #Do not promote down to Integers
+ end
+ else
+ @eval function zeta(s::ComplexOrReal{$T1}, z::ComplexOrReal{$T2})
+ ofpromotedtype((s, z), zeta(f64(s), f64(z)))
+ end
+ end
end
+# this is the one definition that is skipped
+function zeta(s::Complex{Int}, z::ComplexOrReal{Float64})::Complex{Float64}
+ zeta(f64(s), f64(z))
+end
function zeta(s::Integer, z::Number)
x = float(z)
- t = Int(s)
+ t = Int(s) # One could worry here about converting a BigInteger into a Int32/Int64
if typeof(x) === typeof(z) && typeof(t) === typeof(s)
# There is nothing to fallback to, since this didn't work
throw(MethodError(zeta,(s,t)))
end
- ofpromotedtype((x,y), zeta(t, x))
+ zeta(t, x)
end
+
function zeta(s::Number, z::Number)
x = float(z)
t = float(s)
@@ -647,7 +687,7 @@ function zeta(s::Number, z::Number)
# There is nothing to fallback to, since this didn't work
throw(MethodError(zeta,(s,t)))
end
- ofpromotedtype((x,t), zeta(t, x))
+ zeta(t, x)
end
@@ -655,35 +695,5 @@ function polygamma(m::Integer, z::Number)
x = float(z)
typeof(x) == typeof(z) && throw(MethodError(polygamma, (m,z)))
# There is nothing to fallback to, since this didn't work
- oftype(x, polygamma(m, x))
-end
-
-
-for f in (:digamma, :trigamma, :zeta, :eta, :invdigamma)
- @eval begin
- $f(z::Base.BitInteger) = $f(Float64(z))
- $f(z::Float32) = Float32($f(Float64(z)))
- $f(z::Float16) = Float16($f(Float64(z)))
-
- function $f(z::Number)
- x = float(z)
- typeof(x) == typeof(z) && throw(MethodError($f, (z,)))
- # There is nothing to fallback to, since this didn't work
- oftype(x, $f(x))
- end
- end
-end
-
-for f in (:zeta, :eta)
- @eval begin
- $f{T<:Union{Base.BitInteger,Float32,Float16}}(z::Complex{T}) = oftype(float(z), $f(Complex128(z)))
-
- function $f(z::Complex)
- x = float(z)
- typeof(x) == typeof(z) && throw(MethodError($f, (z,)))
- # There is nothing to fallback to, since this didn't work
- oftype(x, $f(x))
- end
- end
+ polygamma(m, x)
end
-
diff --git a/test/math.jl b/test/math.jl
index e1f85f8..b22fb5f 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -936,8 +936,10 @@ end
end
end
-@test Base.Math.f32(complex(1.0,1.0)) == complex(Float32(1.),Float32(1.))
-@test Base.Math.f16(complex(1.0,1.0)) == complex(Float16(1.),Float16(1.))
+
+@test Base.Math.f64(complex(1f0,1f0)) == complex(1.0, 1.0)
+@test Base.Math.f64(1f0) == 1.0
+@test Base.Math.ofpromotedtype((complex(1f0, 1f0), 1.0), 5.0) == complex(5.0)
# no domain error is thrown for negative values
@test invoke(cbrt, Tuple{AbstractFloat}, -1.0) == -1.0
From 63fc253e70d900a9fa3a43605e0bc00ac5616b8f Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Wed, 5 Oct 2016 18:58:28 +0800
Subject: [PATCH 04/14] =revert changes to sysimg (rebase these out after)
---
base/sysimg.jl | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/base/sysimg.jl b/base/sysimg.jl
index 6821e37..01ac383 100644
--- a/base/sysimg.jl
+++ b/base/sysimg.jl
@@ -376,7 +376,7 @@ function __init__()
end
include = include_from_node1
-#include("precompile.jl") #Don't commit me. Speed up testing l
+include("precompile.jl")
end # baremodule Base
From f62b0491109b3ef5fb319c0f118dc9c7fe830f9e Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Wed, 5 Oct 2016 19:23:46 +0800
Subject: [PATCH 05/14] =Fix whitespace and line wrap
---
base/special/gamma.jl | 39 ++++++++++++++++++++++-----------------
1 file changed, 22 insertions(+), 17 deletions(-)
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 7605c5c..2a3e4bb 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -73,7 +73,8 @@ function lgamma(z::Complex{Float64})
else
return Complex(NaN, NaN)
end
- elseif x > 7 || yabs > 7 # use the Stirling asymptotic series for sufficiently large x or |y|
+ elseif x > 7 || yabs > 7
+ # use the Stirling asymptotic series for sufficiently large x or |y|
return lgamma_asymptotic(z)
elseif x < 0.1 # use reflection formula to transform to x > 0
if x == 0 && y == 0 # return Inf with the correct imaginary part for z == 0
@@ -107,7 +108,8 @@ function lgamma(z::Complex{Float64})
-2.2315475845357937976132853e-04,9.9457512781808533714662972e-05,
-4.4926236738133141700224489e-05,2.0507212775670691553131246e-05)
end
- # use recurrence relation lgamma(z) = lgamma(z+1) - log(z) to shift to x > 7 for asymptotic series
+ # use recurrence relation lgamma(z) = lgamma(z+1) - log(z )
+ # to shift to x > 7 for asymptotic series
shiftprod = Complex(x,yabs)
x += 1
sb = false # == signbit(imag(shiftprod)) == signbit(yabs)
@@ -394,7 +396,8 @@ function zeta(s::Union{Int, ComplexOrReal{Float64}},
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
+ ζ += 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,
@@ -447,8 +450,8 @@ end
"""
polygamma(m, x)
-Compute the polygamma function of order `m` of argument `x` (the `(m+1)th` derivative of the
-logarithm of `gamma(x)`)
+Compute the polygamma function of order `m` of argument `x`
+(the `(m+1)th` derivative of the logarithm of `gamma(x)`)
"""
function polygamma(m::Integer, z::ComplexOrReal{Float64})
m == 0 && return digamma(z)
@@ -511,7 +514,8 @@ end
"""
beta(x, y)
-Euler integral of the first kind ``\\operatorname{B}(x,y) = \\Gamma(x)\\Gamma(y)/\\Gamma(x+y)``.
+Euler integral of the first kind
+``\\operatorname{B}(x,y) = \\Gamma(x)\\Gamma(y)/\\Gamma(x+y)``.
"""
function beta(x::Number, w::Number)
yx, sx = lgamma_r(x)
@@ -611,10 +615,11 @@ end
# and if we really cared about half precision, we could make a faster
# Float16 version, by using a precomputed table look-up.
-# Float16, and Float32 and their Complex equivalents can be cast to Float64
-# and results cast back. Similar for BitIntegers (eg Int32), but no casting back
+# Float16, and Float32 and their Complex equivalents can be converted to Float64
+# and results converted back. Similar for BitIntegers (eg Int32), but no converting back
# Otherwise, we need to make things use their own `float` converting methods
-# and in those cases, we do not cast back either.
+# and in those cases, we do not convert back either as we assume
+# they also implement there own versions of the functions
ofpromotedtype(as::Tuple, c) = convert(promote_type(typeof.(as)...), c)
@@ -632,19 +637,19 @@ for f in (:digamma, :trigamma, :zeta, :eta, :invdigamma)
@eval begin
$f(z::ComplexOrRealUnion(Base.BitInteger.types...)) = $f(f64(z))
$f(z::ComplexOrRealUnion(Float16,Float32)) = oftype(z, $f(f64(z)))
-
+
function $f(z::Number)
x = float(z)
typeof(x) == typeof(z) && throw(MethodError($f, (z,)))
- # There is nothing to fallback to, since this didn't work
+ # There is nothing to fallback to, since this didn't change the argument types
$f(x)
end
end
end
-
-polygamma(m::Integer, z::ComplexOrRealUnion(Float16,Float32)) = oftype(z, polygamma(m, f64(z)))
-
+function polygamma(m::Integer, z::ComplexOrRealUnion(Float16,Float32))
+ oftype(z, polygamma(m, f64(z)))
+end
for T1 in types_le_Float64, T2 in types_le_Float64
if (T1 == T2 == Float64) || (T1 == Int && T2 == Float64)
@@ -664,12 +669,12 @@ for T1 in types_le_Float64, T2 in types_le_Float64
end
end
-# this is the one definition that is skipped
+# this is the one definition that is skipped
function zeta(s::Complex{Int}, z::ComplexOrReal{Float64})::Complex{Float64}
zeta(f64(s), f64(z))
end
-function zeta(s::Integer, z::Number)
+function zeta(s::Integer, z::Number)
x = float(z)
t = Int(s) # One could worry here about converting a BigInteger into a Int32/Int64
if typeof(x) === typeof(z) && typeof(t) === typeof(s)
@@ -680,7 +685,7 @@ function zeta(s::Integer, z::Number)
end
-function zeta(s::Number, z::Number)
+function zeta(s::Number, z::Number)
x = float(z)
t = float(s)
if typeof(x) === typeof(z) && typeof(t) === typeof(s)
From dca825686fa9c2e0ba6ad240a98f6b0a16e2c89b Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Wed, 5 Oct 2016 20:11:43 +0800
Subject: [PATCH 06/14] =haul in master (something went wrong in rebase?)
---
base/sysimg.jl | 22 +++++++++++++++++-----
1 file changed, 17 insertions(+), 5 deletions(-)
diff --git a/base/sysimg.jl b/base/sysimg.jl
index 01ac383..8e5dd92 100644
--- a/base/sysimg.jl
+++ b/base/sysimg.jl
@@ -4,7 +4,19 @@ baremodule Base
using Core.Intrinsics
ccall(:jl_set_istopmod, Void, (Bool,), true)
-include = Core.include
+function include(path::AbstractString)
+ local result
+ if INCLUDE_STATE === 1
+ result = Core.include(path)
+ elseif INCLUDE_STATE === 2
+ result = _include(path)
+ elseif INCLUDE_STATE === 3
+ result = include_from_node1(path)
+ end
+ result
+end
+INCLUDE_STATE = 1 # include = Core.include
+
include("coreio.jl")
eval(x) = Core.eval(Base,x)
@@ -68,9 +80,9 @@ end
Symbol(x...) = Symbol(string(x...))
# array structures
+include("array.jl")
include("abstractarray.jl")
include("subarray.jl")
-include("array.jl")
# Array convenience converting constructors
(::Type{Array{T}}){T}(m::Integer) = Array{T,1}(Int(m))
@@ -102,7 +114,6 @@ include("multinverses.jl")
using .MultiplicativeInverses
include("abstractarraymath.jl")
include("arraymath.jl")
-include("float16.jl")
# SIMD loops
include("simdloop.jl")
@@ -205,7 +216,7 @@ include("permuteddimsarray.jl")
using .PermutedDimsArrays
let SOURCE_PATH = ""
- global include = function(path)
+ global function _include(path)
prev = SOURCE_PATH
path = joinpath(dirname(prev),path)
SOURCE_PATH = path
@@ -213,6 +224,7 @@ let SOURCE_PATH = ""
SOURCE_PATH = prev
end
end
+INCLUDE_STATE = 2 # include = _include (from lines above)
# reduction along dims
include("reducedim.jl") # macros in this file relies on string.jl
@@ -375,7 +387,7 @@ function __init__()
init_threadcall()
end
-include = include_from_node1
+INCLUDE_STATE = 3 # include = include_from_node1
include("precompile.jl")
end # baremodule Base
From c96f57888416b4b8fdf7e0a9378576bc803f8038 Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Wed, 5 Oct 2016 20:57:38 +0800
Subject: [PATCH 07/14] =correct MethodErrors for #17474
---
base/special/gamma.jl | 16 ++++++++--------
test/math.jl | 11 +++++++++++
2 files changed, 19 insertions(+), 8 deletions(-)
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 2a3e4bb..a322d92 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -640,7 +640,7 @@ for f in (:digamma, :trigamma, :zeta, :eta, :invdigamma)
function $f(z::Number)
x = float(z)
- typeof(x) == typeof(z) && throw(MethodError($f, (z,)))
+ typeof(x) === typeof(z) && throw(MethodError($f, (z,)))
# There is nothing to fallback to, since this didn't change the argument types
$f(x)
end
@@ -675,22 +675,22 @@ function zeta(s::Complex{Int}, z::ComplexOrReal{Float64})::Complex{Float64}
end
function zeta(s::Integer, z::Number)
- x = float(z)
t = Int(s) # One could worry here about converting a BigInteger into a Int32/Int64
- if typeof(x) === typeof(z) && typeof(t) === typeof(s)
+ x = float(z)
+ if typeof(t) === typeof(s) && typeof(x) === typeof(z)
# There is nothing to fallback to, since this didn't work
- throw(MethodError(zeta,(s,t)))
+ throw(MethodError(zeta,(s,z)))
end
zeta(t, x)
end
function zeta(s::Number, z::Number)
- x = float(z)
t = float(s)
- if typeof(x) === typeof(z) && typeof(t) === typeof(s)
+ x = float(z)
+ if typeof(t) === typeof(s) && typeof(x) === typeof(z)
# There is nothing to fallback to, since this didn't work
- throw(MethodError(zeta,(s,t)))
+ throw(MethodError(zeta,(s,z)))
end
zeta(t, x)
end
@@ -698,7 +698,7 @@ end
function polygamma(m::Integer, z::Number)
x = float(z)
- typeof(x) == typeof(z) && throw(MethodError(polygamma, (m,z)))
+ typeof(x) === typeof(z) && throw(MethodError(polygamma, (m,z)))
# There is nothing to fallback to, since this didn't work
polygamma(m, x)
end
diff --git a/test/math.jl b/test/math.jl
index b22fb5f..bbc6b8f 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -948,11 +948,22 @@ end
@test typeof(eta(big"2")) == BigFloat
@test typeof(zeta(big"2")) == BigFloat
@test typeof(digamma(big"2")) == BigFloat
+
@test_throws MethodError trigamma(big"2")
@test_throws MethodError trigamma(big"2.0")
@test_throws MethodError invdigamma(big"2")
@test_throws MethodError invdigamma(big"2.0")
+
@test_throws MethodError eta(Complex(big"2"))
@test_throws MethodError eta(Complex(big"2.0"))
@test_throws MethodError zeta(Complex(big"2"))
@test_throws MethodError zeta(Complex(big"2.0"))
+@test_throws MethodError zeta(1.0,big"2")
+@test_throws MethodError zeta(1.0,big"2.0")
+@test_throws MethodError zeta(big"1.0",2.0)
+
+@test typeof(zeta(complex(1),2.0)) == Complex{Float64}
+@test typeof(polygamma(3, 0x2)) == Float64
+@test typeof(polygamma(big"3", 2f0)) == Float32
+@test typeof(zeta(big"1",2.0)) == Float64 #Is this really desirable behavour?
+
From 706773b022133ba7db1a6b8d99fd8a99c1ec25a5 Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Sat, 8 Oct 2016 22:44:13 +0800
Subject: [PATCH 08/14] =fix whitespace
---
base/special/gamma.jl | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index a322d92..102b167 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -108,7 +108,7 @@ function lgamma(z::Complex{Float64})
-2.2315475845357937976132853e-04,9.9457512781808533714662972e-05,
-4.4926236738133141700224489e-05,2.0507212775670691553131246e-05)
end
- # use recurrence relation lgamma(z) = lgamma(z+1) - log(z )
+ # use recurrence relation lgamma(z) = lgamma(z+1) - log(z)
# to shift to x > 7 for asymptotic series
shiftprod = Complex(x,yabs)
x += 1
From 5a11d1d1ba89f47bc8c1a775b56e55312f967206 Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Mon, 31 Oct 2016 20:02:48 +0800
Subject: [PATCH 09/14] =More comments explaining convert strat
---
base/special/gamma.jl | 16 +++++++++++-----
1 file changed, 11 insertions(+), 5 deletions(-)
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 102b167..253f5a8 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -610,16 +610,22 @@ end
# Converting types that we can convert, and not ones we can not
+# Float16, and Float32 and their Complex equivalents can be converted to Float64
+# and results converted back. Similar for BitIntegers (eg Int32), but no converting back.
+# Otherwise, we need to make things use their own `float` converting methods
+# and in those cases, we do not convert back either as we assume
+# they also implement their own versions of the functions, with the correct return types.
+# Otherwise, if they do not implement their version of the functions we
+# manually throw a `MethodError`.
+# This case occurs, when calling `float` on a type does not change its type,
+# as it is already a `float`, and would have hit own method, if one had existed.
+
+
# If we really cared about single precision, we could make a faster
# Float32 version by truncating the Stirling series at a smaller cutoff.
# and if we really cared about half precision, we could make a faster
# Float16 version, by using a precomputed table look-up.
-# Float16, and Float32 and their Complex equivalents can be converted to Float64
-# and results converted back. Similar for BitIntegers (eg Int32), but no converting back
-# Otherwise, we need to make things use their own `float` converting methods
-# and in those cases, we do not convert back either as we assume
-# they also implement there own versions of the functions
ofpromotedtype(as::Tuple, c) = convert(promote_type(typeof.(as)...), c)
From f0573352016d8729c9c74fee156745a56d3a51dd Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Tue, 1 Nov 2016 22:56:25 +0800
Subject: [PATCH 10/14] =whitespace [ci skip]
---
base/special/gamma.jl | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 253f5a8..c38bee0e 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -615,7 +615,7 @@ end
# Otherwise, we need to make things use their own `float` converting methods
# and in those cases, we do not convert back either as we assume
# they also implement their own versions of the functions, with the correct return types.
-# Otherwise, if they do not implement their version of the functions we
+# Otherwise, if they do not implement their version of the functions we
# manually throw a `MethodError`.
# This case occurs, when calling `float` on a type does not change its type,
# as it is already a `float`, and would have hit own method, if one had existed.
From 6334ddac892f5c998cf8d8a2421a9fe3c6b0b422 Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Mon, 21 Nov 2016 11:28:30 +0800
Subject: [PATCH 11/14] =remove unrequired special casing for bitintegers
---
base/special/gamma.jl | 75 ++++++++++++++++++++-------------------------------
test/math.jl | 11 +++++---
2 files changed, 36 insertions(+), 50 deletions(-)
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index c38bee0e..5adc4d8 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -344,10 +344,7 @@ 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, ComplexOrReal{Float64}},
- z::ComplexOrReal{Float64})
+function zeta(s::ComplexOrReal{Float64}, z::ComplexOrReal{Float64})
ζ = zero(promote_type(typeof(s), typeof(z)))
(z == 1 || z == 0) && return oftype(ζ, zeta(s))
@@ -611,10 +608,11 @@ end
# Converting types that we can convert, and not ones we can not
# Float16, and Float32 and their Complex equivalents can be converted to Float64
-# and results converted back. Similar for BitIntegers (eg Int32), but no converting back.
+# and results converted back.
# Otherwise, we need to make things use their own `float` converting methods
# and in those cases, we do not convert back either as we assume
# they also implement their own versions of the functions, with the correct return types.
+# This is the case for BitIntegers (which become `Float64` when `float`ed).
# Otherwise, if they do not implement their version of the functions we
# manually throw a `MethodError`.
# This case occurs, when calling `float` on a type does not change its type,
@@ -627,13 +625,7 @@ end
# Float16 version, by using a precomputed table look-up.
-
-ofpromotedtype(as::Tuple, c) = convert(promote_type(typeof.(as)...), c)
-ComplexOrRealUnion(TS...) = Union{(ComplexOrReal{T} for T in TS)...}
-
-const types_le_Float64 = (Float16, Float32, Float64, Base.BitInteger.types...)
-
-for T in types_le_Float64
+for T in (Float16, Float32, Float64)
@eval f64(x::Complex{$T}) = Complex128(x)
@eval f64(x::$T) = Float64(x)
end
@@ -641,8 +633,9 @@ end
for f in (:digamma, :trigamma, :zeta, :eta, :invdigamma)
@eval begin
- $f(z::ComplexOrRealUnion(Base.BitInteger.types...)) = $f(f64(z))
- $f(z::ComplexOrRealUnion(Float16,Float32)) = oftype(z, $f(f64(z)))
+ function $f(z::Union{ComplexOrReal{Float16}, ComplexOrReal{Float32}})
+ oftype(z, $f(f64(z)))
+ end
function $f(z::Number)
x = float(z)
@@ -653,41 +646,19 @@ for f in (:digamma, :trigamma, :zeta, :eta, :invdigamma)
end
end
-function polygamma(m::Integer, z::ComplexOrRealUnion(Float16,Float32))
- oftype(z, polygamma(m, f64(z)))
-end
-for T1 in types_le_Float64, T2 in types_le_Float64
- if (T1 == T2 == Float64) || (T1 == Int && T2 == Float64)
- continue # Avoid redefining base definition
- # However this skips `zeta(::Complex{Int}, ::ComplexOrReal{Float64})`
- # so that will need to be added back after
- end
-
- if T1<:Integer && T2<:Integer
- @eval function zeta(s::ComplexOrReal{$T1}, z::ComplexOrReal{$T2})
- zeta(f64(s), f64(z)) #Do not promote down to Integers
- end
- else
- @eval function zeta(s::ComplexOrReal{$T1}, z::ComplexOrReal{$T2})
- ofpromotedtype((s, z), zeta(f64(s), f64(z)))
- end
- end
-end
+for T1 in (Float16, Float32, Float64), T2 in (Float16, Float32, Float64)
+ (T1 == T2 == Float64) && continue # Avoid redefining base definition
-# this is the one definition that is skipped
-function zeta(s::Complex{Int}, z::ComplexOrReal{Float64})::Complex{Float64}
- zeta(f64(s), f64(z))
-end
+ @eval function zeta(s::ComplexOrReal{$T1}, z::ComplexOrReal{$T2})
+ ζ = zeta(f64(s), f64(z))
-function zeta(s::Integer, z::Number)
- t = Int(s) # One could worry here about converting a BigInteger into a Int32/Int64
- x = float(z)
- if typeof(t) === typeof(s) && typeof(x) === typeof(z)
- # There is nothing to fallback to, since this didn't work
- throw(MethodError(zeta,(s,z)))
- end
- zeta(t, x)
+ if typeof(s) <: Complex || typeof(z) <: Complex
+ convert(Complex{$T2}, ζ)
+ else
+ convert(typeof(z), ζ)
+ end
+ end
end
@@ -701,6 +672,18 @@ function zeta(s::Number, z::Number)
zeta(t, x)
end
+# It is safe to convert `s` to a float as underflow will occur before precision issues
+zeta(s::Integer, z::Number) = zeta(oftype(z, s), z)
+function zeta(s::Integer, z::Integer)
+ x=float(z)
+ zeta(oftype(x, s), x)
+end
+
+
+function polygamma(m::Integer, z::Union{ComplexOrReal{Float16}, ComplexOrReal{Float32}})
+ oftype(z, polygamma(m, f64(z)))
+end
+
function polygamma(m::Integer, z::Number)
x = float(z)
diff --git a/test/math.jl b/test/math.jl
index bbc6b8f..1599302 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -939,7 +939,6 @@ end
@test Base.Math.f64(complex(1f0,1f0)) == complex(1.0, 1.0)
@test Base.Math.f64(1f0) == 1.0
-@test Base.Math.ofpromotedtype((complex(1f0, 1f0), 1.0), 5.0) == complex(5.0)
# no domain error is thrown for negative values
@test invoke(cbrt, Tuple{AbstractFloat}, -1.0) == -1.0
@@ -962,8 +961,12 @@ end
@test_throws MethodError zeta(1.0,big"2.0")
@test_throws MethodError zeta(big"1.0",2.0)
-@test typeof(zeta(complex(1),2.0)) == Complex{Float64}
+
+@test typeof(zeta(complex(1), 2.0)) == Complex{Float64}
@test typeof(polygamma(3, 0x2)) == Float64
@test typeof(polygamma(big"3", 2f0)) == Float32
-@test typeof(zeta(big"1",2.0)) == Float64 #Is this really desirable behavour?
-
+@test typeof(zeta(big"1", 2.0)) == Float64
+@test typeof(zeta(big"1", 2f0)) == Float32
+@test typeof(zeta(2f0, 2f0)) == Float32
+@test typeof(zeta(2f0, complex(2f0,0f0))) == Complex{Float32}
+@test typeof(zeta(complex(1,1), 2f0)) == Complex{Float32}
From e81314ffa0364d8e0926b8febc3fcaf8abed4755 Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Thu, 24 Nov 2016 10:28:14 +0800
Subject: [PATCH 12/14] =remove last of the special casing for Integers
---
base/special/gamma.jl | 32 +++++++++++---------------------
test/math.jl | 10 +++++-----
2 files changed, 16 insertions(+), 26 deletions(-)
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 5adc4d8..d03185f 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -468,7 +468,9 @@ function polygamma(m::Integer, z::ComplexOrReal{Float64})
# constants. We throw a DomainError() since the definition is unclear.
real(m) < 0 && throw(DomainError())
- s = m+1
+ s = Float64(m+1)
+ # It is safe to convert any integer (including `BigInt`) to a float here
+ # as underflow occurs before precision issues.
if real(z) <= 0 # reflection formula
(zeta(s, 1-z) + signflip(m, cotderiv(m,z))) * (-gamma(s))
else
@@ -633,14 +635,14 @@ end
for f in (:digamma, :trigamma, :zeta, :eta, :invdigamma)
@eval begin
- function $f(z::Union{ComplexOrReal{Float16}, ComplexOrReal{Float32}})
- oftype(z, $f(f64(z)))
- end
+ function $f(z::Union{ComplexOrReal{Float16}, ComplexOrReal{Float32}})
+ oftype(z, $f(f64(z)))
+ end
function $f(z::Number)
x = float(z)
typeof(x) === typeof(z) && throw(MethodError($f, (z,)))
- # There is nothing to fallback to, since this didn't change the argument types
+ # There is nothing to fallback to, as this didn't change the argument types
$f(x)
end
end
@@ -650,15 +652,10 @@ end
for T1 in (Float16, Float32, Float64), T2 in (Float16, Float32, Float64)
(T1 == T2 == Float64) && continue # Avoid redefining base definition
- @eval function zeta(s::ComplexOrReal{$T1}, z::ComplexOrReal{$T2})
- ζ = zeta(f64(s), f64(z))
-
- if typeof(s) <: Complex || typeof(z) <: Complex
- convert(Complex{$T2}, ζ)
- else
- convert(typeof(z), ζ)
- end
- end
+ @eval function zeta(s::ComplexOrReal{$T1}, z::ComplexOrReal{$T2})
+ ζ = zeta(f64(s), f64(z))
+ convert(promote_type(typeof(s), typeof(z)), ζ)
+ end
end
@@ -672,13 +669,6 @@ function zeta(s::Number, z::Number)
zeta(t, x)
end
-# It is safe to convert `s` to a float as underflow will occur before precision issues
-zeta(s::Integer, z::Number) = zeta(oftype(z, s), z)
-function zeta(s::Integer, z::Integer)
- x=float(z)
- zeta(oftype(x, s), x)
-end
-
function polygamma(m::Integer, z::Union{ComplexOrReal{Float16}, ComplexOrReal{Float32}})
oftype(z, polygamma(m, f64(z)))
diff --git a/test/math.jl b/test/math.jl
index 1599302..9e4f8ce 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -960,13 +960,13 @@ end
@test_throws MethodError zeta(1.0,big"2")
@test_throws MethodError zeta(1.0,big"2.0")
@test_throws MethodError zeta(big"1.0",2.0)
+@test_throws MethodError zeta(big"1",2.0)
-@test typeof(zeta(complex(1), 2.0)) == Complex{Float64}
@test typeof(polygamma(3, 0x2)) == Float64
@test typeof(polygamma(big"3", 2f0)) == Float32
-@test typeof(zeta(big"1", 2.0)) == Float64
-@test typeof(zeta(big"1", 2f0)) == Float32
-@test typeof(zeta(2f0, 2f0)) == Float32
+@test typeof(zeta(1, 2.0)) == Float64
+@test typeof(zeta(1, 2f0)) == Float64 # BitIntegers result in Float64 returns
@test typeof(zeta(2f0, complex(2f0,0f0))) == Complex{Float32}
-@test typeof(zeta(complex(1,1), 2f0)) == Complex{Float32}
+@test typeof(zeta(complex(1,1), 2f0)) == Complex{Float64}
+@test typeof(zeta(complex(1), 2.0)) == Complex{Float64}
From 0f62f12dc946747a4e116eb5dd0320ca347059db Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Mon, 28 Nov 2016 13:16:08 +0800
Subject: [PATCH 13/14] =fix formatting in docstring
---
base/special/gamma.jl | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index d03185f..e59006f 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -448,7 +448,7 @@ end
polygamma(m, x)
Compute the polygamma function of order `m` of argument `x`
-(the `(m+1)th` derivative of the logarithm of `gamma(x)`)
+(the `(m+1)`th derivative of the logarithm of `gamma(x)`)
"""
function polygamma(m::Integer, z::ComplexOrReal{Float64})
m == 0 && return digamma(z)
From 734ca5426ae533a1d6f8f2f4884ea433655061e6 Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Sun, 5 Feb 2017 07:57:59 +0800
Subject: [PATCH 14/14] =remove whitespace
---
test/math.jl | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/test/math.jl b/test/math.jl
index ff5e688..e12bf9f 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -242,7 +242,7 @@ end
@test hypot(T(Inf), T(x)) === T(Inf)
@test hypot(T(Inf), T(NaN)) === T(Inf)
@test isnan(hypot(T(x), T(NaN)))
-
+
@test invoke(cbrt, Tuple{AbstractFloat}, -1.0) == -1.0 # no domain error is thrown for negative values
end
end
@@ -1001,7 +1001,7 @@ end
@testset "issue #17474" begin
@test Base.Math.f64(complex(1f0,1f0)) == complex(1.0, 1.0)
@test Base.Math.f64(1f0) == 1.0
-
+
@test typeof(eta(big"2")) == BigFloat
@test typeof(zeta(big"2")) == BigFloat
@test typeof(digamma(big"2")) == BigFloat
@@ -1029,7 +1029,7 @@ end
@test typeof(zeta(complex(1,1), 2f0)) == Complex{Float64}
@test typeof(zeta(complex(1), 2.0)) == Complex{Float64}
end
-
+
@testset "promote Float16 irrational #15359" begin
@test typeof(Float16(.5) * pi) == Float16
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment