diff --git a/base/complex.jl b/base/complex.jl index d074e13bcbd12..3501f0970f35c 100644 --- a/base/complex.jl +++ b/base/complex.jl @@ -366,7 +366,7 @@ function /(z::ComplexF64, w::ComplexF64) cd = max(abs(c), abs(d)) ov = realmax(a) un = realmin(a) - ϵ = eps(Float64) + ϵ = ulp(Float64) bs = two/(ϵ*ϵ) s = 1.0 ab >= half*ov && (a=half*a; b=half*b; s=two*s ) # scale down a,b @@ -399,7 +399,7 @@ function inv(w::ComplexF64) cd = max(abs(c), abs(d)) ov = realmax(c) un = realmin(c) - ϵ = eps(Float64) + ϵ = ulp(Float64) bs = two/(ϵ*ϵ) s = 1.0 cd >= half*ov && (c=half*c; d=half*d; s=s*half) # scale down c,d @@ -424,7 +424,7 @@ function ssqs(x::T, y::T) where T<:AbstractFloat ρ = x*x + y*y if !isfinite(ρ) && (isinf(x) || isinf(y)) ρ = convert(T, Inf) - elseif isinf(ρ) || (ρ==0 && (x!=0 || y!=0)) || ρ 1`, so the only value which can be truncated to + # Here `ulp(Tf(typemin(Ti))) > 1`, so the only value which can be truncated to # `Tf(typemin(Ti)` is itself. Similarly, `Tf(typemax(Ti))` is inexact and will # be rounded up. This assumes that `Tf(typemin(Ti)) > -Inf`, which is true for # these types, but not for `Float16` or larger integer types. @@ -732,11 +732,11 @@ end realmax(::Type{Float32}) = $(bitcast(Float32, 0x7f7fffff)) realmax(::Type{Float64}) = $(bitcast(Float64, 0x7fefffffffffffff)) - eps(x::AbstractFloat) = isfinite(x) ? abs(x) >= realmin(x) ? ldexp(eps(typeof(x)), exponent(x)) : nextfloat(zero(x)) : oftype(x, NaN) - eps(::Type{Float16}) = $(bitcast(Float16, 0x1400)) - eps(::Type{Float32}) = $(bitcast(Float32, 0x34000000)) - eps(::Type{Float64}) = $(bitcast(Float64, 0x3cb0000000000000)) - eps() = eps(Float64) + ulp(x::AbstractFloat) = isfinite(x) ? abs(x) >= realmin(x) ? ldexp(ulp(typeof(x)), exponent(x)) : nextfloat(zero(x)) : oftype(x, NaN) + ulp(::Type{Float16}) = $(bitcast(Float16, 0x1400)) + ulp(::Type{Float32}) = $(bitcast(Float32, 0x34000000)) + ulp(::Type{Float64}) = $(bitcast(Float64, 0x3cb0000000000000)) + ulp() = ulp(Float64) end """ @@ -767,74 +767,74 @@ realmin() = realmin(Float64) realmax() = realmax(Float64) """ - eps(::Type{T}) where T<:AbstractFloat - eps() + ulp(::Type{T}) where T<:AbstractFloat + ulp() Return the *machine epsilon* of the floating point type `T` (`T = Float64` by default). This is defined as the gap between 1 and the next largest value representable by -`typeof(one(T))`, and is equivalent to `eps(one(T))`. (Since `eps(T)` is a +`typeof(one(T))`, and is equivalent to `ulp(one(T))`. (Since `ulp(T)` is a bound on the *relative error* of `T`, it is a "dimensionless" quantity like [`one`](@ref).) # Examples ```jldoctest -julia> eps() +julia> ulp() 2.220446049250313e-16 -julia> eps(Float32) +julia> ulp(Float32) 1.1920929f-7 -julia> 1.0 + eps() +julia> 1.0 + ulp() 1.0000000000000002 -julia> 1.0 + eps()/2 +julia> 1.0 + ulp()/2 1.0 ``` """ -eps(::Type{<:AbstractFloat}) +ulp(::Type{<:AbstractFloat}) """ - eps(x::AbstractFloat) + ulp(x::AbstractFloat) Return the *unit in last place* (ulp) of `x`. This is the distance between consecutive representable floating point values at `x`. In most cases, if the distance on either side of `x` is different, then the larger of the two is taken, that is - eps(x) == max(x-prevfloat(x), nextfloat(x)-x) + ulp(x) == max(x-prevfloat(x), nextfloat(x)-x) The exceptions to this rule are the smallest and largest finite values (e.g. `nextfloat(-Inf)` and `prevfloat(Inf)` for [`Float64`](@ref)), which round to the smaller of the values. -The rationale for this behavior is that `eps` bounds the floating point rounding +The rationale for this behavior is that `ulp` bounds the floating point rounding error. Under the default `RoundNearest` rounding mode, if ``y`` is a real number and ``x`` is the nearest floating point number to ``y``, then ```math -|y-x| \\leq \\operatorname{eps}(x)/2. +|y-x| \\leq \\operatorname{ulp}(x)/2. ``` # Examples ```jldoctest -julia> eps(1.0) +julia> ulp(1.0) 2.220446049250313e-16 -julia> eps(prevfloat(2.0)) +julia> ulp(prevfloat(2.0)) 2.220446049250313e-16 -julia> eps(2.0) +julia> ulp(2.0) 4.440892098500626e-16 julia> x = prevfloat(Inf) # largest finite Float64 1.7976931348623157e308 -julia> x + eps(x)/2 # rounds up +julia> x + ulp(x)/2 # rounds up Inf -julia> x + prevfloat(eps(x)/2) # rounds down +julia> x + prevfloat(ulp(x)/2) # rounds down 1.7976931348623157e308 ``` """ -eps(::AbstractFloat) +ulp(::AbstractFloat) ## byte order swaps for arbitrary-endianness serialization/deserialization ## diff --git a/base/floatfuncs.jl b/base/floatfuncs.jl index 68dad13bcdfa8..e4c9a7e7267d6 100644 --- a/base/floatfuncs.jl +++ b/base/floatfuncs.jl @@ -219,14 +219,14 @@ end # isapprox: approximate equality of numbers """ - isapprox(x, y; rtol::Real=atol>0 ? 0 : √eps, atol::Real=0, nans::Bool=false, norm::Function) + isapprox(x, y; rtol::Real=atol>0 ? 0 : √ulp, atol::Real=0, nans::Bool=false, norm::Function) Inexact equality comparison: `true` if `norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))`. The default `atol` is zero and the default `rtol` depends on the types of `x` and `y`. The keyword argument `nans` determines whether or not NaN values are considered equal (defaults to false). For real or complex floating-point values, if an `atol > 0` is not specified, `rtol` defaults to -the square root of [`eps`](@ref) of the type of `x` or `y`, whichever is bigger (least precise). +the square root of [`ulp`](@ref) of the type of `x` or `y`, whichever is bigger (least precise). This corresponds to requiring equality of about half of the significand digits. Otherwise, e.g. for integer arguments or if an `atol > 0` is supplied, `rtol` defaults to zero. @@ -281,7 +281,7 @@ This is equivalent to `!isapprox(x,y)` (see [`isapprox`](@ref)). ≉(args...; kws...) = !≈(args...; kws...) # default tolerance arguments -rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T)) +rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(ulp(T)) rtoldefault(::Type{<:Real}) = 0 function rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}, atol::Real) where {T<:Number,S<:Number} rtol = max(rtoldefault(real(T)),rtoldefault(real(S))) diff --git a/base/gmp.jl b/base/gmp.jl index 2063fd9328e60..af686fcae7119 100644 --- a/base/gmp.jl +++ b/base/gmp.jl @@ -368,7 +368,7 @@ function (::Type{T})(n::BigInt, ::RoundingMode{:Nearest}) where T<:CdoubleMax x = T(n,RoundToZero) if maxintfloat(T) <= abs(x) < T(Inf) r = n-BigInt(x) - h = eps(x)/2 + h = ulp(x)/2 if iseven(reinterpret(Unsigned,x)) # check if last bit is odd/even if r < -h return prevfloat(x) @@ -623,7 +623,7 @@ function ndigits0zpb(x::BigInt, b::Integer) lb = log2(b) # assumed accurate to <1ulp (true for openlibm) q,r = divrem(n,lb) iq = Int(q) - maxerr = q*eps(lb) # maximum error in remainder + maxerr = q*ulp(lb) # maximum error in remainder if r-1.0 < maxerr abs(x) >= big(b)^iq ? iq+1 : iq elseif lb-r < maxerr diff --git a/base/int.jl b/base/int.jl index e6610e9af4774..6ff2163392ad7 100644 --- a/base/int.jl +++ b/base/int.jl @@ -209,10 +209,10 @@ julia> mod(9, 3) julia> mod(8.9, 3) 2.9000000000000004 -julia> mod(eps(), 3) +julia> mod(ulp(), 3) 2.220446049250313e-16 -julia> mod(-eps(), 3) +julia> mod(-ulp(), 3) 3.0 ``` """ diff --git a/base/irrationals.jl b/base/irrationals.jl index dd47c402f1ce6..8226a7d602aad 100644 --- a/base/irrationals.jl +++ b/base/irrationals.jl @@ -36,7 +36,7 @@ Complex{T}(x::AbstractIrrational) where {T<:Real} = Complex{T}(T(x)) setprecision(BigFloat, p) bx = BigFloat(x) r = rationalize(T, bx, tol=0) - if abs(BigFloat(r) - bx) > eps(bx) + if abs(BigFloat(r) - bx) > ulp(bx) setprecision(BigFloat, o) return r end diff --git a/base/mpfr.jl b/base/mpfr.jl index ef149e7f8bfcf..c4cf8b67f5e60 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -13,7 +13,7 @@ import nextfloat, prevfloat, promote_rule, rem, rem2pi, round, show, float, sum, sqrt, string, print, trunc, precision, exp10, expm1, log1p, - eps, signbit, sin, cos, sincos, tan, sec, csc, cot, acos, asin, atan, + ulp, signbit, sin, cos, sincos, tan, sec, csc, cot, acos, asin, atan, cosh, sinh, tanh, sech, csch, coth, acosh, asinh, atanh, cbrt, typemax, typemin, unsafe_trunc, realmin, realmax, rounding, setrounding, maxintfloat, widen, significand, frexp, tryparse, iszero, @@ -879,7 +879,7 @@ function prevfloat(x::BigFloat) return z end -eps(::Type{BigFloat}) = nextfloat(BigFloat(1)) - BigFloat(1) +ulp(::Type{BigFloat}) = nextfloat(BigFloat(1)) - BigFloat(1) realmin(::Type{BigFloat}) = nextfloat(zero(BigFloat)) realmax(::Type{BigFloat}) = prevfloat(BigFloat(Inf)) diff --git a/base/rational.jl b/base/rational.jl index 6a05b18b7eb3b..25e987261cd21 100644 --- a/base/rational.jl +++ b/base/rational.jl @@ -108,7 +108,7 @@ promote_rule(::Type{Rational{T}}, ::Type{S}) where {T<:Integer,S<:AbstractFloat} widen(::Type{Rational{T}}) where {T} = Rational{widen(T)} """ - rationalize([T<:Integer=Int,] x; tol::Real=eps(x)) + rationalize([T<:Integer=Int,] x; tol::Real=ulp(x)) Approximate floating point number `x` as a [`Rational`](@ref) number with components of the given integer type. The result will differ from `x` by no more than `tol`. @@ -184,7 +184,7 @@ function rationalize(::Type{T}, x::AbstractFloat, tol::Real) where T<:Integer return p // q end end -rationalize(::Type{T}, x::AbstractFloat; tol::Real = eps(x)) where {T<:Integer} = rationalize(T, x, tol)::Rational{T} +rationalize(::Type{T}, x::AbstractFloat; tol::Real = ulp(x)) where {T<:Integer} = rationalize(T, x, tol)::Rational{T} rationalize(x::AbstractFloat; kvs...) = rationalize(Int, x; kvs...) """ diff --git a/base/special/rem_pio2.jl b/base/special/rem_pio2.jl index 67cb067955f7a..75d25d1fcea2b 100644 --- a/base/special/rem_pio2.jl +++ b/base/special/rem_pio2.jl @@ -64,7 +64,7 @@ end fn = round(x*(2/pi)) # round to integer # on older systems, the above could be faster with - # rf = 1.5/eps(Float64) + # rf = 1.5/ulp(Float64) # fn = (x*(2/pi)+rf)-rf r = muladd(-fn, pio2_1, x) # x - fn*pio2_1 diff --git a/base/special/trig.jl b/base/special/trig.jl index 8fa13415c1739..b141307c9251a 100644 --- a/base/special/trig.jl +++ b/base/special/trig.jl @@ -29,7 +29,7 @@ end function sin(x::T) where T<:Union{Float32, Float64} absx = abs(x) if absx < T(pi)/4 #|x| ~<= pi/4, no need for reduction - if absx < sqrt(eps(T)) + if absx < sqrt(ulp(T)) return x end return sin_kernel(x) @@ -99,7 +99,7 @@ end function cos(x::T) where T<:Union{Float32, Float64} absx = abs(x) if absx < T(pi)/4 - if absx < sqrt(eps(T)/T(2.0)) + if absx < sqrt(ulp(T)/T(2.0)) return T(1.0) end return cos_kernel(x) @@ -214,7 +214,7 @@ sincos(x) = _sincos(float(x)) function tan(x::T) where T<:Union{Float32, Float64} absx = abs(x) if absx < T(pi)/4 - if absx < sqrt(eps(T))/2 # first order dominates, but also allows tan(-0)=-0 + if absx < sqrt(ulp(T))/2 # first order dominates, but also allows tan(-0)=-0 return x end return tan_kernel(x) @@ -360,7 +360,7 @@ sincos_kernel(x::Real) = sincos(x) # Inverse trigonometric functions # asin methods ASIN_X_MIN_THRESHOLD(::Type{Float32}) = 2.0f0^-12 -ASIN_X_MIN_THRESHOLD(::Type{Float64}) = sqrt(eps(Float64)) +ASIN_X_MIN_THRESHOLD(::Type{Float64}) = sqrt(ulp(Float64)) arc_p(t::Float64) = t*@horner(t, diff --git a/doc/src/base/base.md b/doc/src/base/base.md index 36e47cc53aae2..e09600c9c1a8f 100644 --- a/doc/src/base/base.md +++ b/doc/src/base/base.md @@ -165,8 +165,8 @@ Base.typemax Base.realmin Base.realmax Base.maxintfloat -Base.eps(::Type{<:AbstractFloat}) -Base.eps(::AbstractFloat) +Base.ulp(::Type{<:AbstractFloat}) +Base.ulp(::AbstractFloat) Base.instances ``` diff --git a/doc/src/manual/integers-and-floating-point-numbers.md b/doc/src/manual/integers-and-floating-point-numbers.md index b218ab8d582e4..db932f2fa5012 100644 --- a/doc/src/manual/integers-and-floating-point-numbers.md +++ b/doc/src/manual/integers-and-floating-point-numbers.md @@ -420,47 +420,50 @@ julia> (typemin(Float64),typemax(Float64)) ### Machine epsilon Most real numbers cannot be represented exactly with floating-point numbers, and so for many purposes -it is important to know the distance between two adjacent representable floating-point numbers, -which is often known as [machine epsilon](https://en.wikipedia.org/wiki/Machine_epsilon). +it is important to know the distance between two adjacent representable floating-point numbers. -Julia provides [`eps`](@ref), which gives the distance between `1.0` and the next larger representable -floating-point value: - -```jldoctest -julia> eps(Float32) -1.1920929f-7 - -julia> eps(Float64) -2.220446049250313e-16 - -julia> eps() # same as eps(Float64) -2.220446049250313e-16 -``` - -These values are `2.0^-23` and `2.0^-52` as [`Float32`](@ref) and [`Float64`](@ref) values, -respectively. The [`eps`](@ref) function can also take a floating-point value as an +The [`ulp`](@ref) function takes a floating-point value as an argument, and gives the absolute difference between that value and the next representable -floating point value. That is, `eps(x)` yields a value of the same type as `x` such that -`x + eps(x)` is the next representable floating-point value larger than `x`: +floating point value. That is, `ulp(x)` yields a value of the same type as `x` such that +`x + ulp(x)` is the next representable floating-point value larger than `x`: ```jldoctest -julia> eps(1.0) +julia> ulp(1.0) 2.220446049250313e-16 -julia> eps(1000.) +julia> ulp(1000.) 1.1368683772161603e-13 -julia> eps(1e-27) +julia> ulp(1e-27) 1.793662034335766e-43 -julia> eps(0.0) +julia> ulp(0.0) 5.0e-324 ``` +For convenience, if the argument value is omitted, it defaults to `1.0`, +and if the argument is a type `T`, it defaults to `one(T)`. For a given floating +point type, this value which is commonly known as the +[machine epsilon](https://en.wikipedia.org/wiki/Machine_epsilon). + +```jldoctest +julia> ulp(Float32) # same as ulp(1.0f0) +1.1920929f-7 + +julia> ulp(Float64) # same as ulp(1.0) +2.220446049250313e-16 + +julia> ulp() # same as ulp(Float64) == ulp(1.0) +2.220446049250313e-16 +``` + +These values are `2.0^-23` and `2.0^-52` as [`Float32`](@ref) and [`Float64`](@ref) values, +respectively. + The distance between two adjacent representable floating-point numbers is not constant, but is smaller for smaller values and larger for larger values. In other words, the representable floating-point numbers are densest in the real number line near zero, and grow sparser exponentially as one moves -farther away from zero. By definition, `eps(1.0)` is the same as `eps(Float64)` since `1.0` is +farther away from zero. By definition, `ulp(1.0)` is the same as `ulp(Float64)` since `1.0` is a 64-bit floating-point value. Julia also provides the [`nextfloat`](@ref) and [`prevfloat`](@ref) functions which return diff --git a/stdlib/Dates/docs/src/index.md b/stdlib/Dates/docs/src/index.md index dc2e7bd7d8a80..d67dbcfd45d71 100644 --- a/stdlib/Dates/docs/src/index.md +++ b/stdlib/Dates/docs/src/index.md @@ -665,7 +665,7 @@ Dates.Time(::Function, ::Any...) Dates.Time(::Dates.DateTime) Dates.now() Dates.now(::Type{Dates.UTC}) -Base.eps +Base.ulp ``` ### Accessor Functions diff --git a/stdlib/Dates/src/types.jl b/stdlib/Dates/src/types.jl index 783e4f9b53b7f..49eaf5bac9276 100644 --- a/stdlib/Dates/src/types.jl +++ b/stdlib/Dates/src/types.jl @@ -321,17 +321,17 @@ calendar(dt::DateTime) = ISOCalendar calendar(dt::Date) = ISOCalendar """ - eps(::DateTime) -> Millisecond - eps(::Date) -> Day - eps(::Time) -> Nanosecond + ulp(::DateTime) -> Millisecond + ulp(::Date) -> Day + ulp(::Time) -> Nanosecond Returns `Millisecond(1)` for `DateTime` values, `Day(1)` for `Date` values, and `Nanosecond(1)` for `Time` values. """ -Base.eps +Base.ulp -Base.eps(dt::DateTime) = Millisecond(1) -Base.eps(dt::Date) = Day(1) -Base.eps(t::Time) = Nanosecond(1) +Base.ulp(dt::DateTime) = Millisecond(1) +Base.ulp(dt::Date) = Day(1) +Base.ulp(t::Time) = Nanosecond(1) Base.typemax(::Union{DateTime, Type{DateTime}}) = DateTime(146138512, 12, 31, 23, 59, 59) Base.typemin(::Union{DateTime, Type{DateTime}}) = DateTime(-146138511, 1, 1, 0, 0, 0) diff --git a/stdlib/Dates/test/types.jl b/stdlib/Dates/test/types.jl index 1cda5b51e6101..f13244faa7793 100644 --- a/stdlib/Dates/test/types.jl +++ b/stdlib/Dates/test/types.jl @@ -181,9 +181,9 @@ c = Dates.Time(0) @testset "DateTime traits" begin @test Dates.calendar(a) == Dates.ISOCalendar @test Dates.calendar(b) == Dates.ISOCalendar - @test eps(a) == Dates.Millisecond(1) - @test eps(b) == Dates.Day(1) - @test eps(c) == Dates.Nanosecond(1) + @test ulp(a) == Dates.Millisecond(1) + @test ulp(b) == Dates.Day(1) + @test ulp(c) == Dates.Nanosecond(1) @test string(typemax(Dates.DateTime)) == "146138512-12-31T23:59:59" @test string(typemin(Dates.DateTime)) == "-146138511-01-01T00:00:00" @test typemax(Dates.DateTime) - typemin(Dates.DateTime) == Dates.Millisecond(9223372017043199000) diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index ee707efe523fb..12b40712dfaa0 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -1248,10 +1248,10 @@ the pseudoinverse by inverting only singular values above a given threshold, The optimal choice of `tol` varies both with the value of `M` and the intended application of the pseudoinverse. The default value of `tol` is -`eps(real(float(one(eltype(M)))))*minimum(size(M))`, which is essentially machine epsilon +`ulp(real(float(one(eltype(M)))))*minimum(size(M))`, which is essentially machine epsilon for the real part of a matrix element multiplied by the larger matrix dimension. For inverting dense ill-conditioned matrices in a least-squares sense, -`tol = sqrt(eps(real(float(one(eltype(M))))))` is recommended. +`tol = sqrt(ulp(real(float(one(eltype(M))))))` is recommended. For more information, see [^issue8859], [^B96], [^S84], [^KY88]. @@ -1311,7 +1311,7 @@ function pinv(A::StridedMatrix{T}, tol::Real) where T return SVD.Vt' * (Diagonal(Sinv) * SVD.U') end function pinv(A::StridedMatrix{T}) where T - tol = eps(real(float(one(T))))*min(size(A)...) + tol = ulp(real(float(one(T))))*min(size(A)...) return pinv(A, tol) end function pinv(x::Number) @@ -1327,7 +1327,7 @@ end Computes a basis for the nullspace of `M` by including the singular vectors of A whose singular have magnitude are greater than `tol*σ₁`, where `σ₁` is `A`'s largest singular values. By default, the value of -`tol` is the smallest dimension of `A` multiplied by the [`eps`](@ref) +`tol` is the smallest dimension of `A` multiplied by the [`ulp`](@ref) of the [`eltype`](@ref) of `A`. # Examples @@ -1351,14 +1351,14 @@ julia> nullspace(M, 2) 0.0 0.0 1.0 ``` """ -function nullspace(A::StridedMatrix, tol::Real = min(size(A)...)*eps(real(float(one(eltype(A)))))) +function nullspace(A::StridedMatrix, tol::Real = min(size(A)...)*ulp(real(float(one(eltype(A)))))) m, n = size(A) (m == 0 || n == 0) && return Matrix{T}(I, n, n) SVD = svd(A, full=true) indstart = sum(SVD.S .> SVD.S[1]*tol) + 1 return copy(SVD.Vt[indstart:end,:]') end -nullspace(a::StridedVector, tol::Real = min(size(a)...)*eps(real(float(one(eltype(a)))))) = nullspace(reshape(a, length(a), 1), tol) +nullspace(a::StridedVector, tol::Real = min(size(a)...)*ulp(real(float(one(eltype(a)))))) = nullspace(reshape(a, length(a), 1), tol) """ cond(M, p::Real=2) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index ac55b9ed696ae..80f96989fda5a 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -720,7 +720,7 @@ end Compute the rank of a matrix by counting how many singular values of `A` have magnitude greater than `tol*σ₁` where `σ₁` is `A`'s largest singular values. By default, the value of `tol` is the smallest -dimension of `A` multiplied by the [`eps`](@ref) +dimension of `A` multiplied by the [`ulp`](@ref) of the [`eltype`](@ref) of `A`. # Examples @@ -738,7 +738,7 @@ julia> rank(diagm(0 => [1, 0.001, 2]), 0.00001) 3 ``` """ -function rank(A::AbstractMatrix, tol::Real = min(size(A)...)*eps(real(float(one(eltype(A)))))) +function rank(A::AbstractMatrix, tol::Real = min(size(A)...)*ulp(real(float(one(eltype(A)))))) s = svdvals(A) count(x -> x > tol*s[1], s) end @@ -1366,7 +1366,7 @@ end rmul!(v, invnrm) else # scale elements to avoid overflow - εδ = eps(one(nrm))/δ + εδ = ulp(one(nrm))/δ rmul!(v, εδ) rmul!(v, inv(nrm*εδ)) end diff --git a/stdlib/LinearAlgebra/src/givens.jl b/stdlib/LinearAlgebra/src/givens.jl index 9113c2e2e1949..becda8050c863 100644 --- a/stdlib/LinearAlgebra/src/givens.jl +++ b/stdlib/LinearAlgebra/src/givens.jl @@ -63,7 +63,7 @@ Base.copy(aR::Adjoint{<:Any,Rotation{T}}) where {T} = Rotation{T}(reverse!([r' f realmin2(::Type{Float32}) = reinterpret(Float32, 0x26000000) realmin2(::Type{Float64}) = reinterpret(Float64, 0x21a0000000000000) -realmin2(::Type{T}) where {T} = (twopar = 2one(T); twopar^trunc(Integer,log(realmin(T)/eps(T))/log(twopar)/twopar)) +realmin2(::Type{T}) where {T} = (twopar = 2one(T); twopar^trunc(Integer,log(realmin(T)/ulp(T))/log(twopar)/twopar)) # derived from LAPACK's dlartg # Copyright: @@ -201,11 +201,11 @@ function givensAlgorithm(f::Complex{T}, g::Complex{T}) where T<:AbstractFloat # g2 is at least safmin, and g2s is at least safmn2 g2s = sqrt(g2) # error in cs from underflow in f2s is at most - # unfl / safmn2 .lt. sqrt(unfl*eps) .lt. eps + # unfl / safmn2 .lt. sqrt(unfl*ulp()) .lt. ulp() # if max(g2,one)=g2, then f2 .lt. g2*safmin, # and so cs .lt. sqrt(safmin) # if max(g2,one)=one, then f2 .lt. safmin - # and so cs .lt. sqrt(safmin)/safmn2 = sqrt(eps) + # and so cs .lt. sqrt(safmin)/safmn2 = sqrt(ulp()) # therefore, cs = f2s/g2s / sqrt( 1 + (f2s/g2s)**2 ) = f2s/g2s cs = f2s/g2s # make sure abs(ff) = 1 diff --git a/stdlib/LinearAlgebra/src/lapack.jl b/stdlib/LinearAlgebra/src/lapack.jl index aa7590f07565b..5cedc321e8baf 100644 --- a/stdlib/LinearAlgebra/src/lapack.jl +++ b/stdlib/LinearAlgebra/src/lapack.jl @@ -1265,7 +1265,7 @@ for (gelsd, gelsy, elty) in # * .. Array Arguments .. # INTEGER JPVT( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) - function gelsy!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=eps($elty)) + function gelsy!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=ulp($elty)) @assert !has_offset_axes(A, B) chkstride1(A) m = size(A, 1) @@ -1363,7 +1363,7 @@ for (gelsd, gelsy, elty, relty) in # INTEGER JPVT( * ) # DOUBLE PRECISION RWORK( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) - function gelsy!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=eps($relty)) + function gelsy!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=ulp($relty)) @assert !has_offset_axes(A, B) chkstride1(A, B) m, n = size(A) diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index 9f32d7c8d530e..a97ac11dd6759 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -775,7 +775,7 @@ end ldiv!(A::QRPivoted{T}, B::StridedVector{T}) where {T<:BlasFloat} = vec(ldiv!(A,reshape(B,length(B),1))) ldiv!(A::QRPivoted{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = - ldiv!(A, B, min(size(A)...)*eps(real(float(one(eltype(B))))))[1] + ldiv!(A, B, min(size(A)...)*ulp(real(float(one(eltype(B))))))[1] function ldiv!(A::QR{T}, B::StridedMatrix{T}) where T m, n = size(A) minmn = min(m,n) diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl index b41c2ab9e4aee..deb868a3bdc75 100644 --- a/stdlib/LinearAlgebra/src/svd.jl +++ b/stdlib/LinearAlgebra/src/svd.jl @@ -220,7 +220,7 @@ svdvals(S::SVD{<:Any,T}) where {T} = (S.S)::Vector{T} # SVD least squares function ldiv!(A::SVD{T}, B::StridedVecOrMat) where T - k = searchsortedlast(A.S, eps(real(T))*A.S[1], rev=true) + k = searchsortedlast(A.S, ulp(real(T))*A.S[1], rev=true) view(A.Vt,1:k,:)' * (view(A.S,1:k) .\ (view(A.U,:,1:k)' * B)) end diff --git a/stdlib/LinearAlgebra/test/bidiag.jl b/stdlib/LinearAlgebra/test/bidiag.jl index f90a13caddaa3..5b50213648366 100644 --- a/stdlib/LinearAlgebra/test/bidiag.jl +++ b/stdlib/LinearAlgebra/test/bidiag.jl @@ -173,15 +173,15 @@ srand(1) if relty != BigFloat x = transpose(T)\transpose(c) tx = transpose(Tfull) \ transpose(c) - elty <: AbstractFloat && @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(promty)*norm(x,Inf)) + elty <: AbstractFloat && @test norm(x-tx,Inf) <= 4*condT*max(ulp()*norm(tx,Inf), ulp(promty)*norm(x,Inf)) @test_throws DimensionMismatch transpose(T)\transpose(b) x = T'\copy(transpose(c)) tx = Tfull'\copy(transpose(c)) - @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(promty)*norm(x,Inf)) + @test norm(x-tx,Inf) <= 4*condT*max(ulp()*norm(tx,Inf), ulp(promty)*norm(x,Inf)) @test_throws DimensionMismatch T'\copy(transpose(b)) x = T\transpose(c) tx = Tfull\transpose(c) - @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(promty)*norm(x,Inf)) + @test norm(x-tx,Inf) <= 4*condT*max(ulp()*norm(tx,Inf), ulp(promty)*norm(x,Inf)) @test_throws DimensionMismatch T\transpose(b) end offsizemat = Matrix{elty}(undef, n+1, 2) @@ -202,7 +202,7 @@ srand(1) x = T \ b tx = Tfull \ b @test_throws DimensionMismatch LinearAlgebra.naivesub!(T,Vector{elty}(undef,n+1)) - @test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(promty)*norm(x,Inf)) + @test norm(x-tx,Inf) <= 4*condT*max(ulp()*norm(tx,Inf), ulp(promty)*norm(x,Inf)) @testset "Generic Mat-vec ops" begin @test T*b ≈ Tfull*b @test T'*b ≈ Tfull'*b @@ -256,7 +256,7 @@ srand(1) test_approx_eq_modphase(u1, u2) test_approx_eq_modphase(copy(v1), copy(v2)) end - @test 0 ≈ norm(u2*Diagonal(d2)*v2'-Tfull) atol=n*max(n^2*eps(relty),norm(u1*Diagonal(d1)*v1'-Tfull)) + @test 0 ≈ norm(u2*Diagonal(d2)*v2'-Tfull) atol=n*max(n^2*ulp(relty),norm(u1*Diagonal(d1)*v1'-Tfull)) @inferred svdvals(T) @inferred svd(T) end diff --git a/stdlib/LinearAlgebra/test/bunchkaufman.jl b/stdlib/LinearAlgebra/test/bunchkaufman.jl index 5b815b10eb4d7..446c882aeb947 100644 --- a/stdlib/LinearAlgebra/test/bunchkaufman.jl +++ b/stdlib/LinearAlgebra/test/bunchkaufman.jl @@ -32,7 +32,7 @@ bimg = randn(n,2)/2 view(a2, 1:n, 1:n), view(aher, 1:n, 1:n), view(apd , 1:n, 1:n))) - ε = εa = eps(abs(float(one(eltya)))) + ε = εa = ulp(abs(float(one(eltya)))) # check that factorize gives a Bunch-Kaufman @test isa(factorize(asym), LinearAlgebra.BunchKaufman) @@ -77,7 +77,7 @@ bimg = randn(n,2)/2 @testset "$eltyb argument B" for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int) b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal) for b in (b, view(b, 1:n, 1:2)) - εb = eps(abs(float(one(eltyb)))) + εb = ulp(abs(float(one(eltyb)))) ε = max(εa,εb) @testset "$uplo Bunch-Kaufman factor of indefinite matrix" for uplo in (:L, :U) @@ -98,7 +98,7 @@ bimg = randn(n,2)/2 @test logabsdet(bc2)[1] ≈ log(abs(det(bc2))) @test logabsdet(bc2)[2] == sign(det(bc2)) @test inv(bc2)*apd ≈ Matrix(I, n, n) - @test apd*(bc2\b) ≈ b rtol=eps(cond(apd)) + @test apd*(bc2\b) ≈ b rtol=ulp(cond(apd)) @test ishermitian(bc2) == !issymmetric(bc2) end end diff --git a/stdlib/LinearAlgebra/test/cholesky.jl b/stdlib/LinearAlgebra/test/cholesky.jl index 4f869dc836ade..413a5c132d085 100644 --- a/stdlib/LinearAlgebra/test/cholesky.jl +++ b/stdlib/LinearAlgebra/test/cholesky.jl @@ -52,7 +52,7 @@ end a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal) a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real) - ε = εa = eps(abs(float(one(eltya)))) + ε = εa = ulp(abs(float(one(eltya)))) # Test of symmetric pos. def. strided matrix apd = a'*a @@ -136,7 +136,7 @@ end for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int) b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal) - εb = eps(abs(float(one(eltyb)))) + εb = ulp(abs(float(one(eltyb)))) ε = max(εa,εb) for b in (b, view(b, 1:n, 1)) # Array and SubArray @@ -215,8 +215,8 @@ end U = Matrix(LinearAlgebra._chol!(X*X', UpperTriangular)[1]) XX = Matrix(X*X') - @test sum(sum(norm, L*L' - XX)) < eps() - @test sum(sum(norm, U'*U - XX)) < eps() + @test sum(sum(norm, L*L' - XX)) < ulp() + @test sum(sum(norm, U'*U - XX)) < ulp() end @@ -260,7 +260,7 @@ end cholesky(Hermitian(apd, :L), Val(true)) \ b r = factorize(apd).U E = abs.(apd - r'*r) - ε = eps(abs(float(one(ComplexF32)))) + ε = ulp(abs(float(one(ComplexF32)))) n = 10 for i=1:n, j=1:n @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) diff --git a/stdlib/LinearAlgebra/test/dense.jl b/stdlib/LinearAlgebra/test/dense.jl index ed192251a958a..58020683f691e 100644 --- a/stdlib/LinearAlgebra/test/dense.jl +++ b/stdlib/LinearAlgebra/test/dense.jl @@ -41,7 +41,7 @@ bimg = randn(n,2)/2 @testset "For A containing $eltya" for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int) ainit = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal) ainit2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real) - ε = εa = eps(abs(float(one(eltya)))) + ε = εa = ulp(abs(float(one(eltya)))) apd = ainit'*ainit # symmetric positive-definite @testset "Positive definiteness" begin @@ -54,7 +54,7 @@ bimg = randn(n,2)/2 end @testset "For b containing $eltyb" for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int) binit = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal) - εb = eps(abs(float(one(eltyb)))) + εb = ulp(abs(float(one(eltyb)))) ε = max(εa,εb) for (a, b) in ((copy(ainit), copy(binit)), (view(ainit, 1:n, 1:n), view(binit, 1:n, 1:2))) @testset "Solve square general system of equations" begin @@ -389,7 +389,7 @@ end end @testset "Additional tests for $elty" for elty in (Float64, Complex{Float64}) - A4 = convert(Matrix{elty}, [1/2 1/3 1/4 1/5+eps(); + A4 = convert(Matrix{elty}, [1/2 1/3 1/4 1/5+ulp(); 1/3 1/4 1/5 1/6; 1/4 1/5 1/6 1/7; 1/5 1/6 1/7 1/8]) @@ -577,7 +577,7 @@ end @test coth(acoth(coth(A))) ≈ coth(A) # Definition of principal values (Aprahamian & Higham, 2016, pp. 4-5) - abstol = sqrt(eps(real(elty))) * norm(acosh(A)) + abstol = sqrt(ulp(real(elty))) * norm(acosh(A)) @test all(z -> (0 < real(z) < π || abs(real(z)) < abstol && imag(z) >= 0 || abs(real(z) - π) < abstol && imag(z) <= 0), @@ -634,7 +634,7 @@ end @test log(A1) ≈ logA1 @test exp(log(A1)) ≈ A1 - A4 = convert(Matrix{elty}, [1/2 1/3 1/4 1/5+eps(); + A4 = convert(Matrix{elty}, [1/2 1/3 1/4 1/5+ulp(); 1/3 1/4 1/5 1/6; 1/4 1/5 1/6 1/7; 1/5 1/6 1/7 1/8]) diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index 028413ef3d318..9e89afcc2d71a 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -75,34 +75,34 @@ srand(1) end for func in (det, tr) - @test func(D) ≈ func(DM) atol=n^2*eps(relty)*(1+(elty<:Complex)) + @test func(D) ≈ func(DM) atol=n^2*ulp(relty)*(1+(elty<:Complex)) end if relty <: BlasFloat for func in (exp, sinh, cosh, tanh, sech, csch, coth) - @test func(D) ≈ func(DM) atol=n^3*eps(relty) + @test func(D) ≈ func(DM) atol=n^3*ulp(relty) end - @test log(Diagonal(abs.(D.diag))) ≈ log(abs.(DM)) atol=n^3*eps(relty) + @test log(Diagonal(abs.(D.diag))) ≈ log(abs.(DM)) atol=n^3*ulp(relty) end if elty <: BlasComplex for func in (logdet, sqrt, sin, cos, tan, sec, csc, cot, asin, acos, atan, asec, acsc, acot, asinh, acosh, atanh, asech, acsch, acoth) - @test func(D) ≈ func(DM) atol=n^2*eps(relty)*2 + @test func(D) ≈ func(DM) atol=n^2*ulp(relty)*2 end end end @testset "Linear solve" begin for (v, U) in ((vv, UU), (view(vv, 1:n), view(UU, 1:n, 1:2))) - @test D*v ≈ DM*v atol=n*eps(relty)*(1+(elty<:Complex)) - @test D*U ≈ DM*U atol=n^2*eps(relty)*(1+(elty<:Complex)) + @test D*v ≈ DM*v atol=n*ulp(relty)*(1+(elty<:Complex)) + @test D*U ≈ DM*U atol=n^2*ulp(relty)*(1+(elty<:Complex)) @test transpose(U)*D ≈ transpose(U)*Array(D) @test U'*D ≈ U'*Array(D) if relty != BigFloat - atol_two = 2n^2 * eps(relty) * (1 + (elty <: Complex)) - atol_three = 2n^3 * eps(relty) * (1 + (elty <: Complex)) + atol_two = 2n^2 * ulp(relty) * (1 + (elty <: Complex)) + atol_three = 2n^3 * ulp(relty) * (1 + (elty <: Complex)) @test D\v ≈ DM\v atol=atol_two @test D\U ≈ DM\U atol=atol_three @test ldiv!(D, copy(v)) ≈ DM\v atol=atol_two @@ -442,7 +442,7 @@ end for (transform1, transform2) in ((identity, identity), (identity, adjoint ), (adjoint, identity ), (adjoint, adjoint ), (identity, transpose), (transpose, identity ), (transpose, transpose) ) - @test *(transform1(D), transform2(B))::typeof(D) ≈ *(transform1(Matrix(D)), transform2(Matrix(B))) atol=2 * eps() + @test *(transform1(D), transform2(B))::typeof(D) ≈ *(transform1(Matrix(D)), transform2(Matrix(B))) atol=2 * ulp() @test *(transform1(DD), transform2(BB))::typeof(DD) == *(transform1(fullDD), transform2(fullBB)) end end diff --git a/stdlib/LinearAlgebra/test/eigen.jl b/stdlib/LinearAlgebra/test/eigen.jl index 7862ae8d78822..3a0a6d04387ed 100644 --- a/stdlib/LinearAlgebra/test/eigen.jl +++ b/stdlib/LinearAlgebra/test/eigen.jl @@ -24,7 +24,7 @@ aimg = randn(n,n)/2 (view(aa, 1:n, 1:n), view(asym, 1:n, 1:n), view(apd, 1:n, 1:n))) - ε = εa = eps(abs(float(one(eltya)))) + ε = εa = ulp(abs(float(one(eltya)))) α = rand(eltya) β = rand(eltya) diff --git a/stdlib/LinearAlgebra/test/generic.jl b/stdlib/LinearAlgebra/test/generic.jl index 306f1d89d4b02..cfcb039e30ee3 100644 --- a/stdlib/LinearAlgebra/test/generic.jl +++ b/stdlib/LinearAlgebra/test/generic.jl @@ -235,9 +235,9 @@ end v = convert(Vector{T}, vr) @test norm(v) == 5.0 w = normalize(v) - @test norm(w - [0.6, 0.8], Inf) < eps(Tr) + @test norm(w - [0.6, 0.8], Inf) < ulp(Tr) @test norm(w) == 1.0 - @test norm(normalize!(copy(v)) - w, Inf) < eps(Tr) + @test norm(normalize!(copy(v)) - w, Inf) < ulp(Tr) @test isempty(normalize!(T[])) end end @@ -251,7 +251,7 @@ end w = normalize(v) @test w ≈ [1/√2, -1/√2] @test norm(w) === 1.0 - @test norm(normalize!(v) - w, Inf) < eps() + @test norm(normalize!(v) - w, Inf) < ulp() end @testset "Issue 14657" begin @@ -269,8 +269,8 @@ end @test LinearAlgebra.promote_leaf_eltypes([[1,2, [3,4]], 5.0, [6im, [7.0, 8.0]]]) == ComplexF64 @test [1, 2, 3] ≈ [1, 2, 3] @test [[1, 2], [3, 4]] ≈ [[1, 2], [3, 4]] - @test [[1, 2], [3, 4]] ≈ [[1.0-eps(), 2.0+eps()], [3.0+2eps(), 4.0-1e8eps()]] - @test [[1, 2], [3, 4]] ≉ [[1.0-eps(), 2.0+eps()], [3.0+2eps(), 4.0-1e9eps()]] + @test [[1, 2], [3, 4]] ≈ [[1.0-ulp(), 2.0+ulp()], [3.0+2ulp(), 4.0-1e8ulp()]] + @test [[1, 2], [3, 4]] ≉ [[1.0-ulp(), 2.0+ulp()], [3.0+2ulp(), 4.0-1e9ulp()]] @test [[1,2, [3,4]], 5.0, [6im, [7.0, 8.0]]] ≈ [[1,2, [3,4]], 5.0, [6im, [7.0, 8.0]]] end diff --git a/stdlib/LinearAlgebra/test/givens.jl b/stdlib/LinearAlgebra/test/givens.jl index e803931cb26d3..dc11e619c49ce 100644 --- a/stdlib/LinearAlgebra/test/givens.jl +++ b/stdlib/LinearAlgebra/test/givens.jl @@ -55,17 +55,17 @@ using LinearAlgebra: rmul!, lmul! end G, r = givens(x[2], x[4], 2, 4) @test (G*x)[2] ≈ r - @test abs((G*x)[4]) < eps(real(elty)) + @test abs((G*x)[4]) < ulp(real(elty)) @inferred givens(x[2], x[4], 2, 4) G, r = givens(x, 2, 4) @test (G*x)[2] ≈ r - @test abs((G*x)[4]) < eps(real(elty)) + @test abs((G*x)[4]) < ulp(real(elty)) @inferred givens(x, 2, 4) G, r = givens(x, 4, 2) @test (G*x)[4] ≈ r - @test abs((G*x)[2]) < eps(real(elty)) + @test abs((G*x)[2]) < ulp(real(elty)) end end end diff --git a/stdlib/LinearAlgebra/test/lapack.jl b/stdlib/LinearAlgebra/test/lapack.jl index 9289a0b3efbc9..26de3c36eff6c 100644 --- a/stdlib/LinearAlgebra/test/lapack.jl +++ b/stdlib/LinearAlgebra/test/lapack.jl @@ -147,7 +147,7 @@ end B = rand(elty, 10, 10) C, j = LAPACK.gelsd!(copy(A),copy(B)) D, k = LAPACK.gelsy!(copy(A),copy(B)) - @test C ≈ D rtol=4*eps(cond(A)) + @test C ≈ D rtol=4*ulp(cond(A)) @test_throws DimensionMismatch LAPACK.gelsd!(A,rand(elty,12,10)) @test_throws DimensionMismatch LAPACK.gelsy!(A,rand(elty,12,10)) end @@ -202,7 +202,7 @@ end Bvs = eigvecs(B) Avs = eigvecs(A) Bvs = LAPACK.gebak!('S','R',ilo,ihi,scale,Bvs) - @test norm(diff(Avs ./ Bvs, dims=1)) < 100 * eps(abs(float(one(elty)))) + @test norm(diff(Avs ./ Bvs, dims=1)) < 100 * ulp(abs(float(one(elty)))) end end @@ -359,7 +359,7 @@ end A = A + transpose(A) #symmetric! B = copy(A) B,ipiv = LAPACK.sytrf!('U',B) - @test triu(inv(A)) ≈ triu(LAPACK.sytri!('U',B,ipiv)) rtol=eps(cond(A)) + @test triu(inv(A)) ≈ triu(LAPACK.sytri!('U',B,ipiv)) rtol=ulp(cond(A)) @test_throws DimensionMismatch LAPACK.sytrs!('U',B,ipiv,rand(elty,11,5)) @test LAPACK.sytrf!('U',zeros(elty,0,0)) == (zeros(elty,0,0),zeros(BlasInt,0),zero(BlasInt)) end @@ -370,7 +370,7 @@ end A = A + transpose(A) #symmetric! B = copy(A) B,ipiv = LAPACK.sytrf_rook!('U', B) - @test triu(inv(A)) ≈ triu(LAPACK.sytri_rook!('U', B, ipiv)) rtol=eps(cond(A)) + @test triu(inv(A)) ≈ triu(LAPACK.sytri_rook!('U', B, ipiv)) rtol=ulp(cond(A)) @test_throws DimensionMismatch LAPACK.sytrs_rook!('U', B, ipiv, rand(elty, 11, 5)) @test LAPACK.sytrf_rook!('U',zeros(elty, 0, 0)) == (zeros(elty, 0, 0),zeros(BlasInt, 0),zero(BlasInt)) A = rand(elty, 10, 10) @@ -379,7 +379,7 @@ end c = A \ b cnd = cond(A) b,A = LAPACK.sysv_rook!('U', A, b) - @test b ≈ c rtol=eps(cnd) + @test b ≈ c rtol=ulp(cnd) @test_throws DimensionMismatch LAPACK.sysv_rook!('U',A,rand(elty,11)) # syconvf_rook error handling @@ -546,7 +546,7 @@ end C = copy(A) D = copy(B) X, rcond, f, b, r = LAPACK.gesvx!(C,D) - @test X ≈ A\B rtol=inv(rcond)*eps(real(elty)) + @test X ≈ A\B rtol=inv(rcond)*ulp(real(elty)) end end diff --git a/stdlib/LinearAlgebra/test/lq.jl b/stdlib/LinearAlgebra/test/lq.jl index d9a7f552de63e..2ce98e220dabe 100644 --- a/stdlib/LinearAlgebra/test/lq.jl +++ b/stdlib/LinearAlgebra/test/lq.jl @@ -29,11 +29,11 @@ rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q) a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real) asym = a' + a # symmetric indefinite apd = a' * a # symmetric positive-definite - ε = εa = eps(abs(float(one(eltya)))) + ε = εa = ulp(abs(float(one(eltya)))) @testset for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int) b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal) - εb = eps(abs(float(one(eltyb)))) + εb = ulp(abs(float(one(eltyb)))) ε = max(εa,εb) α = rand(eltya) diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl index 7cd747c4c9a21..b82bd4499b8a6 100644 --- a/stdlib/LinearAlgebra/test/lu.jl +++ b/stdlib/LinearAlgebra/test/lu.jl @@ -37,7 +37,7 @@ dimg = randn(n)/2 else convert(Tridiagonal{eltya}, Tridiagonal(dlreal, dreal, dureal)) end - ε = εa = eps(abs(float(one(eltya)))) + ε = εa = ulp(abs(float(one(eltya)))) if eltya <: BlasFloat @testset "LU factorization for Number" begin @@ -46,9 +46,9 @@ dimg = randn(n)/2 @test convert(Array, lu(num)) ≈ eltya[num] end @testset "Balancing in eigenvector calculations" begin - A = convert(Matrix{eltya}, [ 3.0 -2.0 -0.9 2*eps(real(one(eltya))); - -2.0 4.0 1.0 -eps(real(one(eltya))); - -eps(real(one(eltya)))/4 eps(real(one(eltya)))/2 -1.0 0; + A = convert(Matrix{eltya}, [ 3.0 -2.0 -0.9 2*ulp(real(one(eltya))); + -2.0 4.0 1.0 -ulp(real(one(eltya))); + -ulp(real(one(eltya)))/4 ulp(real(one(eltya)))/2 -1.0 0; -0.5 -0.5 0.1 1.0]) F = eigen(A, permute=false, scale=false) @test F.vectors*Diagonal(F.values)/F.vectors ≈ A @@ -92,7 +92,7 @@ dimg = randn(n)/2 convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal) c = eltyb == Int ? rand(1:5, n) : convert(Vector{eltyb}, eltyb <: Complex ? complex.(creal, cimg) : creal) - εb = eps(abs(float(one(eltyb)))) + εb = ulp(abs(float(one(eltyb)))) ε = max(εa,εb) @testset "(Automatic) Square LU decomposition" begin lua = factorize(a) diff --git a/stdlib/LinearAlgebra/test/pinv.jl b/stdlib/LinearAlgebra/test/pinv.jl index 489f7c0c8e421..5854299bdb4f7 100644 --- a/stdlib/LinearAlgebra/test/pinv.jl +++ b/stdlib/LinearAlgebra/test/pinv.jl @@ -92,7 +92,7 @@ function test_pinv(a,m,n,tol1,tol2,tol3) @test norm(a*apinv*a-a)/norm(a) ≈ 0 atol=tol1 x0 = randn(n); b = a*x0; x = apinv*b @test norm(a*x-b)/norm(b) ≈ 0 atol=tol1 - apinv = pinv(a,sqrt(eps(real(one(eltype(a)))))) + apinv = pinv(a,sqrt(ulp(real(one(eltype(a)))))) @test norm(a*apinv*a-a)/norm(a) ≈ 0 atol=tol2 x0 = randn(n); b = a*x0; x = apinv*b diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl index 91106676c52f6..ddf933256cb19 100644 --- a/stdlib/LinearAlgebra/test/qr.jl +++ b/stdlib/LinearAlgebra/test/qr.jl @@ -29,11 +29,11 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q) raw_a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real) asym = raw_a' + raw_a # symmetric indefinite apd = raw_a' * raw_a # symmetric positive-definite - ε = εa = eps(abs(float(one(eltya)))) + ε = εa = ulp(abs(float(one(eltya)))) @testset for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int) raw_b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal) - εb = eps(abs(float(one(eltyb)))) + εb = ulp(abs(float(one(eltyb)))) ε = max(εa, εb) tab = promote_type(eltya, eltyb) @@ -174,7 +174,7 @@ end @testset "Issue 7304" begin A = [-√.5 -√.5; -√.5 √.5] Q = rectangularQ(qr(A).Q) - @test norm(A-Q) < eps() + @test norm(A-Q) < ulp() end @testset "qr on AbstractVector" begin @@ -183,7 +183,7 @@ end for T in (Tr, Complex{Tr}) v = convert(Vector{T}, vr) nv, nm = qr(v) - @test norm(nv - [-0.6 -0.8; -0.8 0.6], Inf) < eps(Tr) + @test norm(nv - [-0.6 -0.8; -0.8 0.6], Inf) < ulp(Tr) @test nm == fill(-5.0, 1, 1) end end diff --git a/stdlib/LinearAlgebra/test/schur.jl b/stdlib/LinearAlgebra/test/schur.jl index 488b4b86f15d3..a1b7da8c77fdf 100644 --- a/stdlib/LinearAlgebra/test/schur.jl +++ b/stdlib/LinearAlgebra/test/schur.jl @@ -24,7 +24,7 @@ aimg = randn(n,n)/2 (view(a, 1:n, 1:n), view(asym, 1:n, 1:n), view(apd, 1:n, 1:n))) - ε = εa = eps(abs(float(one(eltya)))) + ε = εa = ulp(abs(float(one(eltya)))) d,v = eigen(a) f = schur(a) diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl index d3377d6412bf1..3a1bd29cc225f 100644 --- a/stdlib/LinearAlgebra/test/svd.jl +++ b/stdlib/LinearAlgebra/test/svd.jl @@ -50,7 +50,7 @@ a2img = randn(n,n)/2 asym = aa' + aa # symmetric indefinite apd = aa' * aa # symmetric positive-definite for (a, a2) in ((aa, aa2), (view(aa, 1:n, 1:n), view(aa2, 1:n, 1:n))) - ε = εa = eps(abs(float(one(eltya)))) + ε = εa = ulp(abs(float(one(eltya)))) usv = svd(a) @testset "singular value decomposition" begin diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 2d4329fab100c..23c13b74f26a6 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -38,7 +38,7 @@ end aherm = a' + a # Hermitian indefinite apos = a' * a # Hermitian positive definite aposs = apos + transpose(apos) # Symmetric positive definite - ε = εa = eps(abs(float(one(eltya)))) + ε = εa = ulp(abs(float(one(eltya)))) x = randn(n) y = randn(n) diff --git a/stdlib/LinearAlgebra/test/testutils.jl b/stdlib/LinearAlgebra/test/testutils.jl index 33eff29765c70..bd21fb3c46ba3 100644 --- a/stdlib/LinearAlgebra/test/testutils.jl +++ b/stdlib/LinearAlgebra/test/testutils.jl @@ -13,12 +13,12 @@ # # Inputs: # a, b:: StridedVecOrMat to be compared -# err :: Default: m^3*(eps(S)+eps(T)), where m is the number of rows +# err :: Default: m^3*(ulp(S)+ulp(T)), where m is the number of rows # # Raises an error if any columnwise vector norm exceeds err. Otherwise, returns # nothing. function test_approx_eq_modphase(a::StridedVecOrMat{S}, b::StridedVecOrMat{T}, - err = length(axes(a,1))^3*(eps(S)+eps(T))) where {S<:Real,T<:Real} + err = length(axes(a,1))^3*(ulp(S)+ulp(T))) where {S<:Real,T<:Real} @test axes(a,1) == axes(b,1) && axes(a,2) == axes(b,2) for i in axes(a,2) v1, v2 = a[:, i], b[:, i] diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl index 13e2aeba68e7a..54688f089f8f6 100644 --- a/stdlib/LinearAlgebra/test/triangular.jl +++ b/stdlib/LinearAlgebra/test/triangular.jl @@ -234,12 +234,12 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo end # Determinant - @test det(A1) ≈ det(lu(Matrix(A1))) atol=sqrt(eps(real(float(one(elty1)))))*n*n - @test logdet(A1) ≈ logdet(lu(Matrix(A1))) atol=sqrt(eps(real(float(one(elty1)))))*n*n + @test det(A1) ≈ det(lu(Matrix(A1))) atol=sqrt(ulp(real(float(one(elty1)))))*n*n + @test logdet(A1) ≈ logdet(lu(Matrix(A1))) atol=sqrt(ulp(real(float(one(elty1)))))*n*n lada, ladb = logabsdet(A1) flada, fladb = logabsdet(lu(Matrix(A1))) - @test lada ≈ flada atol=sqrt(eps(real(float(one(elty1)))))*n*n - @test ladb ≈ fladb atol=sqrt(eps(real(float(one(elty1)))))*n*n + @test lada ≈ flada atol=sqrt(ulp(real(float(one(elty1)))))*n*n + @test ladb ≈ fladb atol=sqrt(ulp(real(float(one(elty1)))))*n*n # Matrix square root @test sqrt(A1) |> t -> t*t ≈ A1 @@ -251,7 +251,7 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo if !(elty1 in (BigFloat, Complex{BigFloat})) # Not handled yet vals, vecs = eigen(A1) if (t1 == UpperTriangular || t1 == LowerTriangular) && elty1 != Int # Cannot really handle degenerate eigen space and Int matrices will probably have repeated eigenvalues. - @test vecs*diagm(0 => vals)/vecs ≈ A1 atol=sqrt(eps(float(real(one(vals[1])))))*(opnorm(A1,Inf)*n)^2 + @test vecs*diagm(0 => vals)/vecs ≈ A1 atol=sqrt(ulp(float(real(one(vals[1])))))*(opnorm(A1,Inf)*n)^2 end end @@ -415,10 +415,10 @@ A2img = randn(n, n)/2 for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int) A = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(Areal, Aimg) : Areal) # a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real) - εa = eps(abs(float(one(eltya)))) + εa = ulp(abs(float(one(eltya)))) for eltyb in (Float32, Float64, ComplexF32, ComplexF64) - εb = eps(abs(float(one(eltyb)))) + εb = ulp(abs(float(one(eltyb)))) ε = max(εa,εb) debug && println("\ntype of A: ", eltya, " type of b: ", eltyb, "\n") diff --git a/stdlib/LinearAlgebra/test/tridiag.jl b/stdlib/LinearAlgebra/test/tridiag.jl index 29e5942ea248f..4dea063c7248a 100644 --- a/stdlib/LinearAlgebra/test/tridiag.jl +++ b/stdlib/LinearAlgebra/test/tridiag.jl @@ -10,7 +10,7 @@ include("testutils.jl") # test_approx_eq_modphase function test_approx_eq_vecs(a::StridedVecOrMat{S}, b::StridedVecOrMat{T}, error=nothing) where {S<:Real,T<:Real} n = size(a, 1) @test n==size(b,1) && size(a,2)==size(b,2) - error===nothing && (error=n^3*(eps(S)+eps(T))) + error===nothing && (error=n^3*(ulp(S)+ulp(T))) for i=1:n ev1, ev2 = a[:,i], b[:,i] deviation = min(abs(norm(ev1-ev2)),abs(norm(ev1+ev2))) @@ -55,7 +55,7 @@ end if elty != Int @testset "issue #1490" begin - @test det(fill(elty(1),3,3)) ≈ zero(elty) atol=3*eps(real(one(elty))) + @test det(fill(elty(1),3,3)) ≈ zero(elty) atol=3*ulp(real(one(elty))) @test det(SymTridiagonal(elty[],elty[])) == one(elty) end end @@ -203,7 +203,7 @@ end if elty != Int @testset "Simple unary functions" begin for func in (det, inv) - @test func(A) ≈ func(fA) atol=n^2*sqrt(eps(real(one(elty)))) + @test func(A) ≈ func(fA) atol=n^2*sqrt(ulp(real(one(elty)))) end end end diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index e73becb43976a..db37fdeb4cd77 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -184,7 +184,7 @@ end for i = 1:5 a = sprand(10, 5, 0.5) b = rand(5) - @test maximum(abs.(a*b - Array(a)*b)) < 100*eps() + @test maximum(abs.(a*b - Array(a)*b)) < 100*ulp() end end @@ -206,86 +206,86 @@ end d = randn(5) + im*randn(5) α = rand(ComplexF64) β = rand(ComplexF64) - @test (maximum(abs.(a*b - Array(a)*b)) < 100*eps()) - @test (maximum(abs.(mul!(similar(b), a, b) - Array(a)*b)) < 100*eps()) # for compatibility with present matmul API. Should go away eventually. - @test (maximum(abs.(mul!(similar(c), a, c) - Array(a)*c)) < 100*eps()) # for compatibility with present matmul API. Should go away eventually. - @test (maximum(abs.(mul!(similar(b), transpose(a), b) - transpose(Array(a))*b)) < 100*eps()) # for compatibility with present matmul API. Should go away eventually. - @test (maximum(abs.(mul!(similar(c), transpose(a), c) - transpose(Array(a))*c)) < 100*eps()) # for compatibility with present matmul API. Should go away eventually. - @test (maximum(abs.(a'b - Array(a)'b)) < 100*eps()) - @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps()) - @test (maximum(abs.(a\b - Array(a)\b)) < 1000*eps()) - @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*eps()) - @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*eps()) - @test (maximum(abs.((a'*c + d) - (Array(a)'*c + d))) < 1000*eps()) - @test (maximum(abs.((α*transpose(a)*c + β*d) - (α*transpose(Array(a))*c + β*d))) < 1000*eps()) - @test (maximum(abs.((transpose(a)*c + d) - (transpose(Array(a))*c + d))) < 1000*eps()) + @test (maximum(abs.(a*b - Array(a)*b)) < 100*ulp()) + @test (maximum(abs.(mul!(similar(b), a, b) - Array(a)*b)) < 100*ulp()) # for compatibility with present matmul API. Should go away eventually. + @test (maximum(abs.(mul!(similar(c), a, c) - Array(a)*c)) < 100*ulp()) # for compatibility with present matmul API. Should go away eventually. + @test (maximum(abs.(mul!(similar(b), transpose(a), b) - transpose(Array(a))*b)) < 100*ulp()) # for compatibility with present matmul API. Should go away eventually. + @test (maximum(abs.(mul!(similar(c), transpose(a), c) - transpose(Array(a))*c)) < 100*ulp()) # for compatibility with present matmul API. Should go away eventually. + @test (maximum(abs.(a'b - Array(a)'b)) < 100*ulp()) + @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*ulp()) + @test (maximum(abs.(a\b - Array(a)\b)) < 1000*ulp()) + @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*ulp()) + @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*ulp()) + @test (maximum(abs.((a'*c + d) - (Array(a)'*c + d))) < 1000*ulp()) + @test (maximum(abs.((α*transpose(a)*c + β*d) - (α*transpose(Array(a))*c + β*d))) < 1000*ulp()) + @test (maximum(abs.((transpose(a)*c + d) - (transpose(Array(a))*c + d))) < 1000*ulp()) c = randn(6) + im*randn(6) @test_throws DimensionMismatch α*transpose(a)*c + β*c @test_throws DimensionMismatch α*transpose(a)*fill(1.,5) + β*c a = I + 0.1*sprandn(5, 5, 0.2) + 0.1*im*sprandn(5, 5, 0.2) b = randn(5,3) - @test (maximum(abs.(a*b - Array(a)*b)) < 100*eps()) - @test (maximum(abs.(a'b - Array(a)'b)) < 100*eps()) - @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps()) - @test (maximum(abs.(a\b - Array(a)\b)) < 1000*eps()) - @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*eps()) - @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*eps()) + @test (maximum(abs.(a*b - Array(a)*b)) < 100*ulp()) + @test (maximum(abs.(a'b - Array(a)'b)) < 100*ulp()) + @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*ulp()) + @test (maximum(abs.(a\b - Array(a)\b)) < 1000*ulp()) + @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*ulp()) + @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*ulp()) a = I + tril(0.1*sprandn(5, 5, 0.2)) b = randn(5,3) + im*randn(5,3) - @test (maximum(abs.(a*b - Array(a)*b)) < 100*eps()) - @test (maximum(abs.(a'b - Array(a)'b)) < 100*eps()) - @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps()) - @test (maximum(abs.(a\b - Array(a)\b)) < 1000*eps()) - @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*eps()) - @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*eps()) + @test (maximum(abs.(a*b - Array(a)*b)) < 100*ulp()) + @test (maximum(abs.(a'b - Array(a)'b)) < 100*ulp()) + @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*ulp()) + @test (maximum(abs.(a\b - Array(a)\b)) < 1000*ulp()) + @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*ulp()) + @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*ulp()) a = I + tril(0.1*sprandn(5, 5, 0.2) + 0.1*im*sprandn(5, 5, 0.2)) b = randn(5,3) - @test (maximum(abs.(a*b - Array(a)*b)) < 100*eps()) - @test (maximum(abs.(a'b - Array(a)'b)) < 100*eps()) - @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps()) - @test (maximum(abs.(a\b - Array(a)\b)) < 1000*eps()) - @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*eps()) - @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*eps()) + @test (maximum(abs.(a*b - Array(a)*b)) < 100*ulp()) + @test (maximum(abs.(a'b - Array(a)'b)) < 100*ulp()) + @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*ulp()) + @test (maximum(abs.(a\b - Array(a)\b)) < 1000*ulp()) + @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*ulp()) + @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*ulp()) a = I + triu(0.1*sprandn(5, 5, 0.2)) b = randn(5,3) + im*randn(5,3) - @test (maximum(abs.(a*b - Array(a)*b)) < 100*eps()) - @test (maximum(abs.(a'b - Array(a)'b)) < 100*eps()) - @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps()) - @test (maximum(abs.(a\b - Array(a)\b)) < 1000*eps()) - @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*eps()) - @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*eps()) + @test (maximum(abs.(a*b - Array(a)*b)) < 100*ulp()) + @test (maximum(abs.(a'b - Array(a)'b)) < 100*ulp()) + @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*ulp()) + @test (maximum(abs.(a\b - Array(a)\b)) < 1000*ulp()) + @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*ulp()) + @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*ulp()) a = I + triu(0.1*sprandn(5, 5, 0.2) + 0.1*im*sprandn(5, 5, 0.2)) b = randn(5,3) - @test (maximum(abs.(a*b - Array(a)*b)) < 100*eps()) - @test (maximum(abs.(a'b - Array(a)'b)) < 100*eps()) - @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps()) - @test (maximum(abs.(a\b - Array(a)\b)) < 1000*eps()) - @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*eps()) - @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*eps()) + @test (maximum(abs.(a*b - Array(a)*b)) < 100*ulp()) + @test (maximum(abs.(a'b - Array(a)'b)) < 100*ulp()) + @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*ulp()) + @test (maximum(abs.(a\b - Array(a)\b)) < 1000*ulp()) + @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*ulp()) + @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*ulp()) a = I + triu(0.1*sprandn(5, 5, 0.2)) b = randn(5,3) + im*randn(5,3) - @test (maximum(abs.(a*b - Array(a)*b)) < 100*eps()) - @test (maximum(abs.(a'b - Array(a)'b)) < 100*eps()) - @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps()) - @test (maximum(abs.(a\b - Array(a)\b)) < 1000*eps()) - @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*eps()) - @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*eps()) + @test (maximum(abs.(a*b - Array(a)*b)) < 100*ulp()) + @test (maximum(abs.(a'b - Array(a)'b)) < 100*ulp()) + @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*ulp()) + @test (maximum(abs.(a\b - Array(a)\b)) < 1000*ulp()) + @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*ulp()) + @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*ulp()) # UpperTriangular/LowerTriangular solve a = UpperTriangular(I + triu(0.1*sprandn(5, 5, 0.2))) b = sprandn(5, 5, 0.2) - @test (maximum(abs.(a\b - Array(a)\Array(b))) < 1000*eps()) + @test (maximum(abs.(a\b - Array(a)\Array(b))) < 1000*ulp()) # test error throwing for bwdTrisolve @test_throws DimensionMismatch a\Matrix{Float64}(I, 6, 6) a = LowerTriangular(I + tril(0.1*sprandn(5, 5, 0.2))) b = sprandn(5, 5, 0.2) - @test (maximum(abs.(a\b - Array(a)\Array(b))) < 1000*eps()) + @test (maximum(abs.(a\b - Array(a)\Array(b))) < 1000*ulp()) # test error throwing for fwdTrisolve @test_throws DimensionMismatch a\Matrix{Float64}(I, 6, 6) @@ -293,20 +293,20 @@ end a = sparse(Diagonal(randn(5) + im*randn(5))) b = randn(5,3) - @test (maximum(abs.(a*b - Array(a)*b)) < 100*eps()) - @test (maximum(abs.(a'b - Array(a)'b)) < 100*eps()) - @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps()) - @test (maximum(abs.(a\b - Array(a)\b)) < 1000*eps()) - @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*eps()) - @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*eps()) + @test (maximum(abs.(a*b - Array(a)*b)) < 100*ulp()) + @test (maximum(abs.(a'b - Array(a)'b)) < 100*ulp()) + @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*ulp()) + @test (maximum(abs.(a\b - Array(a)\b)) < 1000*ulp()) + @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*ulp()) + @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*ulp()) b = randn(5,3) + im*randn(5,3) - @test (maximum(abs.(a*b - Array(a)*b)) < 100*eps()) - @test (maximum(abs.(a'b - Array(a)'b)) < 100*eps()) - @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps()) - @test (maximum(abs.(a\b - Array(a)\b)) < 1000*eps()) - @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*eps()) - @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*eps()) + @test (maximum(abs.(a*b - Array(a)*b)) < 100*ulp()) + @test (maximum(abs.(a'b - Array(a)'b)) < 100*ulp()) + @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*ulp()) + @test (maximum(abs.(a\b - Array(a)\b)) < 1000*ulp()) + @test (maximum(abs.(a'\b - Array(a')\b)) < 1000*ulp()) + @test (maximum(abs.(transpose(a)\b - Array(transpose(a))\b)) < 1000*ulp()) end end end @@ -315,9 +315,9 @@ end for i = 1:5 a = sprand(10, 5, 0.7) b = sprand(5, 15, 0.3) - @test maximum(abs.(a*b - Array(a)*Array(b))) < 100*eps() - @test maximum(abs.(SparseArrays.spmatmul(a,b,sortindices=:sortcols) - Array(a)*Array(b))) < 100*eps() - @test maximum(abs.(SparseArrays.spmatmul(a,b,sortindices=:doubletranspose) - Array(a)*Array(b))) < 100*eps() + @test maximum(abs.(a*b - Array(a)*Array(b))) < 100*ulp() + @test maximum(abs.(SparseArrays.spmatmul(a,b,sortindices=:sortcols) - Array(a)*Array(b))) < 100*ulp() + @test maximum(abs.(SparseArrays.spmatmul(a,b,sortindices=:doubletranspose) - Array(a)*Array(b))) < 100*ulp() f = Diagonal(rand(5)) @test Array(a*f) == Array(a)*f @test Array(f*b) == f*Array(b) @@ -1724,7 +1724,7 @@ end @testset "sparse matrix cond" begin local A = sparse(reshape([1.0], 1, 1)) Ac = sprandn(20, 20,.5) + im*sprandn(20, 20,.5) - Ar = sprandn(20, 20,.5) + eps()*I + Ar = sprandn(20, 20,.5) + ulp()*I @test cond(A, 1) == 1.0 # For a discussion of the tolerance, see #14778 if Base.USE_GPL_LIBS diff --git a/stdlib/SuiteSparse/src/spqr.jl b/stdlib/SuiteSparse/src/spqr.jl index 9aa68e6c5ddea..86524d00b24ce 100644 --- a/stdlib/SuiteSparse/src/spqr.jl +++ b/stdlib/SuiteSparse/src/spqr.jl @@ -137,7 +137,7 @@ Base.axes(Q::QRSparseQ) = map(Base.OneTo, size(Q)) # From SPQR manual p. 6 _default_tol(A::SparseMatrixCSC) = - 20*sum(size(A))*eps(real(eltype(A)))*maximum(norm(view(A, :, i)) for i in 1:size(A, 2)) + 20*sum(size(A))*ulp(real(eltype(A)))*maximum(norm(view(A, :, i)) for i in 1:size(A, 2)) function LinearAlgebra.qr(A::SparseMatrixCSC{Tv}; tol = _default_tol(A)) where {Tv <: CHOLMOD.VTypes} R = Ref{Ptr{CHOLMOD.C_Sparse{Tv}}}() diff --git a/stdlib/Test/src/Test.jl b/stdlib/Test/src/Test.jl index 9bba070655481..0c2be1f4a9c16 100644 --- a/stdlib/Test/src/Test.jl +++ b/stdlib/Test/src/Test.jl @@ -1614,8 +1614,8 @@ begin end end - array_eps(a::AbstractArray{Complex{T}}) where {T} = eps(float(maximum(x->(isfinite(x) ? abs(x) : T(NaN)), a))) - array_eps(a) = eps(float(maximum(x->(isfinite(x) ? abs(x) : oftype(x,NaN)), a))) + array_eps(a::AbstractArray{Complex{T}}) where {T} = ulp(float(maximum(x->(isfinite(x) ? abs(x) : T(NaN)), a))) + array_eps(a) = ulp(float(maximum(x->(isfinite(x) ? abs(x) : oftype(x,NaN)), a))) test_approx_eq(va, vb, astr, bstr) = test_approx_eq(va, vb, 1E4*length(LinearIndices(va))*max(array_eps(va), array_eps(vb)), astr, bstr) diff --git a/test/arrayops.jl b/test/arrayops.jl index e5ace0234d2a9..b9473829b60bd 100644 --- a/test/arrayops.jl +++ b/test/arrayops.jl @@ -1796,7 +1796,7 @@ end # issue #9648 let x = fill(1.5f0, 10^7) - @test abs(1.5f7 - cumsum(x)[end]) < 3*eps(1.5f7) + @test abs(1.5f7 - cumsum(x)[end]) < 3*ulp(1.5f7) @test cumsum(x) == cumsum!(similar(x), x) end @@ -2343,7 +2343,7 @@ Float64(x::F21666) = Float64(x.x) # test that cumsum uses more stable algorithm # for types with unknown/rounding arithmetic # we make v pretty large, because stable algorithm may have a large base case - v = zeros(300); v[1] = 2; v[200:end] .= eps(Float32) + v = zeros(300); v[1] = 2; v[200:end] .= ulp(Float32) f_rounds = Float64.(cumsum(F21666{Base.ArithmeticRounds}.(v))) f_unknown = Float64.(cumsum(F21666{Base.ArithmeticUnknown}.(v))) diff --git a/test/complex.jl b/test/complex.jl index 67be8f2a5fa4d..95d565c0a4f3d 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -103,7 +103,7 @@ end @test exp(x) ≈ exp(big(x)) @test exp10(x) ≈ exp10(big(x)) @test exp2(x) ≈ exp2(big(x)) - @test expm1(x) ≈ expm1(big(x)) atol=eps(T) + @test expm1(x) ≈ expm1(big(x)) atol=ulp(T) @test log(x) ≈ log(big(x)) @test log10(x) ≈ log10(big(x)) @test log1p(x) ≈ log1p(big(x)) diff --git a/test/fastmath.jl b/test/fastmath.jl index 8c82df9adf8bb..294083d3af8b6 100644 --- a/test/fastmath.jl +++ b/test/fastmath.jl @@ -11,13 +11,13 @@ @test macroexpand(Main, :(@fastmath sincos(x))) == :(Base.FastMath.sincos_fast(x)) end const one32 = one(Float32) -const eps32 = eps(Float32) +const eps32 = ulp(Float32) const eps32_2 = eps32/2 # Note: Cannot use local functions since these are not yet optimized fm_ieee_32(x) = x + eps32_2 + eps32_2 fm_fast_32(x) = @fastmath x + eps32_2 + eps32_2 const one64 = one(Float64) -const eps64 = eps(Float64) +const eps64 = ulp(Float64) const eps64_2 = eps64/2 # Note: Cannot use local functions since these are not yet optimized fm_ieee_64(x) = x + eps64_2 + eps64_2 @@ -45,7 +45,7 @@ fm_fast_64_upd(x) = @fastmath (r=x; r+=eps64_2; r+=eps64_2) for T in (Float32, Float64, BigFloat) zero = convert(T, 0) - one = convert(T, 1) + eps(T) + one = convert(T, 1) + ulp(T) two = convert(T, 2) + 1//10 three = convert(T, 3) + 1//100 @@ -74,7 +74,7 @@ fm_fast_64_upd(x) = @fastmath (r=x; r+=eps64_2; r+=eps64_2) for T in (ComplexF32, ComplexF64, Complex{BigFloat}) zero = convert(T,0) - one = convert(T,1) + im*eps(real(convert(T,1))) + one = convert(T,1) + im*ulp(real(convert(T,1))) two = convert(T,2) + im//10 three = convert(T,3) + im//100 diff --git a/test/math.jl b/test/math.jl index f4dc7f874a736..e059aa7016378 100644 --- a/test/math.jl +++ b/test/math.jl @@ -173,38 +173,38 @@ end @test isequal(T(1//4)^2, T(1//16)) @test isequal(acos(T(1)), T(0)) @test isequal(acosh(T(1)), T(0)) - @test asin(T(1)) ≈ T(pi)/2 atol=eps(T) - @test atan(T(1)) ≈ T(pi)/4 atol=eps(T) - @test atan(T(1),T(1)) ≈ T(pi)/4 atol=eps(T) + @test asin(T(1)) ≈ T(pi)/2 atol=ulp(T) + @test atan(T(1)) ≈ T(pi)/4 atol=ulp(T) + @test atan(T(1),T(1)) ≈ T(pi)/4 atol=ulp(T) @test isequal(cbrt(T(0)), T(0)) @test isequal(cbrt(T(1)), T(1)) @test isequal(cbrt(T(1000000000)), T(1000)) @test isequal(cos(T(0)), T(1)) - @test cos(T(pi)/2) ≈ T(0) atol=eps(T) + @test cos(T(pi)/2) ≈ T(0) atol=ulp(T) @test isequal(cos(T(pi)), T(-1)) - @test exp(T(1)) ≈ T(ℯ) atol=10*eps(T) + @test exp(T(1)) ≈ T(ℯ) atol=10*ulp(T) @test isequal(exp10(T(1)), T(10)) @test isequal(exp2(T(1)), T(2)) @test isequal(expm1(T(0)), T(0)) - @test expm1(T(1)) ≈ T(ℯ)-1 atol=10*eps(T) + @test expm1(T(1)) ≈ T(ℯ)-1 atol=10*ulp(T) @test isequal(hypot(T(3),T(4)), T(5)) @test isequal(log(T(1)), T(0)) @test isequal(log(ℯ,T(1)), T(0)) - @test log(T(ℯ)) ≈ T(1) atol=eps(T) + @test log(T(ℯ)) ≈ T(1) atol=ulp(T) @test isequal(log10(T(1)), T(0)) @test isequal(log10(T(10)), T(1)) @test isequal(log1p(T(0)), T(0)) - @test log1p(T(ℯ)-1) ≈ T(1) atol=eps(T) + @test log1p(T(ℯ)-1) ≈ T(1) atol=ulp(T) @test isequal(log2(T(1)), T(0)) @test isequal(log2(T(2)), T(1)) @test isequal(sin(T(0)), T(0)) @test isequal(sin(T(pi)/2), T(1)) - @test sin(T(pi)) ≈ T(0) atol=eps(T) + @test sin(T(pi)) ≈ T(0) atol=ulp(T) @test isequal(sqrt(T(0)), T(0)) @test isequal(sqrt(T(1)), T(1)) @test isequal(sqrt(T(100000000)), T(10000)) @test isequal(tan(T(0)), T(0)) - @test tan(T(pi)/4) ≈ T(1) atol=eps(T) + @test tan(T(pi)/4) ≈ T(1) atol=ulp(T) end @testset "Inverses" begin @test acos(cos(x)) ≈ x @@ -263,7 +263,7 @@ end X = map(T, vcat(-10:0.0002:10, -80:0.001:80, 2.0^-27, 2.0^-28, 2.0^-14, 2.0^-13)) for x in X y, yb = exp(x), exp(big(x)) - @test abs(y-yb) <= 1.0*eps(T(yb)) + @test abs(y-yb) <= 1.0*ulp(T(yb)) end end @testset "$T edge cases" begin @@ -281,12 +281,12 @@ end X = map(Float64, vcat(-10:0.00021:10, -35:0.0023:100, -300:0.001:300)) for x in X y, yb = exp10(x), exp10(big(x)) - @test abs(y-yb) <= 1.2*eps(Float64(yb)) + @test abs(y-yb) <= 1.2*ulp(Float64(yb)) end X = map(Float32, vcat(-10:0.00021:10, -35:0.0023:35, -35:0.001:35)) for x in X y, yb = exp10(x), exp10(big(x)) - @test abs(y-yb) <= 1.2*eps(Float32(yb)) + @test abs(y-yb) <= 1.2*ulp(Float32(yb)) end end @testset "$T edge cases" for T in (Float64, Float32) @@ -346,8 +346,8 @@ end @testset "$T" for T = (Float32,Float64,Rational{Int}) fT = typeof(float(one(T))) for x = -400:40:400 - @test sind(convert(T,x))::fT ≈ convert(fT,sin(pi/180*x)) atol=eps(deg2rad(convert(fT,x))) - @test cosd(convert(T,x))::fT ≈ convert(fT,cos(pi/180*x)) atol=eps(deg2rad(convert(fT,x))) + @test sind(convert(T,x))::fT ≈ convert(fT,sin(pi/180*x)) atol=ulp(deg2rad(convert(fT,x))) + @test cosd(convert(T,x))::fT ≈ convert(fT,cos(pi/180*x)) atol=ulp(deg2rad(convert(fT,x))) end @testset "sind" begin @test sind(convert(T,0.0))::fT === zero(fT) @@ -366,8 +366,8 @@ end @testset "sinpi and cospi" begin for x = -3:0.3:3 - @test sinpi(convert(T,x))::fT ≈ convert(fT,sin(pi*x)) atol=eps(pi*convert(fT,x)) - @test cospi(convert(T,x))::fT ≈ convert(fT,cos(pi*x)) atol=eps(pi*convert(fT,x)) + @test sinpi(convert(T,x))::fT ≈ convert(fT,sin(pi*x)) atol=ulp(pi*convert(fT,x)) + @test cospi(convert(T,x))::fT ≈ convert(fT,cos(pi*x)) atol=ulp(pi*convert(fT,x)) end @test sinpi(convert(T,0.0))::fT === zero(fT) @@ -499,16 +499,16 @@ end y = log(xt) yb = log(big(xt)) - @test abs(y-yb) <= 0.56*eps(T(yb)) + @test abs(y-yb) <= 0.56*ulp(T(yb)) y = log1p(xt) yb = log1p(big(xt)) - @test abs(y-yb) <= 0.56*eps(T(yb)) + @test abs(y-yb) <= 0.56*ulp(T(yb)) if n <= 0 y = log1p(-xt) yb = log1p(big(-xt)) - @test abs(y-yb) <= 0.56*eps(T(yb)) + @test abs(y-yb) <= 0.56*ulp(T(yb)) end end end @@ -609,8 +609,8 @@ end # #22742: updated isapprox semantics @test !isapprox(1.0, 1.0+1e-12, atol=1e-14) -@test isapprox(1.0, 1.0+0.5*sqrt(eps(1.0))) -@test !isapprox(1.0, 1.0+1.5*sqrt(eps(1.0)), atol=sqrt(eps(1.0))) +@test isapprox(1.0, 1.0+0.5*sqrt(ulp(1.0))) +@test !isapprox(1.0, 1.0+1.5*sqrt(ulp(1.0)), atol=sqrt(ulp(1.0))) # test AbstractFloat fallback pr22716 struct Float22716{T<:AbstractFloat} <: AbstractFloat @@ -632,9 +632,9 @@ end @test asin(-one(T)) === -T(pi)/2 for x in (0.45, 0.6, 0.98) by = asin(big(T(x))) - @test T(abs(asin(T(x)) - by))/eps(T(abs(by))) <= 1 + @test T(abs(asin(T(x)) - by))/ulp(T(abs(by))) <= 1 bym = asin(big(T(-x))) - @test T(abs(asin(T(-x)) - bym))/eps(T(abs(bym))) <= 1 + @test T(abs(asin(T(-x)) - bym))/ulp(T(abs(bym))) <= 1 end @test_throws DomainError asin(-T(Inf)) @test_throws DomainError asin(T(Inf)) @@ -655,9 +655,9 @@ end for x in (0.1, 0.45, 0.6, 0.75, 0.79, 0.98) for op in (sin, cos, tan) by = T(op(big(x))) - @test abs(op(T(x)) - by)/eps(by) <= one(T) + @test abs(op(T(x)) - by)/ulp(by) <= one(T) bym = T(op(big(-x))) - @test abs(op(T(-x)) - bym)/eps(bym) <= one(T) + @test abs(op(T(-x)) - bym)/ulp(bym) <= one(T) end end @test_throws DomainError sin(-T(Inf)) @@ -699,13 +699,13 @@ end (T(39/16)+T(2)^23)/2, T(2)^23) x = T(7/16) by = T(atan(big(x))) - @test abs(atan(x) - by)/eps(by) <= one(T) + @test abs(atan(x) - by)/ulp(by) <= one(T) x = prevfloat(T(7/16)) by = T(atan(big(x))) - @test abs(atan(x) - by)/eps(by) <= one(T) + @test abs(atan(x) - by)/ulp(by) <= one(T) x = nextfloat(T(7/16)) by = T(atan(big(x))) - @test abs(atan(x) - by)/eps(by) <= one(T) + @test abs(atan(x) - by)/ulp(by) <= one(T) end # This case was used to find a bug, but it isn't special in itself @test atan(1.7581305072934137) ≈ 1.053644580517088 @@ -814,9 +814,9 @@ end @test acos(-one(T)) === T(pi) for x in (0.45, 0.6, 0.98) by = acos(big(T(x))) - @test T((acos(T(x)) - by))/eps(abs(T(by))) <= 1 + @test T((acos(T(x)) - by))/ulp(abs(T(by))) <= 1 bym = acos(big(T(-x))) - @test T(abs(acos(T(-x)) - bym))/eps(abs(T(bym))) <= 1 + @test T(abs(acos(T(-x)) - bym))/ulp(abs(T(bym))) <= 1 end @test_throws DomainError acos(-T(Inf)) @test_throws DomainError acos(T(Inf)) @@ -838,8 +838,8 @@ import Base.Math: COSH_SMALL_X, H_SMALL_X, H_MEDIUM_X, H_LARGE_X @test sinh(-T(1000)) === -T(Inf) @test isnan_type(T, sinh(T(NaN))) for x in Iterators.flatten(pcnfloat.([H_SMALL_X(T), H_MEDIUM_X(T), H_LARGE_X(T)])) - @test sinh(x) ≈ sinh(big(x)) rtol=eps(T) - @test sinh(-x) ≈ sinh(big(-x)) rtol=eps(T) + @test sinh(x) ≈ sinh(big(x)) rtol=ulp(T) + @test sinh(-x) ≈ sinh(big(-x)) rtol=ulp(T) end end end @@ -854,8 +854,8 @@ end @test cosh(-T(1000)) === T(Inf) @test isnan_type(T, cosh(T(NaN))) for x in Iterators.flatten(pcnfloat.([COSH_SMALL_X(T), H_MEDIUM_X(T), H_LARGE_X(T)])) - @test cosh(x) ≈ cosh(big(x)) rtol=eps(T) - @test cosh(-x) ≈ cosh(big(-x)) rtol=eps(T) + @test cosh(x) ≈ cosh(big(x)) rtol=ulp(T) + @test cosh(-x) ≈ cosh(big(-x)) rtol=ulp(T) end end end @@ -870,8 +870,8 @@ end @test tanh(-T(1000)) === -one(T) @test isnan_type(T, tanh(T(NaN))) for x in Iterators.flatten(pcnfloat.([H_SMALL_X(T), T(1.0), H_MEDIUM_X(T)])) - @test tanh(x) ≈ tanh(big(x)) rtol=eps(T) - @test tanh(-x) ≈ tanh(big(-x)) rtol=eps(T) + @test tanh(x) ≈ tanh(big(x)) rtol=ulp(T) + @test tanh(-x) ≈ tanh(big(-x)) rtol=ulp(T) end end end @@ -884,8 +884,8 @@ end @test asinh(prevfloat(zero(T))) === prevfloat(zero(T)) @test isnan_type(T, asinh(T(NaN))) for x in Iterators.flatten(pcnfloat.([T(2)^-28,T(2),T(2)^28])) - @test asinh(x) ≈ asinh(big(x)) rtol=eps(T) - @test asinh(-x) ≈ asinh(big(-x)) rtol=eps(T) + @test asinh(x) ≈ asinh(big(x)) rtol=ulp(T) + @test asinh(-x) ≈ asinh(big(-x)) rtol=ulp(T) end end end @@ -896,7 +896,7 @@ end @test acosh(one(T)) === zero(T) @test isnan_type(T, acosh(T(NaN))) for x in Iterators.flatten(pcnfloat.([nextfloat(T(1.0)), T(2), T(2)^28])) - @test acosh(x) ≈ acosh(big(x)) rtol=eps(T) + @test acosh(x) ≈ acosh(big(x)) rtol=ulp(T) end end end @@ -912,8 +912,8 @@ end @test atanh(prevfloat(zero(T))) === prevfloat(zero(T)) @test isnan_type(T, atanh(T(NaN))) for x in Iterators.flatten(pcnfloat.([T(2.0)^-28, T(0.5)])) - @test atanh(x) ≈ atanh(big(x)) rtol=eps(T) - @test atanh(-x) ≈ atanh(big(-x)) rtol=eps(T) + @test atanh(x) ≈ atanh(big(x)) rtol=ulp(T) + @test atanh(-x) ≈ atanh(big(-x)) rtol=ulp(T) end end end @@ -971,9 +971,9 @@ end map(x->x^3, 1.0:1.0:1024.0)..., nextfloat(-T(Inf)), prevfloat(T(Inf))) by = cbrt(big(T(x))) - @test cbrt(T(x)) ≈ by rtol=eps(T) + @test cbrt(T(x)) ≈ by rtol=ulp(T) bym = cbrt(big(T(-x))) - @test cbrt(T(-x)) ≈ bym rtol=eps(T) + @test cbrt(T(-x)) ≈ bym rtol=ulp(T) end end end diff --git a/test/mod2pi.jl b/test/mod2pi.jl index 5b0cb906bcef2..230639e66fd57 100644 --- a/test/mod2pi.jl +++ b/test/mod2pi.jl @@ -244,14 +244,14 @@ testModPi() # negative argument n, ret = Base.Math.rem_pio2_kernel(-case) ret_sum = ret.hi+ret.lo - ulp_error = (ret_sum-ieee754_rem_pio2_return[1, i])/eps(ieee754_rem_pio2_return[1, i]) + ulp_error = (ret_sum-ieee754_rem_pio2_return[1, i])/ulp(ieee754_rem_pio2_return[1, i]) @test ulp_error <= 0.5 diff = Float64(mod(big(-case), big(pi)/2))-(ret.hi+ret.lo) @test abs(diff) in (0.0, 1.5707963267948966, 1.5707963267948968) # positive argument n, ret = Base.Math.rem_pio2_kernel(case) ret_sum = ret.hi+ret.lo - ulp_error = (ret_sum-ieee754_rem_pio2_return[2, i])/eps(ieee754_rem_pio2_return[2, i]) + ulp_error = (ret_sum-ieee754_rem_pio2_return[2, i])/ulp(ieee754_rem_pio2_return[2, i]) @test ulp_error <= 0.5 diff = Float64(mod(big(case), big(pi)/2))-(ret.hi+ret.lo) @test abs(diff) in (0.0, 1.5707963267948966, 1.5707963267948968) diff --git a/test/mpfr.jl b/test/mpfr.jl index 1011b4522a6ec..ae8d2974eaa9c 100644 --- a/test/mpfr.jl +++ b/test/mpfr.jl @@ -643,10 +643,10 @@ end @test string(parse(BigFloat, "-9.9")) == "-9.8999999999999999999999999999999999997" end end -@testset "eps" begin - x = eps(BigFloat) +@testset "ulp" begin + x = ulp(BigFloat) @test BigFloat(1) + x == BigFloat(1) + prevfloat(x) - @test eps(BigFloat) == eps(BigFloat(1)) + @test ulp(BigFloat) == ulp(BigFloat(1)) end @testset "realmin/realmax" begin x = realmin(BigFloat) diff --git a/test/numbers.jl b/test/numbers.jl index 620ef4f09fb0d..bbc4901a81caa 100644 --- a/test/numbers.jl +++ b/test/numbers.jl @@ -1577,26 +1577,26 @@ end end # things related to floating-point epsilon -@test eps() == eps(Float64) -@test eps(Float64) == eps(1.0) -@test eps(Float64) == eps(1.5) -@test eps(Float32) == eps(1f0) -@test eps(float(0)) == 5e-324 -@test eps(-float(0)) == 5e-324 -@test eps(nextfloat(float(0))) == 5e-324 -@test eps(-nextfloat(float(0))) == 5e-324 -@test eps(realmin()) == 5e-324 -@test eps(-realmin()) == 5e-324 -@test eps(realmax()) == 2.0^(1023-52) -@test eps(-realmax()) == 2.0^(1023-52) -@test isnan(eps(NaN)) -@test isnan(eps(Inf)) -@test isnan(eps(-Inf)) +@test ulp() == ulp(Float64) +@test ulp(Float64) == ulp(1.0) +@test ulp(Float64) == ulp(1.5) +@test ulp(Float32) == ulp(1f0) +@test ulp(float(0)) == 5e-324 +@test ulp(-float(0)) == 5e-324 +@test ulp(nextfloat(float(0))) == 5e-324 +@test ulp(-nextfloat(float(0))) == 5e-324 +@test ulp(realmin()) == 5e-324 +@test ulp(-realmin()) == 5e-324 +@test ulp(realmax()) == 2.0^(1023-52) +@test ulp(-realmax()) == 2.0^(1023-52) +@test isnan(ulp(NaN)) +@test isnan(ulp(Inf)) +@test isnan(ulp(-Inf)) @test .1+.1+.1 != .3 @test .1+.1+.1 ≈ .3 @test .1+.1+.1-.3 ≉ 0 -@test .1+.1+.1-.3 ≈ 0 atol=eps(.3) +@test .1+.1+.1-.3 ≈ 0 atol=ulp(.3) @test 1.1 ≈ 1.1f0 @test div(1e50,1) == 1e50 @@ -1782,8 +1782,8 @@ end @test 0xf.fP1 === 31.875 @test -0x1.0p2 === -4.0 end -@testset "eps / realmin / realmax" begin - @test 0x1p-52 == eps() +@testset "ulp / realmin / realmax" begin + @test 0x1p-52 == ulp() @test 0x1p-52 + 1 != 1 @test 0x1p-53 + 1 == 1 @test 0x1p-1022 == realmin() @@ -1985,8 +1985,8 @@ for F in (Float16,Float32,Float64) @test reinterpret(Signed,one(F)) === signed(Base.exponent_one(F)) end -@test eps(realmax(Float64)) == 1.99584030953472e292 -@test eps(-realmax(Float64)) == 1.99584030953472e292 +@test ulp(realmax(Float64)) == 1.99584030953472e292 +@test ulp(-realmax(Float64)) == 1.99584030953472e292 # modular multiplicative inverses of odd numbers via exponentiation diff --git a/test/ranges.jl b/test/ranges.jl index cca10347fa065..6a1aa4b2f6292 100644 --- a/test/ranges.jl +++ b/test/ranges.jl @@ -75,7 +75,7 @@ end end for (I, T) in ((Int16, Float16), (Int32, Float32), (Int64, Float64)) x = T(typemax(I)) - Δi = ceil(I, eps(x)) + Δi = ceil(I, ulp(x)) for i = typemax(I)-2Δi:typemax(I)-Δi hi, lo = Base.splitprec(T, i) @test widen(hi) + widen(lo) == i @@ -645,7 +645,7 @@ end @testset "ranges with very small endpoints for type $T" for T = (Float32, Float64) z = zero(T) - u = eps(z) + u = ulp(z) @test first(range(u, stop=u, length=0)) == u @test last(range(u, stop=u, length=0)) == u @test first(range(-u, stop=u, length=0)) == -u @@ -715,7 +715,7 @@ end Rs = AbstractRange[1:1, 1:1:1, 1:2, 1:1:2, map(Int32,1:3:17), map(Int64,1:3:17), 1:0, 1:-1:0, 17:-3:0, 0.0:0.1:1.0, map(Float32,0.0:0.1:1.0), - 1.0:eps():1.0 .+ 10eps(), 9007199254740990.:1.0:9007199254740994, + 1.0:ulp():1.0 .+ 10ulp(), 9007199254740990.:1.0:9007199254740994, range(0, stop=1, length=20), map(Float32, range(0, stop=1, length=20))] for r in Rs local r @@ -792,7 +792,7 @@ end let r = range(1/3, stop=5/7, length=6) @test length(r) == 6 @test r[1] == 1/3 - @test abs(r[end] - 5/7) <= eps(5/7) + @test abs(r[end] - 5/7) <= ulp(5/7) end let r = range(0.25, stop=0.25, length=1) diff --git a/test/rational.jl b/test/rational.jl index 506136b5cfc5d..510f3467bec74 100644 --- a/test/rational.jl +++ b/test/rational.jl @@ -299,15 +299,15 @@ end @test rationalize(Int64, nextfloat(0.1)) == 300239975158034//3002399751580339 @test rationalize(Int128,nextfloat(0.1)) == 300239975158034//3002399751580339 @test rationalize(BigInt,nextfloat(0.1)) == 300239975158034//3002399751580339 - @test rationalize(Int8, nextfloat(0.1),tol=0.5eps(0.1)) == 1//10 - @test rationalize(Int64, nextfloat(0.1),tol=0.5eps(0.1)) == 379250494936463//3792504949364629 - @test rationalize(Int128,nextfloat(0.1),tol=0.5eps(0.1)) == 379250494936463//3792504949364629 - @test rationalize(BigInt,nextfloat(0.1),tol=0.5eps(0.1)) == 379250494936463//3792504949364629 - @test rationalize(Int8, nextfloat(0.1),tol=1.5eps(0.1)) == 1//10 - @test rationalize(Int64, nextfloat(0.1),tol=1.5eps(0.1)) == 1//10 - @test rationalize(Int128,nextfloat(0.1),tol=1.5eps(0.1)) == 1//10 - @test rationalize(BigInt,nextfloat(0.1),tol=1.5eps(0.1)) == 1//10 - @test rationalize(BigInt,nextfloat(parse(BigFloat,"0.1")),tol=1.5eps(big(0.1))) == 1//10 + @test rationalize(Int8, nextfloat(0.1),tol=0.5ulp(0.1)) == 1//10 + @test rationalize(Int64, nextfloat(0.1),tol=0.5ulp(0.1)) == 379250494936463//3792504949364629 + @test rationalize(Int128,nextfloat(0.1),tol=0.5ulp(0.1)) == 379250494936463//3792504949364629 + @test rationalize(BigInt,nextfloat(0.1),tol=0.5ulp(0.1)) == 379250494936463//3792504949364629 + @test rationalize(Int8, nextfloat(0.1),tol=1.5ulp(0.1)) == 1//10 + @test rationalize(Int64, nextfloat(0.1),tol=1.5ulp(0.1)) == 1//10 + @test rationalize(Int128,nextfloat(0.1),tol=1.5ulp(0.1)) == 1//10 + @test rationalize(BigInt,nextfloat(0.1),tol=1.5ulp(0.1)) == 1//10 + @test rationalize(BigInt,nextfloat(parse(BigFloat,"0.1")),tol=1.5ulp(big(0.1))) == 1//10 @test rationalize(Int64, nextfloat(0.1),tol=0) == 7205759403792795//72057594037927936 @test rationalize(Int128,nextfloat(0.1),tol=0) == 7205759403792795//72057594037927936 @test rationalize(BigInt,nextfloat(0.1),tol=0) == 7205759403792795//72057594037927936 diff --git a/test/rounding.jl b/test/rounding.jl index f9537180f0616..3a88a309dc8ef 100644 --- a/test/rounding.jl +++ b/test/rounding.jl @@ -54,7 +54,7 @@ end @test pn == pd || pn == pu @test v > 0 ? pz == pd : pz == pu - @test pu - pd == eps(pz) + @test pu - pd == ulp(pz) end for T in [Float32,Float64] @@ -73,7 +73,7 @@ end @test pn == pd || pn == pu @test v > 0 ? pz == pd : pz == pu - @test isinf(pu) || pu - pd == eps(pz) + @test isinf(pu) || pu - pd == ulp(pz) end end end