-
-
Save anonymous/84de6e5d71d989c549cb20671929302c to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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