From 3e8fb7a78727def0eea4f260736ce7e86df8607f Mon Sep 17 00:00:00 2001 From: Alex Arslan Date: Fri, 27 Jan 2017 15:55:11 -0500 Subject: [PATCH] Remove some special functions from Base --- base/deprecated.jl | 94 +--- base/docs/helpdb/Base.jl | 66 --- base/exports.jl | 39 -- base/math.jl | 35 +- base/mpfr.jl | 64 +-- base/replutil.jl | 2 +- base/special/bessel.jl | 497 ---------------------- base/special/erf.jl | 217 ---------- base/special/gamma.jl | 475 --------------------- doc/src/manual/mathematical-operations.md | 32 -- doc/src/stdlib/math.md | 40 -- test/bigint.jl | 2 - test/math.jl | 408 +----------------- test/mpfr.jl | 15 - test/random.jl | 3 +- 15 files changed, 33 insertions(+), 1956 deletions(-) delete mode 100644 base/special/bessel.jl delete mode 100644 base/special/erf.jl diff --git a/base/deprecated.jl b/base/deprecated.jl index 9db39ec2b8edaa..5e78d033405a03 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -226,14 +226,9 @@ for f in ( # base/special/log.jl :log, :log1p, # base/special/gamma.jl - :gamma, :lfact, :digamma, :trigamma, :zeta, :eta, - # base/special/erf.jl - :erfcx, :erfi, :dawson, - # base/special/bessel.jl - :airyai, :airyaiprime, :airybi, :airybiprime, - :besselj0, :besselj1, :bessely0, :bessely1, + :gamma, :lfact, # base/math.jl - :cbrt, :sinh, :cosh, :tanh, :atan, :asinh, :exp, :erf, :erfc, :exp2, + :cbrt, :sinh, :cosh, :tanh, :atan, :asinh, :exp, :exp2, :expm1, :exp10, :sin, :cos, :tan, :asin, :acos, :acosh, :atanh, #=:log,=# :log2, :log10, :lgamma, #=:log1p,=# :sqrt, # base/floatfuncs.jl @@ -252,8 +247,6 @@ for f in ( :acos_fast, :acosh_fast, :angle_fast, :asin_fast, :asinh_fast, @eval FastMath Base.@dep_vectorize_1arg Number $f end for f in ( - :invdigamma, # base/special/gamma.jl - :erfinc, :erfcinv, # base/special/erf.jl :trunc, :floor, :ceil, :round, # base/floatfuncs.jl :rad2deg, :deg2rad, :exponent, :significand, # base/math.jl :sind, :cosd, :tand, :asind, :acosd, :atand, :asecd, :acscd, :acotd, # base/special/trig.jl @@ -292,11 +285,7 @@ end # Deprecate @vectorize_2arg-vectorized functions from... for f in ( # base/special/gamma.jl - :polygamma, :zeta, :beta, :lbeta, - # base/special/bessel.jl - :besseli, :besselix, :besselj, :besseljx, - :besselk, :besselkx, :bessely, :besselyx, :besselh, - :besselhx, :hankelh1, :hankelh2, :hankelh1x, :hankelh2x, + :beta, :lbeta, # base/math.jl :log, :hypot, :atan2, ) @@ -672,65 +661,6 @@ end # Deprecate isimag (#19947). @deprecate isimag(z::Number) iszero(real(z)) -@deprecate airy(z::Number) airyai(z) -@deprecate airyx(z::Number) airyaix(z) -@deprecate airyprime(z::Number) airyaiprime(z) -@deprecate airy{T<:Number}(x::AbstractArray{T}) airyai.(x) -@deprecate airyx{T<:Number}(x::AbstractArray{T}) airyaix.(x) -@deprecate airyprime{T<:Number}(x::AbstractArray{T}) airyprime.(x) - -function _airy(k::Integer, z::Complex128) - depwarn("`airy(k,x)` is deprecated, use `airyai(x)`, `airyaiprime(x)`, `airybi(x)` or `airybiprime(x)` instead.",:airy) - id = Int32(k==1 || k==3) - if k == 0 || k == 1 - return Base.Math._airy(z, id, Int32(1)) - elseif k == 2 || k == 3 - return Base.Math._biry(z, id, Int32(1)) - else - throw(ArgumentError("k must be between 0 and 3")) - end -end -function _airyx(k::Integer, z::Complex128) - depwarn("`airyx(k,x)` is deprecated, use `airyaix(x)`, `airyaiprimex(x)`, `airybix(x)` or `airybiprimex(x)` instead.",:airyx) - id = Int32(k==1 || k==3) - if k == 0 || k == 1 - return Base.Math._airy(z, id, Int32(2)) - elseif k == 2 || k == 3 - return Base.Math._biry(z, id, Int32(2)) - else - throw(ArgumentError("k must be between 0 and 3")) - end -end - -for afn in (:airy,:airyx) - _afn = Symbol("_"*string(afn)) - suf = string(afn)[5:end] - @eval begin - function $afn(k::Integer, z::Complex128) - afn = $(QuoteNode(afn)) - suf = $(QuoteNode(suf)) - depwarn("`$afn(k,x)` is deprecated, use `airyai$suf(x)`, `airyaiprime$suf(x)`, `airybi$suf(x)` or `airybiprime$suf(x)` instead.",$(QuoteNode(afn))) - $_afn(k,z) - end - - $afn(k::Integer, z::Complex) = $afn(k, float(z)) - $afn{T<:AbstractFloat}(k::Integer, z::Complex{T}) = throw(MethodError($afn,(k,z))) - $afn(k::Integer, z::Complex64) = Complex64($afn(k, Complex128(z))) - $afn(k::Integer, x::Real) = $afn(k, float(x)) - $afn(k::Integer, x::AbstractFloat) = real($afn(k, complex(x))) - - function $afn{T<:Number}(k::Number, x::AbstractArray{T}) - $afn.(k,x) - end - function $afn{S<:Number}(k::AbstractArray{S}, x::Number) - $afn.(k,x) - end - function $afn{S<:Number,T<:Number}(k::AbstractArray{S}, x::AbstractArray{T}) - $afn.(k,x) - end - end -end - # Deprecate vectorized xor in favor of compact broadcast syntax @deprecate xor(a::Bool, B::BitArray) xor.(a, B) @deprecate xor(A::BitArray, b::Bool) xor.(A, b) @@ -1227,6 +1157,24 @@ for name in ("alnum", "alpha", "cntrl", "digit", "number", "graph", @eval @deprecate ($f)(s::AbstractString) all($f, s) end +# Special functions have been moved to a package +for f in (:airyai, :airyaiprime, :airybi, :airybiprime, :airyaix, :airyaiprimex, :airybix, :airybiprimex, + :besselh, :besselhx, :besseli, :besselix, :besselj, :besselj0, :besselj1, :besseljx, :besselk, + :besselkx, :bessely, :bessely0, :bessely1, :besselyx, + :dawson, :erf, :erfc, :erfcinv, :erfcx, :erfi, :erfinv, + :eta, :zeta, :digamma, :invdigamma, :polygamma, :trigamma, + :hankelh1, :hankelh1x, :hankelh2, :hankelh2x, + :airy, :airyx, :airyprime) + @eval begin + function $f(args...; kwargs...) + error(string($f, args, " has been moved to the package SpecialFunctions.jl.\n", + "Run Pkg.add(\"SpecialFunctions\") to install SpecialFunctions on Julia v0.6 and later,\n", + "and then run `using SpecialFunctions`.")) + end + export $f + end +end + # END 0.6 deprecations # BEGIN 1.0 deprecations diff --git a/base/docs/helpdb/Base.jl b/base/docs/helpdb/Base.jl index b79496acbfee0b..a04d85de8d2e1a 100644 --- a/base/docs/helpdb/Base.jl +++ b/base/docs/helpdb/Base.jl @@ -304,13 +304,6 @@ This would create a 25-by-30000 `BitArray`, linked to the file associated with s """ Mmap.mmap(io, ::BitArray, dims = ?, offset = ?) -""" - bessely0(x) - -Bessel function of the second kind of order 0, ``Y_0(x)``. -""" -bessely0 - """ filter!(function, collection) @@ -647,13 +640,6 @@ for use in `Mmap.mmap`. Used by `SharedArray` for creating shared memory arrays. """ Mmap.Anonymous -""" - erfi(x) - -Compute the imaginary error function of `x`, defined by ``-i \\operatorname{erf}(ix)``. -""" -erfi - """ floor([T,] x, [digits, [base]]) @@ -695,13 +681,6 @@ The item or field is not defined for the given object. """ UndefRefError -""" - bessely1(x) - -Bessel function of the second kind of order 1, ``Y_1(x)``. -""" -bessely1 - """ append!(collection, collection2) -> collection. @@ -798,13 +777,6 @@ julia> getfield(a, :num) """ getfield -""" - besselj1(x) - -Bessel function of the first kind of order 1, ``J_1(x)``. -""" -besselj1 - """ select!(v, k, [by=,] [lt=,] [rev=false]) @@ -933,13 +905,6 @@ behavior, including program corruption or segfaults, at any later time. """ unsafe_convert -""" - erfinv(x) - -Compute the inverse error function of a real `x`, defined by ``\\operatorname{erf}(\\operatorname{erfinv}(x)) = x``. -""" -erfinv - """ seek(s, pos) @@ -947,21 +912,6 @@ Seek a stream to the given position. """ seek -""" - besselj0(x) - -Bessel function of the first kind of order 0, ``J_0(x)``. -""" -besselj0 - -""" - erfcinv(x) - -Compute the inverse error complementary function of a real `x`, defined by -``\\operatorname{erfc}(\\operatorname{erfcinv}(x)) = x``. -""" -erfcinv - """ popdisplay() popdisplay(d::Display) @@ -1605,14 +1555,6 @@ Equivalent to [`readdlm`](@ref) with `delim` set to comma, and type optionally d """ readcsv -""" - erfcx(x) - -Compute the scaled complementary error function of `x`, defined by ``e^{x^2} \\operatorname{erfc}(x)``. -Note also that ``\\operatorname{erfcx}(-ix)`` computes the Faddeeva function ``w(x)``. -""" -erfcx - """ UndefVarError(var::Symbol) @@ -2482,11 +2424,3 @@ seekend Integer division was attempted with a denominator value of 0. """ DivideError - -""" - dawson(x) - -Compute the Dawson function (scaled imaginary error function) of `x`, defined by -``\\frac{\\sqrt{\\pi}}{2} e^{-x^2} \\operatorname{erfi}(x)``. -""" -dawson diff --git a/base/exports.jl b/base/exports.jl index 78b9ee477cc2be..6dc7b4ed413e0f 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -309,19 +309,11 @@ export csc, cscd, csch, - dawson, deg2rad, denominator, - digamma, div, divrem, eps, - erf, - erfc, - erfcinv, - erfcx, - erfi, - erfinv, exp, exp10, exp2, @@ -345,7 +337,6 @@ export hypot, imag, inv, - invdigamma, invmod, isapprox, iseven, @@ -416,7 +407,6 @@ export tanh, trailing_ones, trailing_zeros, - trigamma, trunc, unsafe_trunc, typemax, @@ -430,37 +420,8 @@ export ≉, # specfun - airyai, - airyaiprime, - airybi, - airybiprime, - airyaix, - airyaiprimex, - airybix, - airybiprimex, - besselh, - besselhx, - besseli, - besselix, - besselj, - besselj0, - besselj1, - besseljx, - besselk, - besselkx, - bessely, - bessely0, - bessely1, - besselyx, beta, - eta, - hankelh1, - hankelh1x, - hankelh2, - hankelh2x, lbeta, - polygamma, - zeta, # arrays broadcast!, diff --git a/base/math.jl b/base/math.jl index 247196c1290d02..ea443d81cd4750 100644 --- a/base/math.jl +++ b/base/math.jl @@ -10,18 +10,10 @@ export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan, acosd, acotd, acscd, asecd, asind, atand, atan2, rad2deg, deg2rad, log, log2, log10, log1p, exponent, exp, exp2, exp10, expm1, - cbrt, sqrt, erf, erfc, erfcx, erfi, dawson, - significand, + cbrt, sqrt, significand, lgamma, hypot, gamma, lfact, max, min, minmax, ldexp, frexp, clamp, clamp!, modf, ^, mod2pi, rem2pi, - airyai, airyaiprime, airybi, airybiprime, - airyaix, airyaiprimex, airybix, airybiprimex, - besselj0, besselj1, besselj, besseljx, - bessely0, bessely1, bessely, besselyx, - hankelh1, hankelh2, hankelh1x, hankelh2x, - besseli, besselix, besselk, besselkx, besselh, besselhx, - beta, lbeta, eta, zeta, polygamma, invdigamma, digamma, trigamma, - erfinv, erfcinv, @evalpoly + beta, lbeta, @evalpoly import Base: log, exp, sin, cos, tan, sinh, cosh, tanh, asin, acos, atan, asinh, acosh, atanh, sqrt, log2, log10, @@ -224,28 +216,13 @@ Compute the inverse hyperbolic sine of `x`. """ asinh(x) -""" - erf(x) - -Compute the error function of `x`, defined by ``\\frac{2}{\\sqrt{\\pi}} \\int_0^x e^{-t^2} dt`` -for arbitrary complex `x`. -""" -erf(x) - -""" - erfc(x) - -Compute the complementary error function of `x`, defined by ``1 - \\operatorname{erf}(x)``. -""" -erfc(x) - """ expm1(x) Accurately compute ``e^x-1``. """ expm1(x) -for f in (:cbrt, :sinh, :cosh, :tanh, :atan, :asinh, :erf, :erfc, :exp2, :expm1) +for f in (:cbrt, :sinh, :cosh, :tanh, :atan, :asinh, :exp2, :expm1) @eval begin ($f)(x::Float64) = ccall(($(string(f)),libm), Float64, (Float64,), x) ($f)(x::Float32) = ccall(($(string(f,"f")),libm), Float32, (Float32,), x) @@ -284,7 +261,7 @@ julia> exp2(5) ``` """ exp2(x::AbstractFloat) = 2^x -for f in (:sinh, :cosh, :tanh, :atan, :asinh, :exp, :erf, :erfc, :expm1) +for f in (:sinh, :cosh, :tanh, :atan, :asinh, :exp, :expm1) @eval ($f)(x::AbstractFloat) = error("not implemented for ", typeof(x)) end @@ -918,7 +895,7 @@ muladd(x,y,z) = x*y+z # Float16 definitions for func in (:sin,:cos,:tan,:asin,:acos,:atan,:sinh,:cosh,:tanh,:asinh,:acosh, - :atanh,:exp,:log,:log2,:log10,:sqrt,:lgamma,:log1p,:erf,:erfc) + :atanh,:exp,:log,:log2,:log10,:sqrt,:lgamma,:log1p) @eval begin $func(a::Float16) = Float16($func(Float32(a))) $func(a::Complex32) = Complex32($func(Complex64(a))) @@ -935,8 +912,6 @@ cbrt(a::Float16) = Float16(cbrt(Float32(a))) # More special functions include("special/trig.jl") -include("special/bessel.jl") -include("special/erf.jl") include("special/gamma.jl") module JuliaLibm diff --git a/base/mpfr.jl b/base/mpfr.jl index 024860ef1ef399..603e8fa025cf1d 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -8,13 +8,12 @@ export big_str import - Base: (*), +, -, /, <, <=, ==, >, >=, ^, besselj, besselj0, besselj1, bessely, - bessely0, bessely1, ceil, cmp, convert, copysign, div, + Base: (*), +, -, /, <, <=, ==, >, >=, ^, ceil, cmp, convert, copysign, div, exp, exp2, exponent, factorial, floor, fma, hypot, isinteger, isfinite, isinf, isnan, ldexp, log, log2, log10, max, min, mod, modf, nextfloat, prevfloat, promote_rule, rem, rem2pi, round, show, sum, sqrt, string, print, trunc, precision, exp10, expm1, - gamma, lgamma, digamma, erf, erfc, zeta, eta, log1p, airyai, + gamma, lgamma, log1p, eps, signbit, sin, cos, tan, sec, csc, cot, acos, asin, atan, cosh, sinh, tanh, sech, csch, coth, acosh, asinh, atanh, atan2, cbrt, typemax, typemin, unsafe_trunc, realmin, realmax, rounding, @@ -505,8 +504,7 @@ end ^(x::BigFloat, y::Integer) = typemin(Clong) <= y <= typemax(Clong) ? x^Clong(y) : x^BigInt(y) ^(x::BigFloat, y::Unsigned) = typemin(Culong) <= y <= typemax(Culong) ? x^Culong(y) : x^BigInt(y) -for f in (:exp, :exp2, :exp10, :expm1, :digamma, :erf, :erfc, :zeta, - :cosh,:sinh,:tanh,:sech,:csch,:coth, :cbrt) +for f in (:exp, :exp2, :exp10, :expm1, :cosh, :sinh, :tanh, :sech, :csch, :coth, :cbrt) @eval function $f(x::BigFloat) z = BigFloat() ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[]) @@ -522,17 +520,6 @@ function big_ln2() return c end -function eta(x::BigFloat) - x == 1 && return big_ln2() - return -zeta(x) * expm1(big_ln2()*(1-x)) -end - -function airyai(x::BigFloat) - z = BigFloat() - ccall((:mpfr_ai, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[]) - return z -end - function ldexp(x::BigFloat, n::Clong) z = BigFloat() ccall((:mpfr_mul_2si, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Clong, Int32), &z, &x, n, ROUNDING_MODE[]) @@ -547,51 +534,6 @@ ldexp(x::BigFloat, n::ClongMax) = ldexp(x, convert(Clong, n)) ldexp(x::BigFloat, n::CulongMax) = ldexp(x, convert(Culong, n)) ldexp(x::BigFloat, n::Integer) = x*exp2(BigFloat(n)) -function besselj0(x::BigFloat) - z = BigFloat() - ccall((:mpfr_j0, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[]) - return z -end - -function besselj1(x::BigFloat) - z = BigFloat() - ccall((:mpfr_j1, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[]) - return z -end - -function besselj(n::Integer, x::BigFloat) - z = BigFloat() - ccall((:mpfr_jn, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, n, &x, ROUNDING_MODE[]) - return z -end - -function bessely0(x::BigFloat) - if x < 0 - throw(DomainError()) - end - z = BigFloat() - ccall((:mpfr_y0, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[]) - return z -end - -function bessely1(x::BigFloat) - if x < 0 - throw(DomainError()) - end - z = BigFloat() - ccall((:mpfr_y1, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, ROUNDING_MODE[]) - return z -end - -function bessely(n::Integer, x::BigFloat) - if x < 0 - throw(DomainError()) - end - z = BigFloat() - ccall((:mpfr_yn, :libmpfr), Int32, (Ptr{BigFloat}, Clong, Ptr{BigFloat}, Int32), &z, n, &x, ROUNDING_MODE[]) - return z -end - function factorial(x::BigFloat) if x < 0 || !isinteger(x) throw(DomainError()) diff --git a/base/replutil.jl b/base/replutil.jl index 78f95cbe2951d6..ef2cbd466a9fe0 100644 --- a/base/replutil.jl +++ b/base/replutil.jl @@ -223,7 +223,7 @@ function showerror(io::IO, ex::DomainError, bt; backtrace=true) if !code.from_c if code.func == :nan_dom_err continue - elseif code.func in (:log, :log2, :log10, :sqrt) # TODO add :besselj, :besseli, :bessely, :besselk + elseif code.func in (:log, :log2, :log10, :sqrt) print(io, "\n$(code.func) will only return a complex result if called ", "with a complex argument. Try $(string(code.func))(complex(x)).") elseif (code.func == :^ && code.file == Symbol("intfuncs.jl")) || diff --git a/base/special/bessel.jl b/base/special/bessel.jl deleted file mode 100644 index 17e9b905af5341..00000000000000 --- a/base/special/bessel.jl +++ /dev/null @@ -1,497 +0,0 @@ -# This file is a part of Julia. License is MIT: http://julialang.org/license - -type AmosException <: Exception - info::Int32 -end - -## Airy functions -let - const ai::Array{Float64,1} = Array{Float64}(2) - const ae::Array{Int32,1} = Array{Int32}(2) - global _airy, _biry - function _airy(z::Complex128, id::Int32, kode::Int32) - ccall((:zairy_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), - &real(z), &imag(z), - &id, &kode, - pointer(ai,1), pointer(ai,2), - pointer(ae,1), pointer(ae,2)) - if ae[2] == 0 || ae[2] == 3 - return complex(ai[1],ai[2]) - else - throw(AmosException(ae[2])) - end - end - function _biry(z::Complex128, id::Int32, kode::Int32) - ccall((:zbiry_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}), - &real(z), &imag(z), - &id, &kode, - pointer(ai,1), pointer(ai,2), - pointer(ae,1)) - if ae[1] == 0 || ae[1] == 3 # ignore underflow - return complex(ai[1],ai[2]) - else - throw(AmosException(ae[2])) - end - end -end - - -""" - airyai(x) - -Airy function of the first kind ``\\operatorname{Ai}(x)``. -""" -function airyai end -airyai(z::Complex128) = _airy(z, Int32(0), Int32(1)) - -""" - airyaiprime(x) - -Derivative of the Airy function of the first kind ``\\operatorname{Ai}'(x)``. -""" -function airyaiprime end -airyaiprime(z::Complex128) = _airy(z, Int32(1), Int32(1)) - -""" - airybi(x) - -Airy function of the second kind ``\\operatorname{Bi}(x)``. -""" -function airybi end -airybi(z::Complex128) = _biry(z, Int32(0), Int32(1)) - -""" - airybiprime(x) - -Derivative of the Airy function of the second kind ``\\operatorname{Bi}'(x)``. -""" -function airybiprime end -airybiprime(z::Complex128) = _biry(z, Int32(1), Int32(1)) - -""" - airyaix(x) - -Scaled Airy function of the first kind ``\\operatorname{Ai}(x) e^{\\frac{2}{3} x -\\sqrt{x}}``. Throws [`DomainError`](@ref) for negative `Real` arguments. -""" -function airyaix end -airyaix(z::Complex128) = _airy(z, Int32(0), Int32(2)) - -""" - airyaiprimex(x) - -Scaled derivative of the Airy function of the first kind ``\\operatorname{Ai}'(x) -e^{\\frac{2}{3} x \\sqrt{x}}``. Throws [`DomainError`](@ref) for negative `Real` arguments. -""" -function airyaiprimex end -airyaiprimex(z::Complex128) = _airy(z, Int32(1), Int32(2)) - -""" - airybix(x) - -Scaled Airy function of the second kind ``\\operatorname{Bi}(x) e^{- \\left| \\operatorname{Re} \\left( \\frac{2}{3} x \\sqrt{x} \\right) \\right|}``. -""" -function airybix end -airybix(z::Complex128) = _biry(z, Int32(0), Int32(2)) - -""" - airybiprimex(x) - -Scaled derivative of the Airy function of the second kind ``\\operatorname{Bi}'(x) e^{- \\left| \\operatorname{Re} \\left( \\frac{2}{3} x \\sqrt{x} \\right) \\right|}``. -""" -function airybiprimex end -airybiprimex(z::Complex128) = _biry(z, Int32(1), Int32(2)) - -for afn in (:airyai, :airyaiprime, :airybi, :airybiprime, - :airyaix, :airyaiprimex, :airybix, :airybiprimex) - @eval begin - $afn(z::Complex) = $afn(float(z)) - $afn{T<:AbstractFloat}(z::Complex{T}) = throw(MethodError($afn,(z,))) - $afn(z::Complex64) = Complex64($afn(Complex128(z))) - end - if afn in (:airyaix, :airyaiprimex) - @eval $afn(x::Real) = x < 0 ? throw(DomainError()) : real($afn(complex(float(x)))) - else - @eval $afn(x::Real) = real($afn(complex(float(x)))) - end -end - -## Bessel functions - -# besselj0, besselj1, bessely0, bessely1 -for jy in ("j","y"), nu in (0,1) - jynu = Expr(:quote, Symbol(jy,nu)) - jynuf = Expr(:quote, Symbol(jy,nu,"f")) - bjynu = Symbol("bessel",jy,nu) - if jy == "y" - @eval begin - $bjynu(x::Float64) = nan_dom_err(ccall(($jynu,libm), Float64, (Float64,), x), x) - $bjynu(x::Float32) = nan_dom_err(ccall(($jynuf,libm), Float32, (Float32,), x), x) - end - else - @eval begin - $bjynu(x::Float64) = ccall(($jynu,libm), Float64, (Float64,), x) - $bjynu(x::Float32) = ccall(($jynuf,libm), Float32, (Float32,), x) - end - end - @eval begin - $bjynu(x::Real) = $bjynu(float(x)) - $bjynu(x::Complex) = $(Symbol("bessel",jy))($nu,x) - end -end - - -const cy = Array{Float64}(2) -const ae = Array{Int32}(2) -const wrk = Array{Float64}(2) - -function _besselh(nu::Float64, k::Int32, z::Complex128, kode::Int32) - ccall((:zbesh_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), - &real(z), &imag(z), &nu, &kode, &k, &1, - pointer(cy,1), pointer(cy,2), - pointer(ae,1), pointer(ae,2)) - if ae[2] == 0 || ae[2] == 3 - return complex(cy[1],cy[2]) - else - throw(AmosException(ae[2])) - end -end - -function _besseli(nu::Float64, z::Complex128, kode::Int32) - ccall((:zbesi_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), - &real(z), &imag(z), &nu, &kode, &1, - pointer(cy,1), pointer(cy,2), - pointer(ae,1), pointer(ae,2)) - if ae[2] == 0 || ae[2] == 3 - return complex(cy[1],cy[2]) - else - throw(AmosException(ae[2])) - end -end - -function _besselj(nu::Float64, z::Complex128, kode::Int32) - ccall((:zbesj_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), - &real(z), &imag(z), &nu, &kode, &1, - pointer(cy,1), pointer(cy,2), - pointer(ae,1), pointer(ae,2)) - if ae[2] == 0 || ae[2] == 3 - return complex(cy[1],cy[2]) - else - throw(AmosException(ae[2])) - end -end - -function _besselk(nu::Float64, z::Complex128, kode::Int32) - ccall((:zbesk_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}), - &real(z), &imag(z), &nu, &kode, &1, - pointer(cy,1), pointer(cy,2), - pointer(ae,1), pointer(ae,2)) - if ae[2] == 0 || ae[2] == 3 - return complex(cy[1],cy[2]) - else - throw(AmosException(ae[2])) - end -end - -function _bessely(nu::Float64, z::Complex128, kode::Int32) - ccall((:zbesy_,openspecfun), Void, - (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, - Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, - Ptr{Float64}, Ptr{Float64}, Ptr{Int32}), - &real(z), &imag(z), &nu, &kode, &1, - pointer(cy,1), pointer(cy,2), - pointer(ae,1), pointer(wrk,1), - pointer(wrk,2), pointer(ae,2)) - if ae[2] == 0 || ae[2] == 3 - return complex(cy[1],cy[2]) - else - throw(AmosException(ae[2])) - end -end - -""" - besselh(nu, [k=1,] x) - -Bessel function of the third kind of order `nu` (the Hankel function). `k` is either 1 or 2, -selecting [`hankelh1`](@ref) or [`hankelh2`](@ref), respectively. -`k` defaults to 1 if it is omitted. -(See also [`besselhx`](@ref) for an exponentially scaled variant.) -""" -function besselh end - -function besselh(nu::Float64, k::Integer, z::Complex128) - if nu < 0 - s = (k == 1) ? 1 : -1 - return _besselh(-nu,Int32(k),z,Int32(1)) * complex(cospi(nu),-s*sinpi(nu)) - end - return _besselh(nu,Int32(k),z,Int32(1)) -end - -""" - besselhx(nu, [k=1,] z) - -Compute the scaled Hankel function ``\\exp(∓iz) H_ν^{(k)}(z)``, where -``k`` is 1 or 2, ``H_ν^{(k)}(z)`` is `besselh(nu, k, z)`, and ``∓`` is -``-`` for ``k=1`` and ``+`` for ``k=2``. `k` defaults to 1 if it is omitted. - -The reason for this function is that ``H_ν^{(k)}(z)`` is asymptotically -proportional to ``\\exp(∓iz)/\\sqrt{z}`` for large ``|z|``, and so the -[`besselh`](@ref) function is susceptible to overflow or underflow -when `z` has a large imaginary part. The `besselhx` function cancels this -exponential factor (analytically), so it avoids these problems. -""" -function besselhx end - -function besselhx(nu::Float64, k::Integer, z::Complex128) - if nu < 0 - s = (k == 1) ? 1 : -1 - return _besselh(-nu,Int32(k),z,Int32(2)) * complex(cospi(nu),-s*sinpi(nu)) - end - return _besselh(nu,Int32(k),z,Int32(2)) -end - -function besseli(nu::Float64, z::Complex128) - if nu < 0 - return _besseli(-nu,z,Int32(1)) - 2_besselk(-nu,z,Int32(1))*sinpi(nu)/pi - else - return _besseli(nu,z,Int32(1)) - end -end - -function besselix(nu::Float64, z::Complex128) - if nu < 0 - return _besseli(-nu,z,Int32(2)) - 2_besselk(-nu,z,Int32(2))*exp(-abs(real(z))-z)*sinpi(nu)/pi - else - return _besseli(nu,z,Int32(2)) - end -end - -function besselj(nu::Float64, z::Complex128) - if nu < 0 - return _besselj(-nu,z,Int32(1))*cospi(nu) + _bessely(-nu,z,Int32(1))*sinpi(nu) - else - return _besselj(nu,z,Int32(1)) - end -end - -besselj(nu::Cint, x::Float64) = ccall((:jn, libm), Float64, (Cint, Float64), nu, x) -besselj(nu::Cint, x::Float32) = ccall((:jnf, libm), Float32, (Cint, Float32), nu, x) - - -function besseljx(nu::Float64, z::Complex128) - if nu < 0 - return _besselj(-nu,z,Int32(2))*cospi(nu) + _bessely(-nu,z,Int32(2))*sinpi(nu) - else - return _besselj(nu,z,Int32(2)) - end -end - -besselk(nu::Float64, z::Complex128) = _besselk(abs(nu), z, Int32(1)) - -besselkx(nu::Float64, z::Complex128) = _besselk(abs(nu), z, Int32(2)) - -function bessely(nu::Cint, x::Float64) - if x < 0 - throw(DomainError()) - end - ccall((:yn, libm), Float64, (Cint, Float64), nu, x) -end -function bessely(nu::Cint, x::Float32) - if x < 0 - throw(DomainError()) - end - ccall((:ynf, libm), Float32, (Cint, Float32), nu, x) -end - -function bessely(nu::Float64, z::Complex128) - if nu < 0 - return _bessely(-nu,z,Int32(1))*cospi(nu) - _besselj(-nu,z,Int32(1))*sinpi(nu) - else - return _bessely(nu,z,Int32(1)) - end -end - -function besselyx(nu::Float64, z::Complex128) - if nu < 0 - return _bessely(-nu,z,Int32(2))*cospi(nu) - _besselj(-nu,z,Int32(2))*sinpi(nu) - else - return _bessely(nu,z,Int32(2)) - end -end - -""" - besseli(nu, x) - -Modified Bessel function of the first kind of order `nu`, ``I_\\nu(x)``. -""" -function besseli(nu::Real, x::AbstractFloat) - if x < 0 && !isinteger(nu) - throw(DomainError()) - end - real(besseli(float(nu), complex(x))) -end - -""" - besselix(nu, x) - -Scaled modified Bessel function of the first kind of order `nu`, ``I_\\nu(x) e^{- | \\operatorname{Re}(x) |}``. -""" -function besselix(nu::Real, x::AbstractFloat) - if x < 0 && !isinteger(nu) - throw(DomainError()) - end - real(besselix(float(nu), complex(x))) -end - -""" - besselj(nu, x) - -Bessel function of the first kind of order `nu`, ``J_\\nu(x)``. -""" -function besselj(nu::Real, x::AbstractFloat) - if isinteger(nu) - if typemin(Cint) <= nu <= typemax(Cint) - return besselj(Cint(nu), x) - end - elseif x < 0 - throw(DomainError()) - end - real(besselj(float(nu), complex(x))) -end - -""" - besseljx(nu, x) - -Scaled Bessel function of the first kind of order `nu`, ``J_\\nu(x) e^{- | \\operatorname{Im}(x) |}``. -""" -function besseljx(nu::Real, x::AbstractFloat) - if x < 0 && !isinteger(nu) - throw(DomainError()) - end - real(besseljx(float(nu), complex(x))) -end - -""" - besselk(nu, x) - -Modified Bessel function of the second kind of order `nu`, ``K_\\nu(x)``. -""" -function besselk(nu::Real, x::AbstractFloat) - if x < 0 - throw(DomainError()) - elseif x == 0 - return oftype(x, Inf) - end - real(besselk(float(nu), complex(x))) -end - -""" - besselkx(nu, x) - -Scaled modified Bessel function of the second kind of order `nu`, ``K_\\nu(x) e^x``. -""" -function besselkx(nu::Real, x::AbstractFloat) - if x < 0 - throw(DomainError()) - elseif x == 0 - return oftype(x, Inf) - end - real(besselkx(float(nu), complex(x))) -end - -""" - bessely(nu, x) - -Bessel function of the second kind of order `nu`, ``Y_\\nu(x)``. -""" -function bessely(nu::Real, x::AbstractFloat) - if x < 0 - throw(DomainError()) - elseif isinteger(nu) && typemin(Cint) <= nu <= typemax(Cint) - return bessely(Cint(nu), x) - end - real(bessely(float(nu), complex(x))) -end - -""" - besselyx(nu, x) - -Scaled Bessel function of the second kind of order `nu`, -``Y_\\nu(x) e^{- | \\operatorname{Im}(x) |}``. -""" -function besselyx(nu::Real, x::AbstractFloat) - if x < 0 - throw(DomainError()) - end - real(besselyx(float(nu), complex(x))) -end - -for f in ("i", "ix", "j", "jx", "k", "kx", "y", "yx") - bfn = Symbol("bessel", f) - @eval begin - $bfn(nu::Real, x::Real) = $bfn(nu, float(x)) - function $bfn(nu::Real, z::Complex) - Tf = promote_type(float(typeof(nu)),float(typeof(real(z)))) - $bfn(Tf(nu), Complex{Tf}(z)) - end - $bfn{T<:AbstractFloat}(k::T, z::Complex{T}) = throw(MethodError($bfn,(k,z))) - $bfn(nu::Float32, x::Complex64) = Complex64($bfn(Float64(nu), Complex128(x))) - end -end - - -for bfn in (:besselh, :besselhx) - @eval begin - $bfn(nu, z) = $bfn(nu, 1, z) - $bfn(nu::Real, k::Integer, x::Real) = $bfn(nu, k, float(x)) - $bfn(nu::Real, k::Integer, x::AbstractFloat) = $bfn(float(nu), k, complex(x)) - - function $bfn(nu::Real, k::Integer, z::Complex) - Tf = promote_type(float(typeof(nu)),float(typeof(real(z)))) - $bfn(Tf(nu), k, Complex{Tf}(z)) - end - - $bfn{T<:AbstractFloat}(nu::T, k::Integer, z::Complex{T}) = throw(MethodError($bfn,(nu,k,z))) - $bfn(nu::Float32, k::Integer, x::Complex64) = Complex64($bfn(Float64(nu), k, Complex128(x))) - end -end - -""" - hankelh1(nu, x) - -Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x)``. -""" -hankelh1(nu, z) = besselh(nu, 1, z) - -""" - hankelh2(nu, x) - -Bessel function of the third kind of order `nu`, ``H^{(2)}_\\nu(x)``. -""" -hankelh2(nu, z) = besselh(nu, 2, z) - -""" - hankelh1x(nu, x) - -Scaled Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x) e^{-x i}``. -""" -hankelh1x(nu, z) = besselhx(nu, 1, z) - -""" - hankelh2x(nu, x) - -Scaled Bessel function of the third kind of order `nu`, ``H^{(2)}_\\nu(x) e^{x i}``. -""" -hankelh2x(nu, z) = besselhx(nu, 2, z) diff --git a/base/special/erf.jl b/base/special/erf.jl deleted file mode 100644 index 2e755fde73b4c8..00000000000000 --- a/base/special/erf.jl +++ /dev/null @@ -1,217 +0,0 @@ -# This file is a part of Julia. License is MIT: http://julialang.org/license - -for f in (:erf, :erfc, :erfcx, :erfi, :Dawson) - fname = (f === :Dawson) ? :dawson : f - @eval begin - ($fname)(z::Complex128) = Complex128(ccall(($(string("Faddeeva_",f)),openspecfun), Complex{Float64}, (Complex{Float64}, Float64), z, zero(Float64))) - ($fname)(z::Complex64) = Complex64(ccall(($(string("Faddeeva_",f)),openspecfun), Complex{Float64}, (Complex{Float64}, Float64), Complex128(z), Float64(eps(Float32)))) - ($fname)(z::Complex) = ($fname)(Complex128(z)) - end -end -for f in (:erfcx, :erfi, :Dawson) - fname = (f === :Dawson) ? :dawson : f - @eval begin - ($fname)(x::Float64) = ccall(($(string("Faddeeva_",f,"_re")),openspecfun), Float64, (Float64,), x) - ($fname)(x::Float32) = Float32(ccall(($(string("Faddeeva_",f,"_re")),openspecfun), Float64, (Float64,), Float64(x))) - ($fname)(x::Integer) = ($fname)(float(x)) - end -end - -# Compute the inverse of the error function: erf(erfinv(x)) == x, -# using the rational approximants tabulated in: -# J. M. Blair, C. A. Edwards, and J. H. Johnson, "Rational Chebyshev -# approximations for the inverse of the error function," Math. Comp. 30, -# pp. 827--830 (1976). -# http://dx.doi.org/10.1090/S0025-5718-1976-0421040-7 -# http://www.jstor.org/stable/2005402 -function erfinv(x::Float64) - a = abs(x) - if a >= 1.0 - if x == 1.0 - return Inf - elseif x == -1.0 - return -Inf - end - throw(DomainError()) - elseif a <= 0.75 # Table 17 in Blair et al. - t = x*x - 0.5625 - return x * @horner(t, 0.16030_49558_44066_229311e2, - -0.90784_95926_29603_26650e2, - 0.18644_91486_16209_87391e3, - -0.16900_14273_46423_82420e3, - 0.65454_66284_79448_7048e2, - -0.86421_30115_87247_794e1, - 0.17605_87821_39059_0) / - @horner(t, 0.14780_64707_15138_316110e2, - -0.91374_16702_42603_13936e2, - 0.21015_79048_62053_17714e3, - -0.22210_25412_18551_32366e3, - 0.10760_45391_60551_23830e3, - -0.20601_07303_28265_443e2, - 0.1e1) - elseif a <= 0.9375 # Table 37 in Blair et al. - t = x*x - 0.87890625 - return x * @horner(t, -0.15238_92634_40726_128e-1, - 0.34445_56924_13612_5216, - -0.29344_39867_25424_78687e1, - 0.11763_50570_52178_27302e2, - -0.22655_29282_31011_04193e2, - 0.19121_33439_65803_30163e2, - -0.54789_27619_59831_8769e1, - 0.23751_66890_24448) / - @horner(t, -0.10846_51696_02059_954e-1, - 0.26106_28885_84307_8511, - -0.24068_31810_43937_57995e1, - 0.10695_12997_33870_14469e2, - -0.23716_71552_15965_81025e2, - 0.24640_15894_39172_84883e2, - -0.10014_37634_97830_70835e2, - 0.1e1) - else # Table 57 in Blair et al. - t = 1.0 / sqrt(-log(1.0 - a)) - return @horner(t, 0.10501_31152_37334_38116e-3, - 0.10532_61131_42333_38164_25e-1, - 0.26987_80273_62432_83544_516, - 0.23268_69578_89196_90806_414e1, - 0.71678_54794_91079_96810_001e1, - 0.85475_61182_21678_27825_185e1, - 0.68738_08807_35438_39802_913e1, - 0.36270_02483_09587_08930_02e1, - 0.88606_27392_96515_46814_9) / - (copysign(t, x) * - @horner(t, 0.10501_26668_70303_37690e-3, - 0.10532_86230_09333_27531_11e-1, - 0.27019_86237_37515_54845_553, - 0.23501_43639_79702_53259_123e1, - 0.76078_02878_58012_77064_351e1, - 0.11181_58610_40569_07827_3451e2, - 0.11948_78791_84353_96667_8438e2, - 0.81922_40974_72699_07893_913e1, - 0.40993_87907_63680_15361_45e1, - 0.1e1)) - end -end - -function erfinv(x::Float32) - a = abs(x) - if a >= 1.0f0 - if x == 1.0f0 - return Inf32 - elseif x == -1.0f0 - return -Inf32 - end - throw(DomainError()) - elseif a <= 0.75f0 # Table 10 in Blair et al. - t = x*x - 0.5625f0 - return x * @horner(t, -0.13095_99674_22f2, - 0.26785_22576_0f2, - -0.92890_57365f1) / - @horner(t, -0.12074_94262_97f2, - 0.30960_61452_9f2, - -0.17149_97799_1f2, - 0.1f1) - elseif a <= 0.9375f0 # Table 29 in Blair et al. - t = x*x - 0.87890625f0 - return x * @horner(t, -0.12402_56522_1f0, - 0.10688_05957_4f1, - -0.19594_55607_8f1, - 0.42305_81357f0) / - @horner(t, -0.88276_97997f-1, - 0.89007_43359f0, - -0.21757_03119_6f1, - 0.1f1) - else # Table 50 in Blair et al. - t = 1.0f0 / sqrt(-log(1.0f0 - a)) - return @horner(t, 0.15504_70003_116f0, - 0.13827_19649_631f1, - 0.69096_93488_87f0, - -0.11280_81391_617f1, - 0.68054_42468_25f0, - -0.16444_15679_1f0) / - (copysign(t, x) * - @horner(t, 0.15502_48498_22f0, - 0.13852_28141_995f1, - 0.1f1)) - end -end - -erfinv(x::Integer) = erfinv(float(x)) - -# Inverse complementary error function: use Blair tables for y = 1-x, -# exploiting the greater accuracy of y (vs. x) when y is small. -function erfcinv(y::Float64) - if y > 0.0625 - return erfinv(1.0 - y) - elseif y <= 0.0 - if y == 0.0 - return Inf - end - throw(DomainError()) - elseif y >= 1e-100 # Table 57 in Blair et al. - t = 1.0 / sqrt(-log(y)) - return @horner(t, 0.10501_31152_37334_38116e-3, - 0.10532_61131_42333_38164_25e-1, - 0.26987_80273_62432_83544_516, - 0.23268_69578_89196_90806_414e1, - 0.71678_54794_91079_96810_001e1, - 0.85475_61182_21678_27825_185e1, - 0.68738_08807_35438_39802_913e1, - 0.36270_02483_09587_08930_02e1, - 0.88606_27392_96515_46814_9) / - (t * - @horner(t, 0.10501_26668_70303_37690e-3, - 0.10532_86230_09333_27531_11e-1, - 0.27019_86237_37515_54845_553, - 0.23501_43639_79702_53259_123e1, - 0.76078_02878_58012_77064_351e1, - 0.11181_58610_40569_07827_3451e2, - 0.11948_78791_84353_96667_8438e2, - 0.81922_40974_72699_07893_913e1, - 0.40993_87907_63680_15361_45e1, - 0.1e1)) - else # Table 80 in Blair et al. - t = 1.0 / sqrt(-log(y)) - return @horner(t, 0.34654_29858_80863_50177e-9, - 0.25084_67920_24075_70520_55e-6, - 0.47378_13196_37286_02986_534e-4, - 0.31312_60375_97786_96408_3388e-2, - 0.77948_76454_41435_36994_854e-1, - 0.70045_68123_35816_43868_271e0, - 0.18710_42034_21679_31668_683e1, - 0.71452_54774_31351_45428_3e0) / - (t * @horner(t, 0.34654_29567_31595_11156e-9, - 0.25084_69079_75880_27114_87e-6, - 0.47379_53129_59749_13536_339e-4, - 0.31320_63536_46177_68848_0813e-2, - 0.78073_48906_27648_97214_733e-1, - 0.70715_04479_95337_58619_993e0, - 0.19998_51543_49112_15105_214e1, - 0.15072_90269_27316_80008_56e1, - 0.1e1)) - end -end - -function erfcinv(y::Float32) - if y > 0.0625f0 - return erfinv(1.0f0 - y) - elseif y <= 0.0f0 - if y == 0.0f0 - return Inf32 - end - throw(DomainError()) - else # Table 50 in Blair et al. - t = 1.0f0 / sqrt(-log(y)) - return @horner(t, 0.15504_70003_116f0, - 0.13827_19649_631f1, - 0.69096_93488_87f0, - -0.11280_81391_617f1, - 0.68054_42468_25f0, - -0.16444_15679_1f0) / - (t * - @horner(t, 0.15502_48498_22f0, - 0.13852_28141_995f1, - 0.1f1)) - end -end - -erfcinv(x::Integer) = erfcinv(float(x)) diff --git a/base/special/gamma.jl b/base/special/gamma.jl index f69e548f03db32..074566530ea6c0 100644 --- a/base/special/gamma.jl +++ b/base/special/gamma.jl @@ -134,347 +134,9 @@ lgamma{T<:Union{Float32,Float16}}(z::Complex{T}) = Complex{T}(lgamma(Complex{Flo gamma(z::Complex) = exp(lgamma(z)) -# Bernoulli numbers B_{2k}, using tabulated numerators and denominators from -# the online encyclopedia of integer sequences. (They actually have data -# up to k=249, but we stop here at k=20.) Used for generating the polygamma -# (Stirling series) coefficients below. -# const A000367 = map(BigInt, split("1,1,-1,1,-1,5,-691,7,-3617,43867,-174611,854513,-236364091,8553103,-23749461029,8615841276005,-7709321041217,2577687858367,-26315271553053477373,2929993913841559,-261082718496449122051", ",")) -# const A002445 = [1,6,30,42,30,66,2730,6,510,798,330,138,2730,6,870,14322,510,6,1919190,6,13530] -# const bernoulli = A000367 .// A002445 # even-index Bernoulli numbers - -""" - digamma(x) - -Compute the digamma function of `x` (the logarithmic derivative of [`gamma(x)`](@ref)). -""" -function digamma(z::Union{Float64,Complex{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 - # argument," Computer Phys. Commun. vol. 4, pp. 221–226 (1972). - x = real(z) - if x <= 0 # reflection formula - ψ = -π * cot(π*z) - z = 1 - z - x = real(z) - else - ψ = zero(z) - end - if x < 7 - # shift using recurrence formula - n = 7 - floor(Int,x) - for ν = 1:n-1 - ψ -= inv(z + ν) - end - ψ -= inv(z) - z += n - end - t = inv(z) - ψ += log(z) - 0.5*t - t *= t # 1/z^2 - # the coefficients here are Float64(bernoulli[2:9] .// (2*(1:8))) - ψ -= t * @evalpoly(t,0.08333333333333333,-0.008333333333333333,0.003968253968253968,-0.004166666666666667,0.007575757575757576,-0.021092796092796094,0.08333333333333333,-0.4432598039215686) -end - -""" - trigamma(x) - -Compute the trigamma function of `x` (the logarithmic second derivative of [`gamma(x)`](@ref)). -""" -function trigamma(z::Union{Float64,Complex{Float64}}) - # via the derivative of the Kölbig digamma formulation - x = real(z) - if x <= 0 # reflection formula - return (π * csc(π*z))^2 - trigamma(1 - z) - end - ψ = zero(z) - if x < 8 - # shift using recurrence formula - n = 8 - floor(Int,x) - ψ += inv(z)^2 - for ν = 1:n-1 - ψ += inv(z + ν)^2 - end - z += n - end - t = inv(z) - w = t * t # 1/z^2 - ψ += t + 0.5*w - # the coefficients here are Float64(bernoulli[2:9]) - ψ += t*w * @evalpoly(w,0.16666666666666666,-0.03333333333333333,0.023809523809523808,-0.03333333333333333,0.07575757575757576,-0.2531135531135531,1.1666666666666667,-7.092156862745098) -end - signflip(m::Number, z) = (-1+0im)^m * z signflip(m::Integer, z) = iseven(m) ? z : -z -# (-1)^m d^m/dz^m cot(z) = p_m(cot z), where p_m is a polynomial -# that satisfies the recurrence p_{m+1}(x) = p_m′(x) * (1 + x^2). -# Note that p_m(x) has only even powers of x if m is odd, and -# only odd powers of x if m is even. Therefore, we write p_m(x) -# as p_m(x) = q_m(x^2) m! for m odd and p_m(x) = x q_m(x^2) m! if m is even. -# Hence the polynomials q_m(y) satisfy the recurrence: -# m odd: q_{m+1}(y) = 2 q_m′(y) * (1+y) / (m+1) -# m even: q_{m+1}(y) = [q_m(y) + 2 y q_m′(y)] * (1+y) / (m+1) -# This function computes the coefficients of the polynomial q_m(y), -# returning an array of the coefficients of 1, y, y^2, ..., -function cotderiv_q(m::Int) - m < 0 && throw(ArgumentError("$m < 0 not allowed")) - m == 0 && return [1.0] - m == 1 && return [1.0, 1.0] - q₋ = cotderiv_q(m-1) - d = length(q₋) - 1 # degree of q₋ - if isodd(m-1) - q = Array{Float64}(length(q₋)) - q[end] = d * q₋[end] * 2/m - for i = 1:length(q)-1 - q[i] = ((i-1)*q₋[i] + i*q₋[i+1]) * 2/m - end - else # iseven(m-1) - q = Array{Float64}(length(q₋) + 1) - q[1] = q₋[1] / m - q[end] = (1 + 2d) * q₋[end] / m - for i = 2:length(q)-1 - q[i] = ((1 + 2(i-1))*q₋[i] + (1 + 2(i-2))*q₋[i-1]) / m - end - end - return q -end - -# precompute a table of cot derivative polynomials -const cotderiv_Q = [cotderiv_q(m) for m in 1:100] - -# Evaluate (-1)^m d^m/dz^m [π cot(πz)] / m!. If m is small, we -# use the explicit derivative formula (a polynomial in cot(πz)); -# if m is large, we use the derivative of Euler's harmonic series: -# π cot(πz) = ∑ inv(z + n) -function cotderiv(m::Integer, z) - isinf(imag(z)) && return zero(z) - if m <= 0 - m == 0 && return π * cot(π*z) - throw(DomainError()) - end - if m <= length(cotderiv_Q) - q = cotderiv_Q[m] - x = cot(π*z) - y = x*x - s = q[1] + q[2] * y - t = y - # evaluate q(y) using Horner (TODO: Knuth for complex z?) - for i = 3:length(q) - t *= y - s += q[i] * t - end - return π^(m+1) * (isodd(m) ? s : x*s) - else # m is large, series derivative should converge quickly - p = m+1 - z -= round(real(z)) - s = inv(z^p) - n = 1 - sₒ = zero(s) - while s != sₒ - sₒ = s - a = (z+n)^p - b = (z-n)^p - s += (a + b) / (a * b) - n += 1 - end - return s - end -end - -# Helper macro for polygamma(m, z): -# Evaluate p[1]*c[1] + x*p[2]*c[2] + x^2*p[3]*c[3] + ... -# where c[1] = m + 1 -# c[k] = c[k-1] * (2k+m-1)*(2k+m-2) / ((2k-1)*(2k-2)) = c[k-1] * d[k] -# i.e. d[k] = c[k]/c[k-1] = (2k+m-1)*(2k+m-2) / ((2k-1)*(2k-2)) -# by a modified version of Horner's rule: -# c[1] * (p[1] + d[2]*x * (p[2] + d[3]*x * (p[3] + ...))). -# The entries of p must be literal constants and there must be > 1 of them. -macro pg_horner(x, m, p...) - k = length(p) - me = esc(m) - xe = esc(x) - ex = :(($me + $(2k-1)) * ($me + $(2k-2)) * $(p[end]/((2k-1)*(2k-2)))) - for k = length(p)-1:-1:2 - cdiv = 1 / ((2k-1)*(2k-2)) - ex = :(($cdiv * ($me + $(2k-1)) * ($me + $(2k-2))) * - ($(p[k]) + $xe * $ex)) - end - :(($me + 1) * ($(p[1]) + $xe * $ex)) -end - -# compute oftype(x, y)^p efficiently, choosing the correct branch cut -pow_oftype(x, y, p) = oftype(x, y)^p -pow_oftype(x::Complex, y::Real, p::Complex) = oftype(x, y^p) -function pow_oftype(x::Complex, y::Real, p::Real) - if p >= 0 - # note: this will never be called for y < 0, - # which would throw an error for non-integer p here - return oftype(x, y^p) - else - yp = y^-p # use real power for efficiency - return oftype(x, Complex(yp, -zero(yp))) # get correct sign of zero! - end -end - -# Generalized zeta function, which is related to polygamma -# (at least for integer m > 0 and real(z) > 0) by: -# polygamma(m, z) = (-1)^(m+1) * gamma(m+1) * zeta(m+1, z). -# Our algorithm for the polygamma is just the m-th derivative -# of our digamma approximation, and this derivative process yields -# a function of the form -# (-1)^(m) * gamma(m+1) * (something) -# So identifying the (something) with the -zeta function, we get -# the zeta function for free and might as well export it, especially -# since this is a common generalization of the Riemann zeta function -# (which Julia already exports). Note that this geneneralization -# is equivalent to Mathematica's Zeta[s,z], and is equivalent to the -# Hurwitz zeta function for real(z) > 0. - -""" - zeta(s, z) - -Generalized zeta function ``\\zeta(s, z)``, defined -by the sum ``\\sum_{k=0}^\\infty ((k+z)^2)^{-s/2}``, where -any term with ``k+z=0`` is excluded. For ``\\Re z > 0``, -this definition is equivalent to the Hurwitz zeta function -``\\sum_{k=0}^\\infty (k+z)^{-s}``. For ``z=1``, it yields -the Riemann zeta function ``\\zeta(s)``. -""" -zeta(s,z) - -function zeta(s::Union{Int,Float64,Complex{Float64}}, - z::Union{Float64,Complex{Float64}}) - ζ = zero(promote_type(typeof(s), typeof(z))) - - (z == 1 || z == 0) && return oftype(ζ, zeta(s)) - s == 2 && return oftype(ζ, trigamma(z)) - - x = real(z) - - # annoying s = Inf or NaN cases: - if !isfinite(s) - (isnan(s) || isnan(z)) && return (s*z)^2 # type-stable NaN+Nan*im - if real(s) == Inf - z == 1 && return one(ζ) - if x > 1 || (x >= 0.5 ? abs(z) > 1 : abs(z - round(x)) > 1) - return zero(ζ) # distance to poles is > 1 - end - x > 0 && imag(z) == 0 && imag(s) == 0 && return oftype(ζ, Inf) - end - throw(DomainError()) # nothing clever to return - end - if isnan(x) - if imag(z)==0 && imag(s)==0 - return oftype(ζ, x) - else - return oftype(ζ, Complex(x,x)) - end - end - - m = s - 1 - - # Algorithm is just the m-th derivative of digamma formula above, - # with a modified cutoff of the final asymptotic expansion. - - # Note: we multiply by -(-1)^m m! in polygamma below, so this factor is - # pulled out of all of our derivatives. - - cutoff = 7 + real(m) + abs(imag(m)) # TODO: this cutoff is too conservative? - if x < cutoff - # shift using recurrence formula - xf = floor(x) - nx = Int(xf) - n = ceil(Int,cutoff - nx) - minus_s = -s - if nx < 0 # x < 0 - # need to use (-z)^(-s) recurrence to be correct for real z < 0 - # [the general form of the recurrence term is (z^2)^(-s/2)] - minus_z = -z - ζ += pow_oftype(ζ, minus_z, minus_s) # ν = 0 term - if xf != z - ζ += pow_oftype(ζ, z - nx, minus_s) # real(z - nx) > 0, so use correct branch cut - # otherwise, if xf==z, then the definition skips this term - end - # do loop in different order, depending on the sign of s, - # so that we are looping from largest to smallest summands and - # can halt the loop early if possible; see issue #15946 - # FIXME: still slow for small m, large Im(z) - if real(s) > 0 - for ν in -nx-1:-1:1 - ζₒ= ζ - ζ += pow_oftype(ζ, minus_z - ν, minus_s) - ζ == ζₒ && break # prevent long loop for large -x > 0 - end - else - for ν in 1:-nx-1 - ζₒ= ζ - ζ += pow_oftype(ζ, minus_z - ν, minus_s) - ζ == ζₒ && break # prevent long loop for large -x > 0 - end - end - else # x ≥ 0 && z != 0 - ζ += pow_oftype(ζ, z, minus_s) - end - # loop order depends on sign of s, as above - if real(s) > 0 - for ν in max(1,1-nx):n-1 - ζₒ= ζ - ζ += pow_oftype(ζ, z + ν, minus_s) - ζ == ζₒ && break # prevent long loop for large m - end - else - for ν in n-1:-1:max(1,1-nx) - ζₒ= ζ - ζ += pow_oftype(ζ, z + ν, minus_s) - ζ == ζₒ && break # prevent long loop for large m - end - end - z += n - end - - t = inv(z) - w = isa(t, Real) ? conj(oftype(ζ, t))^m : oftype(ζ, t)^m - ζ += w * (inv(m) + 0.5*t) - - t *= t # 1/z^2 - ζ += w*t * @pg_horner(t,m,0.08333333333333333,-0.008333333333333333,0.003968253968253968,-0.004166666666666667,0.007575757575757576,-0.021092796092796094,0.08333333333333333,-0.4432598039215686,3.0539543302701198) - - return ζ -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)`](@ref)) -""" -function polygamma(m::Integer, z::Union{Float64,Complex{Float64}}) - m == 0 && return digamma(z) - m == 1 && return trigamma(z) - - # In principle, we could support non-integer m here, but the - # extension to complex m seems to be non-unique, the obvious - # extension (just plugging in a complex m below) does not seem to - # be the one used by Mathematica or Maple, and sources do not - # agree on what the "right" extension is (e.g. Mathematica & Maple - # disagree). So, at least for now, we require integer m. - - # real(m) < 0 is valid, but I don't think our asymptotic expansion - # here works in this case. m < 0 polygamma is called a - # "negapolygamma" function in the literature, and there are - # multiple possible definitions arising from different integration - # constants. We throw a DomainError() since the definition is unclear. - real(m) < 0 && throw(DomainError()) - - s = m+1 - if real(z) <= 0 # reflection formula - (zeta(s, 1-z) + signflip(m, cotderiv(m,z))) * (-gamma(s)) - else - signflip(m, zeta(s,z) * (-gamma(s))) - end -end - # TODO: better way to do this f64(x::Real) = Float64(x) f64(z::Complex) = Complex128(z) @@ -483,62 +145,6 @@ 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. -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 - -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) - @eval begin - $f(z::Number) = $f(f64(z)) - 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 -function invdigamma(y::Float64) - # Closed form initial estimates - if y >= -2.22 - x_old = exp(y) + 0.5 - x_new = x_old - else - x_old = -1.0 / (y - digamma(1.0)) - x_new = x_old - end - - # Fixed point algorithm - delta = Inf - iteration = 0 - while delta > 1e-12 && iteration < 25 - iteration += 1 - x_new = x_old - (digamma(x_old) - y) / trigamma(x_old) - delta = abs(x_new - x_old) - x_old = x_new - end - - return x_new -end -invdigamma(x::Float32) = Float32(invdigamma(Float64(x))) - -""" - invdigamma(x) - -Compute the inverse [`digamma`](@ref) function of `x`. -""" -invdigamma(x::Real) = invdigamma(Float64(x)) - """ beta(x, y) @@ -558,84 +164,3 @@ Natural logarithm of the absolute value of the [`beta`](@ref) 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. -function zeta(s::Union{Float64,Complex{Float64}}) - # blows up to ±Inf, but get correct sign of imaginary zero - s == 1 && return NaN + zero(s) * imag(s) - - if !isfinite(s) # annoying NaN and Inf cases - isnan(s) && return imag(s) == 0 ? s : s*s - if isfinite(imag(s)) - real(s) > 0 && return 1.0 - zero(s)*imag(s) - imag(s) == 0 && return NaN + zero(s) - end - return NaN*zero(s) # NaN + NaN*im - elseif real(s) < 0.5 - if abs(real(s)) + abs(imag(s)) < 1e-3 # Taylor series for small |s| - return @evalpoly(s, -0.5, - -0.918938533204672741780329736405617639861, - -1.0031782279542924256050500133649802190, - -1.00078519447704240796017680222772921424, - -0.9998792995005711649578008136558752359121) - end - return zeta(1 - s) * gamma(1 - s) * sinpi(s*0.5) * (2π)^s / π - end - - m = s - 1 - - # shift using recurrence formula: - # n is a semi-empirical cutoff for the Stirling series, based - # on the error term ~ (|m|/n)^18 / n^real(m) - n = ceil(Int,6 + 0.7*abs(imag(s-1))^inv(1 + real(m)*0.05)) - ζ = one(s) - for ν = 2:n - ζₒ= ζ - ζ += inv(ν)^s - ζ == ζₒ && break # prevent long loop for large m - end - z = 1 + n - t = inv(z) - w = t^m - ζ += w * (inv(m) + 0.5*t) - - t *= t # 1/z^2 - ζ += w*t * @pg_horner(t,m,0.08333333333333333,-0.008333333333333333,0.003968253968253968,-0.004166666666666667,0.007575757575757576,-0.021092796092796094,0.08333333333333333,-0.4432598039215686,3.0539543302701198) - - 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))) - -function eta(z::Union{Float64,Complex{Float64}}) - δz = 1 - z - if abs(real(δz)) + abs(imag(δz)) < 7e-3 # Taylor expand around z==1 - return 0.6931471805599453094172321214581765 * - @evalpoly(δz, - 1.0, - -0.23064207462156020589789602935331414700440, - -0.047156357547388879740146103148112380421254, - -0.002263576552598880778433550956278702759143568, - 0.001081837223249910136105931217561387128141157) - else - 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/doc/src/manual/mathematical-operations.md b/doc/src/manual/mathematical-operations.md index 8e4eea9bbbc384..f25778f895dbe9 100644 --- a/doc/src/manual/mathematical-operations.md +++ b/doc/src/manual/mathematical-operations.md @@ -513,40 +513,8 @@ asind acosd atand acotd asecd acscd | Function | Description | |:------------------------------------------------------------- |:--------------------------------------------------------------------------------------------------------------------------------------------------------------- | -| [`erf(x)`](@ref) | [error function](https://en.wikipedia.org/wiki/Error_function) at `x` | -| [`erfc(x)`](@ref) | complementary error function, i.e. the accurate version of `1-erf(x)` for large `x` | -| [`erfinv(x)`](@ref) | inverse function to [`erf()`](@ref) | -| `erfcinv(x)` | inverse function to [`erfc()`](@ref) | -| [`erfi(x)`](@ref) | imaginary error function defined as `-im * erf(x * im)`, where [`im`](@ref) is the imaginary unit | -| [`erfcx(x)`](@ref) | scaled complementary error function, i.e. accurate `exp(x^2) * erfc(x)` for large `x` | -| [`dawson(x)`](@ref) | scaled imaginary error function, a.k.a. Dawson function, i.e. accurate `exp(-x^2) * erfi(x) * sqrt(pi) / 2` for large `x` | | [`gamma(x)`](@ref) | [gamma function](https://en.wikipedia.org/wiki/Gamma_function) at `x` | | [`lgamma(x)`](@ref) | accurate `log(gamma(x))` for large `x` | | [`lfact(x)`](@ref) | accurate `log(factorial(x))` for large `x`; same as `lgamma(x+1)` for `x > 1`, zero otherwise | -| [`digamma(x)`](@ref) | [digamma function](https://en.wikipedia.org/wiki/Digamma_function) (i.e. the derivative of [`lgamma()`](@ref)) at `x` | | [`beta(x,y)`](@ref) | [beta function](https://en.wikipedia.org/wiki/Beta_function) at `x,y` | | [`lbeta(x,y)`](@ref) | accurate `log(beta(x,y))` for large `x` or `y` | -| [`eta(x)`](@ref) | [Dirichlet eta function](https://en.wikipedia.org/wiki/Dirichlet_eta_function) at `x` | -| [`zeta(x)`](@ref) | [Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function) at `x` | -| [`airyai(z)`](@ref) | [Airy Ai function](https://en.wikipedia.org/wiki/Airy_function) at `z` | -| [`airyaiprime(z)`](@ref) | derivative of the Airy Ai function at `z` | -| [`airybi(z)`](@ref) | [Airy Bi function](https://en.wikipedia.org/wiki/Airy_function) at `z` | -| [`airybiprime(z)`](@ref) | derivative of the Airy Bi function at `z` | -| [`airyaix(z)`](@ref), [`airyaiprimex(z)`](@ref), [`airybix(z)`](@ref), [`airybiprimex(z)`](@ref) | scaled Airy AI function and `k` th derivatives at `z` | -| [`besselj(nu,z)`](@ref) | [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the first kind of order `nu` at `z` | -| [`besselj0(z)`](@ref) | `besselj(0,z)` | -| [`besselj1(z)`](@ref) | `besselj(1,z)` | -| [`besseljx(nu,z)`](@ref) | scaled Bessel function of the first kind of order `nu` at `z` | -| [`bessely(nu,z)`](@ref) | [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the second kind of order `nu` at `z` | -| [`bessely0(z)`](@ref) | `bessely(0,z)` | -| [`bessely1(z)`](@ref) | `bessely(1,z)` | -| [`besselyx(nu,z)`](@ref) | scaled Bessel function of the second kind of order `nu` at `z` | -| [`besselh(nu,k,z)`](@ref) | [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the third kind (a.k.a. Hankel function) of order `nu` at `z`; `k` must be either `1` or `2` | -| [`hankelh1(nu,z)`](@ref) | `besselh(nu, 1, z)` | -| [`hankelh1x(nu,z)`](@ref) | scaled `besselh(nu, 1, z)` | -| [`hankelh2(nu,z)`](@ref) | `besselh(nu, 2, z)` | -| [`hankelh2x(nu,z)`](@ref) | scaled `besselh(nu, 2, z)` | -| [`besseli(nu,z)`](@ref) | modified [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the first kind of order `nu` at `z` | -| [`besselix(nu,z)`](@ref) | scaled modified Bessel function of the first kind of order `nu` at `z` | -| [`besselk(nu,z)`](@ref) | modified [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the second kind of order `nu` at `z` | -| [`besselkx(nu,z)`](@ref) | scaled modified Bessel function of the second kind of order `nu` at `z` | diff --git a/doc/src/stdlib/math.md b/doc/src/stdlib/math.md index dd36ec60d3572f..a195a354e9a0f7 100644 --- a/doc/src/stdlib/math.md +++ b/doc/src/stdlib/math.md @@ -153,13 +153,6 @@ Base.flipsign Base.sqrt Base.isqrt Base.Math.cbrt -Base.Math.erf -Base.Math.erfc -Base.Math.erfcx -Base.Math.erfi -Base.Math.dawson -Base.Math.erfinv -Base.Math.erfcinv Base.real(::Complex) Base.imag Base.reim @@ -182,41 +175,8 @@ Base.powermod Base.Math.gamma Base.Math.lgamma Base.Math.lfact -Base.Math.digamma -Base.Math.invdigamma -Base.Math.trigamma -Base.Math.polygamma -Base.Math.airyai -Base.Math.airyaiprime -Base.Math.airyaix -Base.Math.airyaiprimex -Base.Math.airybi -Base.Math.airybiprime -Base.Math.airybix -Base.Math.airybiprimex -Base.Math.besselj0 -Base.Math.besselj1 -Base.Math.besselj -Base.Math.besseljx -Base.Math.bessely0 -Base.Math.bessely1 -Base.Math.bessely -Base.Math.besselyx -Base.Math.hankelh1 -Base.Math.hankelh1x -Base.Math.hankelh2 -Base.Math.hankelh2x -Base.Math.besselh -Base.Math.besselhx -Base.Math.besseli -Base.Math.besselix -Base.Math.besselk -Base.Math.besselkx Base.Math.beta Base.Math.lbeta -Base.Math.eta -Base.Math.zeta(::Complex) -Base.Math.zeta(::Any, ::Any) Base.ndigits Base.widemul Base.Math.@evalpoly diff --git a/test/bigint.jl b/test/bigint.jl index 9c5b4c75865821..aecb932708f756 100644 --- a/test/bigint.jl +++ b/test/bigint.jl @@ -345,8 +345,6 @@ end @test typeof(exp2(a)) == BigFloat @test typeof(exp10(a)) == BigFloat @test typeof(expm1(a)) == BigFloat -@test typeof(erf(a)) == BigFloat -@test typeof(erfc(a)) == BigFloat @test typeof(cosh(a)) == BigFloat @test typeof(sinh(a)) == BigFloat @test typeof(tanh(a)) == BigFloat diff --git a/test/math.jl b/test/math.jl index b83a2bdba876ce..0e59c03123e410 100644 --- a/test/math.jl +++ b/test/math.jl @@ -382,260 +382,6 @@ end end end -@testset "error functions" begin - @test erf(Float16(1)) ≈ 0.84270079294971486934 - @test erf(1) ≈ 0.84270079294971486934 - @test erfc(1) ≈ 0.15729920705028513066 - @test erfc(Float16(1)) ≈ 0.15729920705028513066 - @test erfcx(1) ≈ 0.42758357615580700442 - @test erfcx(Float32(1)) ≈ 0.42758357615580700442 - @test erfcx(Complex64(1)) ≈ 0.42758357615580700442 - @test erfi(1) ≈ 1.6504257587975428760 - @test erfinv(0.84270079294971486934) ≈ 1 - @test erfcinv(0.15729920705028513066) ≈ 1 - @test dawson(1) ≈ 0.53807950691276841914 - - @test erf(1+2im) ≈ -0.53664356577856503399-5.0491437034470346695im - @test erfc(1+2im) ≈ 1.5366435657785650340+5.0491437034470346695im - @test erfcx(1+2im) ≈ 0.14023958136627794370-0.22221344017989910261im - @test erfi(1+2im) ≈ -0.011259006028815025076+1.0036063427256517509im - @test dawson(1+2im) ≈ -13.388927316482919244-11.828715103889593303im - - for elty in [Float32,Float64] - for x in logspace(-200, -0.01) - @test erf(erfinv(x)) ≈ x atol=1e-12*x - @test erf(erfinv(-x)) ≈ -x atol=1e-12*x - @test erfc(erfcinv(2*x)) ≈ 2*x atol=1e-12*x - if x > 1e-20 - xf = Float32(x) - @test erf(erfinv(xf)) ≈ xf atol=1e-5*xf - @test erf(erfinv(-xf)) ≈ -xf atol=1e-5*xf - @test erfc(erfcinv(2xf)) ≈ 2xf atol=1e-5*xf - end - end - @test erfinv(one(elty)) == Inf - @test erfinv(-one(elty)) == -Inf - @test_throws DomainError erfinv(convert(elty,2.0)) - - @test erfcinv(zero(elty)) == Inf - @test_throws DomainError erfcinv(-one(elty)) - end - - @test erfinv(one(Int)) == erfinv(1.0) - @test erfcinv(one(Int)) == erfcinv(1.0) -end - -@testset "airy" begin - @test_throws Base.Math.AmosException airyai(200im) - @test_throws Base.Math.AmosException airybi(200) - - for T in [Float32, Float64, Complex64,Complex128] - @test airyai(T(1.8)) ≈ 0.0470362168668458052247 - @test airyaiprime(T(1.8)) ≈ -0.0685247801186109345638 - @test airybi(T(1.8)) ≈ 2.595869356743906290060 - @test airybiprime(T(1.8)) ≈ 2.98554005084659907283 - end - for T in [Complex64, Complex128] - z = convert(T,1.8 + 1.0im) - @test airyaix(z) ≈ airyai(z) * exp(2/3 * z * sqrt(z)) - @test airyaiprimex(z) ≈ airyaiprime(z) * exp(2/3 * z * sqrt(z)) - @test airybix(z) ≈ airybi(z) * exp(-abs(real(2/3 * z * sqrt(z)))) - @test airybiprimex(z) ≈ airybiprime(z) * exp(-abs(real(2/3 * z * sqrt(z)))) - end - @test_throws MethodError airyai(complex(big(1.0))) - - for x = -3:3 - @test airyai(x) ≈ airyai(complex(x)) - @test airyaiprime(x) ≈ airyaiprime(complex(x)) - @test airybi(x) ≈ airybi(complex(x)) - @test airybiprime(x) ≈ airybiprime(complex(x)) - if x >= 0 - @test airyaix(x) ≈ airyaix(complex(x)) - @test airyaiprimex(x) ≈ airyaiprimex(complex(x)) - else - @test_throws DomainError airyaix(x) - @test_throws DomainError airyaiprimex(x) - end - @test airybix(x) ≈ airybix(complex(x)) - @test airybiprimex(x) ≈ airybiprimex(complex(x)) - end -end - -@testset "bessel functions" begin - bessel_funcs = [(bessely0, bessely1, bessely), (besselj0, besselj1, besselj)] - @testset "$z, $o" for (z, o, f) in bessel_funcs - @test z(Float32(2.0)) ≈ z(Float64(2.0)) - @test o(Float32(2.0)) ≈ o(Float64(2.0)) - @test z(2) ≈ z(2.0) - @test o(2) ≈ o(2.0) - @test z(2.0 + im) ≈ f(0, 2.0 + im) - @test o(2.0 + im) ≈ f(1, 2.0 + im) - end - @testset "besselj error throwing" begin - @test_throws MethodError besselj(1.2,big(1.0)) - @test_throws MethodError besselj(1,complex(big(1.0))) - @test_throws MethodError besseljx(1,big(1.0)) - @test_throws MethodError besseljx(1,complex(big(1.0))) - end - @testset "besselh" begin - true_h133 = 0.30906272225525164362 - 0.53854161610503161800im - @test besselh(3,1,3) ≈ true_h133 - @test besselh(-3,1,3) ≈ -true_h133 - @test besselh(3,2,3) ≈ conj(true_h133) - @test besselh(-3,2,3) ≈ -conj(true_h133) - @testset "Error throwing" begin - @test_throws Base.Math.AmosException besselh(1,0) - @test_throws MethodError besselh(1,big(1.0)) - @test_throws MethodError besselh(1,complex(big(1.0))) - @test_throws MethodError besselhx(1,big(1.0)) - @test_throws MethodError besselhx(1,complex(big(1.0))) - end - end - @testset "besseli" begin - true_i33 = 0.95975362949600785698 - @test besseli(3,3) ≈ true_i33 - @test besseli(-3,3) ≈ true_i33 - @test besseli(3,-3) ≈ -true_i33 - @test besseli(-3,-3) ≈ -true_i33 - @test besseli(Float32(-3),Complex64(-3,0)) ≈ -true_i33 - @testset "Error throwing" begin - @test_throws Base.Math.AmosException besseli(1,1000) - @test_throws DomainError besseli(0.4,-1.0) - @test_throws MethodError besseli(1,big(1.0)) - @test_throws MethodError besseli(1,complex(big(1.0))) - @test_throws MethodError besselix(1,big(1.0)) - @test_throws MethodError besselix(1,complex(big(1.0))) - end - end - @testset "besselj" begin - @test besselj(0,0) == 1 - for i = 1:5 - @test besselj(i,0) == 0 - @test besselj(-i,0) == 0 - @test besselj(-i,Float32(0)) == 0 - @test besselj(-i,Float32(0)) == 0 - end - - j33 = besselj(3,3.) - @test besselj(3,3) == j33 - @test besselj(-3,-3) == j33 - @test besselj(-3,3) == -j33 - @test besselj(3,-3) == -j33 - @test besselj(3,3f0) ≈ j33 - @test besselj(3,complex(3.)) ≈ j33 - @test besselj(3,complex(3f0)) ≈ j33 - @test besselj(3,complex(3)) ≈ j33 - - j43 = besselj(4,3.) - @test besselj(4,3) == j43 - @test besselj(-4,-3) == j43 - @test besselj(-4,3) == j43 - @test besselj(4,-3) == j43 - @test besselj(4,3f0) ≈ j43 - @test besselj(4,complex(3.)) ≈ j43 - @test besselj(4,complex(3f0)) ≈ j43 - @test besselj(4,complex(3)) ≈ j43 - - @test j33 ≈ 0.30906272225525164362 - @test j43 ≈ 0.13203418392461221033 - @test besselj(0.1, complex(-0.4)) ≈ 0.820421842809028916 + 0.266571215948350899im - @test besselj(3.2, 1.3+0.6im) ≈ 0.01135309305831220201 + 0.03927719044393515275im - @test besselj(1, 3im) ≈ 3.953370217402609396im - @test besselj(1.0,3im) ≈ besselj(1,3im) - @testset "Error throwing" begin - @test_throws DomainError besselj(0.1, -0.4) - @test_throws Base.Math.AmosException besselj(20,1000im) - @test_throws MethodError besselj(big(1.0),3im) - end - end - - @testset "besselk" begin - true_k33 = 0.12217037575718356792 - @test besselk(3,3) ≈ true_k33 - @test besselk(-3,3) ≈ true_k33 - true_k3m3 = -0.1221703757571835679 - 3.0151549516807985776im - @test besselk(3,complex(-3)) ≈ true_k3m3 - @test besselk(-3,complex(-3)) ≈ true_k3m3 - # issue #6564 - @test besselk(1.0,0.0) == Inf - @testset "Error throwing" begin - @test_throws Base.Math.AmosException besselk(200,0.01) - @test_throws DomainError besselk(3,-3) - @test_throws MethodError besselk(1,big(1.0)) - @test_throws MethodError besselk(1,complex(big(1.0))) - @test_throws MethodError besselkx(1,big(1.0)) - @test_throws MethodError besselkx(1,complex(big(1.0))) - end - end - - @testset "bessely" begin - y33 = bessely(3,3.) - @test bessely(3,3) == y33 - @test bessely(3.,3.) == y33 - @test bessely(3,Float32(3.)) ≈ y33 - @test bessely(-3,3) ≈ -y33 - @test y33 ≈ -0.53854161610503161800 - @test bessely(3,complex(-3)) ≈ 0.53854161610503161800 - 0.61812544451050328724im - @testset "Error throwing" begin - @test_throws Base.Math.AmosException bessely(200.5,0.1) - @test_throws DomainError bessely(3,-3) - @test_throws DomainError bessely(0.4,-1.0) - @test_throws DomainError bessely(0.4,Float32(-1.0)) - @test_throws DomainError bessely(1,Float32(-1.0)) - @test_throws DomainError bessely(Cint(3),Float32(-3.)) - @test_throws DomainError bessely(Cint(3),Float64(-3.)) - - @test_throws MethodError bessely(1.2,big(1.0)) - @test_throws MethodError bessely(1,complex(big(1.0))) - @test_throws MethodError besselyx(1,big(1.0)) - @test_throws MethodError besselyx(1,complex(big(1.0))) - end - end - - @testset "besselhx" begin - for elty in [Complex64,Complex128] - z = convert(elty, 1.0 + 1.9im) - @test besselhx(1.0, 1, z) ≈ convert(elty,-0.5949634147786144 - 0.18451272807835967im) - @test besselhx(Float32(1.0), 1, z) ≈ convert(elty,-0.5949634147786144 - 0.18451272807835967im) - end - @testset "Error throwing" begin - @test_throws MethodError besselh(1,1,big(1.0)) - @test_throws MethodError besselh(1,1,complex(big(1.0))) - @test_throws MethodError besselhx(1,1,big(1.0)) - @test_throws MethodError besselhx(1,1,complex(big(1.0))) - end - end - @testset "scaled bessel[ijky] and hankelh[12]" begin - for x in (1.0, 0.0, -1.0), y in (1.0, 0.0, -1.0), nu in (1.0, 0.0, -1.0) - z = Complex128(x + y * im) - z == zero(z) || @test hankelh1x(nu, z) ≈ hankelh1(nu, z) * exp(-z * im) - z == zero(z) || @test hankelh2x(nu, z) ≈ hankelh2(nu, z) * exp(z * im) - (nu < 0 && z == zero(z)) || @test besselix(nu, z) ≈ besseli(nu, z) * exp(-abs(real(z))) - (nu < 0 && z == zero(z)) || @test besseljx(nu, z) ≈ besselj(nu, z) * exp(-abs(imag(z))) - z == zero(z) || @test besselkx(nu, z) ≈ besselk(nu, z) * exp(z) - z == zero(z) || @test besselyx(nu, z) ≈ bessely(nu, z) * exp(-abs(imag(z))) - end - @test besselkx(1, 0) == Inf - @testset "Error throwing" begin - @test_throws Base.Math.AmosException hankelh1x(1, 0) - @test_throws Base.Math.AmosException hankelh2x(1, 0) - @test_throws Base.Math.AmosException besselix(-1, 0) - @test_throws Base.Math.AmosException besseljx(-1, 0) - @test_throws Base.Math.AmosException besselyx(1, 0) - @test_throws DomainError besselix(0.4,-1.0) - @test_throws DomainError besseljx(0.4, -1.0) - @test_throws DomainError besselkx(0.4,-1.0) - @test_throws DomainError besselyx(0.4,-1.0) - end - end - @testset "issue #6653" begin - @testset "$f" for f in (besselj,bessely,besseli,besselk,hankelh1,hankelh2) - @test f(0,1) ≈ f(0,Complex128(1)) - @test f(0,1) ≈ f(0,Complex64(1)) - end - end -end - @testset "beta, lbeta" begin @test beta(3/2,7/2) ≈ 5π/128 @test beta(3,5) ≈ 1/105 @@ -716,81 +462,6 @@ relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x))) @test lgamma(-Inf*im) === -Inf - Inf*im @test lgamma(Inf + Inf*im) === lgamma(NaN + 0im) === lgamma(NaN*im) === NaN + NaN*im end - - @testset "digamma" begin - @testset "$elty" for elty in (Float32, Float64) - @test digamma(convert(elty, 9)) ≈ convert(elty, 2.140641477955609996536345) - @test digamma(convert(elty, 2.5)) ≈ convert(elty, 0.7031566406452431872257) - @test digamma(convert(elty, 0.1)) ≈ convert(elty, -10.42375494041107679516822) - @test digamma(convert(elty, 7e-4)) ≈ convert(elty, -1429.147493371120205005198) - @test digamma(convert(elty, 7e-5)) ≈ convert(elty, -14286.29138623969227538398) - @test digamma(convert(elty, 7e-6)) ≈ convert(elty, -142857.7200612932791081972) - @test digamma(convert(elty, 2e-6)) ≈ convert(elty, -500000.5772123750382073831) - @test digamma(convert(elty, 1e-6)) ≈ convert(elty, -1000000.577214019968668068) - @test digamma(convert(elty, 7e-7)) ≈ convert(elty, -1428572.005785942019703646) - @test digamma(convert(elty, -0.5)) ≈ convert(elty, .03648997397857652055902367) - @test digamma(convert(elty, -1.1)) ≈ convert(elty, 10.15416395914385769902271) - - @test digamma(convert(elty, 0.1)) ≈ convert(elty, -10.42375494041108) - @test digamma(convert(elty, 1/2)) ≈ convert(elty, -γ - log(4)) - @test digamma(convert(elty, 1)) ≈ convert(elty, -γ) - @test digamma(convert(elty, 2)) ≈ convert(elty, 1 - γ) - @test digamma(convert(elty, 3)) ≈ convert(elty, 3/2 - γ) - @test digamma(convert(elty, 4)) ≈ convert(elty, 11/6 - γ) - @test digamma(convert(elty, 5)) ≈ convert(elty, 25/12 - γ) - @test digamma(convert(elty, 10)) ≈ convert(elty, 7129/2520 - γ) - end - end - - @testset "trigamma" begin - @testset "$elty" for elty in (Float32, Float64) - @test trigamma(convert(elty, 0.1)) ≈ convert(elty, 101.433299150792758817) - @test trigamma(convert(elty, 0.1)) ≈ convert(elty, 101.433299150792758817) - @test trigamma(convert(elty, 1/2)) ≈ convert(elty, π^2/2) - @test trigamma(convert(elty, 1)) ≈ convert(elty, π^2/6) - @test trigamma(convert(elty, 2)) ≈ convert(elty, π^2/6 - 1) - @test trigamma(convert(elty, 3)) ≈ convert(elty, π^2/6 - 5/4) - @test trigamma(convert(elty, 4)) ≈ convert(elty, π^2/6 - 49/36) - @test trigamma(convert(elty, 5)) ≈ convert(elty, π^2/6 - 205/144) - @test trigamma(convert(elty, 10)) ≈ convert(elty, π^2/6 - 9778141/6350400) - end - end - - @testset "invdigamma" begin - @testset "$elty" for elty in (Float32, Float64) - for val in [0.001, 0.01, 0.1, 1.0, 10.0] - @test abs(invdigamma(digamma(convert(elty, val))) - convert(elty, val)) < 1e-8 - end - end - @test abs(invdigamma(2)) == abs(invdigamma(2.)) - end - - @testset "polygamma" begin - @test polygamma(20, 7.) ≈ -4.644616027240543262561198814998587152547 - @test polygamma(20, Float16(7.)) ≈ -4.644616027240543262561198814998587152547 - end - - @testset "eta" begin - @test eta(1) ≈ log(2) - @test eta(2) ≈ pi^2/12 - @test eta(Float32(2)) ≈ eta(2) - @test eta(Complex64(2)) ≈ eta(2) - end -end - -@testset "zeta" begin - @test zeta(0) ≈ -0.5 - @test zeta(2) ≈ pi^2/6 - @test zeta(Complex64(2)) ≈ zeta(2) - @test zeta(4) ≈ pi^4/90 - @test zeta(1,Float16(2.)) ≈ zeta(1,2.) - @test zeta(1.,Float16(2.)) ≈ zeta(1,2.) - @test zeta(Float16(1.),Float16(2.)) ≈ zeta(1,2.) - @test isnan(zeta(NaN)) - @test isnan(zeta(1.0e0)) - @test isnan(zeta(1.0f0)) - @test isnan(zeta(complex(0,Inf))) - @test isnan(zeta(complex(-Inf,0))) end @testset "subnormal flags" begin @@ -801,81 +472,6 @@ end @test !get_zero_subnormals() end -#(compared to Wolfram Alpha) -@testset "digamma, trigamma, polygamma & zeta" begin - for x in -10.2:0.3456:50 - @test 1e-12 > relerr(digamma(x+0im), digamma(x)) - end - @test digamma(7+0im) ≅ 1.872784335098467139393487909917597568957840664060076401194232 - @test digamma(7im) ≅ 1.94761433458434866917623737015561385331974500663251349960124 + 1.642224898223468048051567761191050945700191089100087841536im - @test digamma(-3.2+0.1im) ≅ 4.65022505497781398615943030397508454861261537905047116427511+2.32676364843128349629415011622322040021960602904363963042380im - @test trigamma(8+0im) ≅ 0.133137014694031425134546685920401606452509991909746283540546 - @test trigamma(8im) ≅ -0.0078125000000000000029194973110119898029284994355721719150 - 0.12467345030312762782439017882063360876391046513966063947im - @test trigamma(-3.2+0.1im) ≅ 15.2073506449733631753218003030676132587307964766963426965699+15.7081038855113567966903832015076316497656334265029416039199im - @test polygamma(2, 8.1+0im) ≅ -0.01723882695611191078960494454602091934457319791968308929600 - @test polygamma(30, 8.1+2im) ≅ -2722.8895150799704384107961215752996280795801958784600407589+6935.8508929338093162407666304759101854270641674671634631058im - @test polygamma(3, 2.1+1im) ≅ 0.00083328137020421819513475400319288216246978855356531898998-0.27776110819632285785222411186352713789967528250214937861im - @test 1e-11 > relerr(polygamma(3, -4.2 + 2im),-0.0037752884324358856340054736472407163991189965406070325067-0.018937868838708874282432870292420046797798431078848805822im) - @test polygamma(13, 5.2 - 2im) ≅ 0.08087519202975913804697004241042171828113370070289754772448-0.2300264043021038366901951197725318713469156789541415899307im - @test 1e-11 > relerr(polygamma(123, -47.2 + 0im), 5.7111648667225422758966364116222590509254011308116701029e291) - @test zeta(4.1+0.3im, -3.2+0.1im) ≅ -281.34474134962502296077659347175501181994490498591796647 + 286.55601240093672668066037366170168712249413003222992205im - @test zeta(4.1+0.3im, 3.2+0.1im) ≅ 0.0121197525131633219465301571139288562254218365173899270675-0.00687228692565614267981577154948499247518236888933925740902im - @test zeta(4.1, 3.2+0.1im) ≅ 0.0137637451187986846516125754047084829556100290057521276517-0.00152194599531628234517456529686769063828217532350810111482im - @test 1e-12 > relerrc(zeta(1.0001, -4.5e2+3.2im), 10003.765660925877260544923069342257387254716966134385170 - 0.31956240712464746491659767831985629577542514145649468090im) - @test zeta(3.1,-4.2) ≅ zeta(3.1,-4.2+0im) ≅ 149.7591329008219102939352965761913107060718455168339040295 - @test 1e-15 > relerrc(zeta(3.1+0im,-4.2), zeta(3.1,-4.2+0im)) - @test zeta(3.1,4.2) ≅ 0.029938344862645948405021260567725078588893266227472565010234 - @test zeta(27, 3.1) ≅ 5.413318813037879056337862215066960774064332961282599376e-14 - @test zeta(27, 2) ≅ 7.4507117898354294919810041706041194547190318825658299932e-9 - @test 1e-12 > relerr(zeta(27, -105.3), 1.3113726525492708826840989036205762823329453315093955e14) - @test polygamma(4, -3.1+Inf*im) == polygamma(4, 3.1+Inf*im) == 0 - @test polygamma(4, -0.0) == Inf == -polygamma(4, +0.0) - @test zeta(4, +0.0) == zeta(4, -0.0) ≅ pi^4 / 90 - @test zeta(5, +0.0) == zeta(5, -0.0) ≅ 1.036927755143369926331365486457034168057080919501912811974 - @test zeta(Inf, 1.) == 1 - @test zeta(Inf, 2.) == 0 - @test isnan(zeta(NaN, 1.)) - @test isa([digamma(x) for x in [1.0]], Vector{Float64}) - @test isa([trigamma(x) for x in [1.0]], Vector{Float64}) - @test isa([polygamma(3,x) for x in [1.0]], Vector{Float64}) - @test zeta(2 + 1im, -1.1) ≅ zeta(2 + 1im, -1.1+0im) ≅ -64.580137707692178058665068045847533319237536295165484548 + 73.992688148809018073371913557697318846844796582012921247im - @test polygamma(3,5) ≈ polygamma(3,5.) - - @test zeta(-3.0, 7.0) ≅ -52919/120 - @test zeta(-3.0, -7.0) ≅ 94081/120 - @test zeta(-3.1, 7.2) ≅ -587.457736596403704429103489502662574345388906240906317350719 - @test zeta(-3.1, -7.2) ≅ 1042.167459863862249173444363794330893294733001752715542569576 - @test zeta(-3.1, 7.0) ≅ -518.431785723446831868686653718848680989961581500352503093748 - @test zeta(-3.1, -7.0) ≅ 935.1284612957581823462429983411337864448020149908884596048161 - @test zeta(-3.1-0.1im, 7.2) ≅ -579.29752287650299181119859268736614017824853865655709516268 - 96.551907752211554484321948972741033127192063648337407683877im - @test zeta(-3.1-0.1im, -7.2) ≅ 1025.17607931184231774568797674684390615895201417983173984531 + 185.732454778663400767583204948796029540252923367115805842138im - @test zeta(-3.1-0.1im, 7.2 + 0.1im) ≅ -571.66133526455569807299410569274606007165253039948889085762 - 131.86744836357808604785415199791875369679879576524477540653im - @test zeta(-3.1-0.1im, -7.2 + 0.1im) ≅ 1035.35760409421020754141207226034191979220047873089445768189 + 130.905870774271320475424492384335798304480814695778053731045im - @test zeta(-3.1-0.1im, -7.0 + 0.1im) ≅ 929.546530292101383210555114424269079830017210969572819344670 + 113.646687807533854478778193456684618838875194573742062527301im - @test zeta(-3.1, 7.2 + 0.1im) ≅ -586.61801005507638387063781112254388285799318636946559637115 - 36.148831292706044180986261734913443701649622026758378669700im - @test zeta(-3.1, -7.2 + 0.1im) ≅ 1041.04241628770682295952302478199641560576378326778432301623 - 55.7154858634145071137760301929537184886497572057171143541058im - @test zeta(-13.4, 4.1) ≅ -3.860040842156185186414774125656116135638705768861917e6 - @test zeta(3.2, -4) ≅ 2.317164896026427640718298719837102378116771112128525719078 - @test zeta(3.2, 0) ≅ 1.166773370984467020452550350896512026772734054324169010977 - @test zeta(-3.2+0.1im, 0.0) ≅ zeta(-3.2+0.1im, 0.0+0im) ≅ 0.0070547946138977701155565365569392198424378109226519905493 + 0.00076891821792430587745535285452496914239014050562476729610im - @test zeta(-3.2, 0.0) ≅ zeta(-3.2, 0.0+0im) ≅ 0.007011972077091051091698102914884052994997144629191121056378 - - @test 1e-14 > relerr(eta(1+1e-9), 0.693147180719814213126976796937244130533478392539154928250926) - @test 1e-14 > relerr(eta(1+5e-3), 0.693945708117842473436705502427198307157819636785324430166786) - @test 1e-13 > relerr(eta(1+7.1e-3), 0.694280602623782381522315484518617968911346216413679911124758) - @test 1e-13 > relerr(eta(1+8.1e-3), 0.694439974969407464789106040237272613286958025383030083792151) - @test 1e-13 > relerr(eta(1 - 2.1e-3 + 2e-3 * im), 0.69281144248566007063525513903467244218447562492555491581+0.00032001240133205689782368277733081683574922990400416791019im) - @test 1e-13 > relerr(eta(1 + 5e-3 + 5e-3 * im), 0.69394652468453741050544512825906295778565788963009705146+0.00079771059614865948716292388790427833787298296229354721960im) - @test 1e-12 > relerrc(zeta(1e-3+1e-3im), -0.5009189365276307665899456585255302329444338284981610162-0.0009209468912269622649423786878087494828441941303691216750im) - @test 1e-13 > relerrc(zeta(1e-4 + 2e-4im), -0.5000918637469642920007659467492165281457662206388959645-0.0001838278317660822408234942825686513084009527096442173056im) - - # Issue #7169: (TODO: better accuracy should be possible?) - @test 1e-9 > relerrc(zeta(0 + 99.69im), 4.67192766128949471267133846066040655597942700322077493021802+3.89448062985266025394674304029984849370377607524207984092848im) - @test 1e-12 > relerrc(zeta(3 + 99.69im), 1.09996958148566565003471336713642736202442134876588828500-0.00948220959478852115901654819402390826992494044787958181148im) - @test 1e-9 > relerrc(zeta(-3 + 99.69im), 10332.6267578711852982128675093428012860119184786399673520976+13212.8641740351391796168658602382583730208014957452167440726im) - @test 1e-13 > relerrc(zeta(2 + 99.69im, 1.3), 0.41617652544777996034143623540420694985469543821307918291931-0.74199610821536326325073784018327392143031681111201859489991im) -end - @testset "evalpoly" begin @test @evalpoly(2,3,4,5,6) == 3+2*(4+2*(5+2*6)) == @evalpoly(2+0im,3,4,5,6) @test let evalcounts=0 @@ -967,9 +563,7 @@ end @testset "vectorization of 2-arg functions" begin binary_math_functions = [ copysign, flipsign, log, atan2, hypot, max, min, - besselh, hankelh1, hankelh2, hankelh1x, hankelh2x, - besseli, besselix, besselj, besseljx, besselk, besselkx, bessely, besselyx, - polygamma, zeta, beta, lbeta, + beta, lbeta, ] @testset "$f" for f in binary_math_functions x = y = 2 diff --git a/test/mpfr.jl b/test/mpfr.jl index 3496396d99c751..aaacfa55d13b38 100644 --- a/test/mpfr.jl +++ b/test/mpfr.jl @@ -420,16 +420,6 @@ setprecision(256) do @test_throws DomainError factorial(BigFloat(331.3)) end -# bessel functions -setprecision(53) do - @test besselj(4, BigFloat(2)) ≈ besselj(4, 2.) - @test besselj0(BigFloat(2)) ≈ besselj0(2.) - @test besselj1(BigFloat(2)) ≈ besselj1(2.) - @test bessely(4, BigFloat(2)) ≈ bessely(4, 2.) - @test bessely0(BigFloat(2)) ≈ bessely0(2.) - @test bessely1(BigFloat(2)) ≈ bessely1(2.) -end - # trigonometric functions setprecision(53) do for f in (:sin,:cos,:tan,:sec,:csc,:cot,:acos,:asin,:atan, @@ -848,11 +838,6 @@ i3 = trunc(Integer,f) @test i3+1 > f @test i3+1 >= f -let err(z, x) = abs(z - x) / abs(x) - @test 1e-60 > err(eta(parse(BigFloat,"1.005")), parse(BigFloat,"0.693945708117842473436705502427198307157819636785324430166786")) - @test 1e-60 > err(exp(eta(big(1.0))), 2.0) -end - # issue #8318 @test convert(Int64,big(500_000_000_000_000.)) == 500_000_000_000_000 diff --git a/test/random.jl b/test/random.jl index 6d2e9ebcfc52e7..90a126a05d996d 100644 --- a/test/random.jl +++ b/test/random.jl @@ -120,7 +120,8 @@ end ziggurat_table_size = 256 nmantissa = Int64(2)^51 # one bit for the sign ziggurat_nor_r = parse(BigFloat,"3.65415288536100879635194725185604664812733315920964488827246397029393565706474") -nor_section_area = ziggurat_nor_r*exp(-ziggurat_nor_r^2/2) + erfc(ziggurat_nor_r/sqrt(BigFloat(2)))*sqrt(big(π)/2) +erfc_zigg_root2 = parse(BigFloat,"2.580324876539008898343885504487203185398584536409033046076029509351995983934371e-04") +nor_section_area = ziggurat_nor_r*exp(-ziggurat_nor_r^2/2) + erfc_zigg_root2*sqrt(big(π)/2) emantissa = Int64(2)^52 ziggurat_exp_r = parse(BigFloat,"7.69711747013104971404462804811408952334296818528283253278834867283241051210533") exp_section_area = (ziggurat_exp_r + 1)*exp(-ziggurat_exp_r)