#37997

Committed 29 Jan 2025 02:08AM UTC coverage: 17.283% (-68.7%) from 85.981%

Build # #37997

Build Type

push

local

Committed by

web-flow

Commit Message

bpart: Start enforcing min_world for global variable definitions (#57150)

This is the analog of #57102 for global variables. Unlike for consants,
there is no automatic global backdate mechanism. The reasoning for this
is that global variables can be declared at any time, unlike constants
which can only be decalared once their value is available. As a result
code patterns using `Core.eval` to declare globals are rarer and likely
incorrect.

Run Details

1 of 22 new or added lines in 3 files covered. (4.55%)

31430 existing lines in 188 files now uncovered.

7903 of 45728 relevant lines covered (17.28%)

98663.7 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

3.7

/base/floatfuncs.jl

# This file is a part of Julia. License is MIT: https://julialang.org/license

## floating-point functions ##

copysign(x::Float64, y::Float64) = copysign_float(x, y)
copysign(x::Float32, y::Float32) = copysign_float(x, y)
copysign(x::Float32, y::Real) = copysign(x, Float32(y))
copysign(x::Float64, y::Real) = copysign(x, Float64(y))

flipsign(x::Float64, y::Float64) = bitcast(Float64, xor_int(bitcast(UInt64, x), and_int(bitcast(UInt64, y), 0x8000000000000000)))
flipsign(x::Float32, y::Float32) = bitcast(Float32, xor_int(bitcast(UInt32, x), and_int(bitcast(UInt32, y), 0x80000000)))
flipsign(x::Float32, y::Real) = flipsign(x, Float32(y))
flipsign(x::Float64, y::Real) = flipsign(x, Float64(y))

signbit(x::Float64) = signbit(bitcast(Int64, x))
signbit(x::Float32) = signbit(bitcast(Int32, x))
signbit(x::Float16) = signbit(bitcast(Int16, x))

"""
    maxintfloat(T=Float64)

The largest consecutive integer-valued floating-point number that is exactly represented in
the given floating-point type `T` (which defaults to `Float64`).

That is, `maxintfloat` returns the smallest positive integer-valued floating-point number
`n` such that `n+1` is *not* exactly representable in the type `T`.

When an `Integer`-type value is needed, use `Integer(maxintfloat(T))`.
"""
maxintfloat(::Type{Float64}) = 9007199254740992.
maxintfloat(::Type{Float32}) = Float32(16777216.)
maxintfloat(::Type{Float16}) = Float16(2048f0)
maxintfloat(x::T) where {T<:AbstractFloat} = maxintfloat(T)

"""
    maxintfloat(T, S)

The largest consecutive integer representable in the given floating-point type `T` that
also does not exceed the maximum integer representable by the integer type `S`.  Equivalently,
it is the minimum of `maxintfloat(T)` and [`typemax(S)`](@ref).
"""
maxintfloat(::Type{S}, ::Type{T}) where {S<:AbstractFloat, T<:Integer} = min(maxintfloat(S), S(typemax(T)))
maxintfloat() = maxintfloat(Float64)

isinteger(x::AbstractFloat) = iszero(x - trunc(x)) # note: x == trunc(x) would be incorrect for x=Inf

# See rounding.jl for docstring.

# NOTE: this relies on the current keyword dispatch behaviour (#9498).
function round(x::Real, r::RoundingMode=RoundNearest;
               digits::Union{Nothing,Integer}=nothing, sigdigits::Union{Nothing,Integer}=nothing, base::Union{Nothing,Integer}=nothing)
    if digits === nothing
        if sigdigits === nothing
            if base === nothing
                # avoid recursive calls
                throw(MethodError(round, (x,r)))
            else
                round(x,r)
                # or throw(ArgumentError("`round` cannot use `base` argument without `digits` or `sigdigits` arguments."))
            end
        else
            isfinite(x) || return float(x)
            _round_sigdigits(x, r, sigdigits, base === nothing ? 10 : base)
        end
    else
        if sigdigits === nothing
            isfinite(x) || return float(x)
            _round_digits(x, r, digits, base === nothing ? 10 : base)
        else
            throw(ArgumentError("`round` cannot use both `digits` and `sigdigits` arguments."))
        end
    end
end

# round x to multiples of 1/invstep
function _round_invstep(x, invstep, r::RoundingMode)
    y = round(x * invstep, r) / invstep
    if !isfinite(y)
        return x
    end
    return y
end

# round x to multiples of 1/(invstepsqrt^2)
# Using square root of step prevents overflowing
function _round_invstepsqrt(x, invstepsqrt, r::RoundingMode)
    y = round((x * invstepsqrt) * invstepsqrt, r) / invstepsqrt / invstepsqrt
    if !isfinite(y)
        return x
    end
    return y
end

# round x to multiples of step
function _round_step(x, step, r::RoundingMode)
    # TODO: use div with rounding mode
    y = round(x / step, r) * step
    if !isfinite(y)
        if x > 0
            return (r == RoundUp ? oftype(x, Inf) : zero(x))
        elseif x < 0
            return (r == RoundDown ? -oftype(x, Inf) : -zero(x))
        else
            return x
        end
    end
    return y
end

function _round_digits(x, r::RoundingMode, digits::Integer, base)
    fx = float(x)
    if digits >= 0
        invstep = oftype(fx, base)^digits
        if isfinite(invstep)
            return _round_invstep(fx, invstep, r)
        else
            invstepsqrt = oftype(fx, base)^oftype(fx, digits/2)
            return _round_invstepsqrt(fx, invstepsqrt, r)
        end
    else
        step = oftype(fx, base)^-digits
        return _round_step(fx, step, r)
    end
end

hidigit(x::Integer, base) = ndigits0z(x, base)
function hidigit(x::AbstractFloat, base)
    iszero(x) && return 0
    if base == 10
        return 1 + floor(Int, log10(abs(x)))
    elseif base == 2
        return 1 + exponent(x)
    else
        return 1 + floor(Int, log(base, abs(x)))
    end
end
hidigit(x::Real, base) = hidigit(float(x), base)

function _round_sigdigits(x, r::RoundingMode, sigdigits::Integer, base)
    h = hidigit(x, base)
    _round_digits(x, r, sigdigits-h, base)
end

# C-style round
function round(x::AbstractFloat, ::RoundingMode{:NearestTiesAway})
    y = trunc(x)
    ifelse(x==y,y,trunc(2*x-y))
end
# Java-style round
function round(x::T, ::RoundingMode{:NearestTiesUp}) where {T <: AbstractFloat}
    copysign(floor((x + (T(0.25) - eps(T(0.5)))) + (T(0.25) + eps(T(0.5)))), x)
end

function Base.round(x::AbstractFloat, ::typeof(RoundFromZero))
    signbit(x) ? round(x, RoundDown) : round(x, RoundUp)
end

# isapprox: approximate equality of numbers
"""
    isapprox(x, y; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps, nans::Bool=false[, norm::Function])

Inexact equality comparison. Two numbers compare equal if their relative distance *or* their
absolute distance is within tolerance bounds: `isapprox` returns `true` if
`norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))`. The default `atol` (absolute tolerance) is zero and the
default `rtol` (relative tolerance) 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).
This corresponds to requiring equality of about half of the significant digits. Otherwise,
e.g. for integer arguments or if an `atol > 0` is supplied, `rtol` defaults to zero.

The `norm` keyword defaults to `abs` for numeric `(x,y)` and to `LinearAlgebra.norm` for
arrays (where an alternative `norm` choice is sometimes useful).
When `x` and `y` are arrays, if `norm(x-y)` is not finite (i.e. `±Inf`
or `NaN`), the comparison falls back to checking whether all elements of `x` and `y` are
approximately equal component-wise.

The binary operator `≈` is equivalent to `isapprox` with the default arguments, and `x ≉ y`
is equivalent to `!isapprox(x,y)`.

Note that `x ≈ 0` (i.e., comparing to zero with the default tolerances) is
equivalent to `x == 0` since the default `atol` is `0`.  In such cases, you should either
supply an appropriate `atol` (or use `norm(x) ≤ atol`) or rearrange your code (e.g.
use `x ≈ y` rather than `x - y ≈ 0`).   It is not possible to pick a nonzero `atol`
automatically because it depends on the overall scaling (the "units") of your problem:
for example, in `x - y ≈ 0`, `atol=1e-9` is an absurdly small tolerance if `x` is the
[radius of the Earth](https://en.wikipedia.org/wiki/Earth_radius) in meters,
but an absurdly large tolerance if `x` is the
[radius of a Hydrogen atom](https://en.wikipedia.org/wiki/Bohr_radius) in meters.

!!! compat "Julia 1.6"
    Passing the `norm` keyword argument when comparing numeric (non-array) arguments
    requires Julia 1.6 or later.

# Examples
```jldoctest
julia> isapprox(0.1, 0.15; atol=0.05)
true

julia> isapprox(0.1, 0.15; rtol=0.34)
true

julia> isapprox(0.1, 0.15; rtol=0.33)
false

julia> 0.1 + 1e-10 ≈ 0.1
true

julia> 1e-10 ≈ 0
false

julia> isapprox(1e-10, 0, atol=1e-8)
true

julia> isapprox([10.0^9, 1.0], [10.0^9, 2.0]) # using `norm`
true
```
"""
function isapprox(x::Number, y::Number;
                  atol::Real=0, rtol::Real=rtoldefault(x,y,atol),
                  nans::Bool=false, norm::Function=abs)
    x′, y′ = promote(x, y) # to avoid integer overflow
    x == y ||
        (isfinite(x) && isfinite(y) && norm(x-y) <= max(atol, rtol*max(norm(x′), norm(y′)))) ||
         (nans && isnan(x) && isnan(y))
end

function isapprox(x::Integer, y::Integer;
                  atol::Real=0, rtol::Real=rtoldefault(x,y,atol),
                  nans::Bool=false, norm::Function=abs)
    if norm === abs && atol < 1 && rtol == 0
        return x == y
    else
        # We need to take the difference `max` - `min` when comparing unsigned integers.
        _x, _y = x < y ? (x, y) : (y, x)
        return norm(_y - _x) <= max(atol, rtol*max(norm(_x), norm(_y)))
    end
end

"""
    isapprox(x; kwargs...) / ≈(x; kwargs...)

Create a function that compares its argument to `x` using `≈`, i.e. a function equivalent to `y -> y ≈ x`.

The keyword arguments supported here are the same as those in the 2-argument `isapprox`.

!!! compat "Julia 1.5"
    This method requires Julia 1.5 or later.
"""
isapprox(y; kwargs...) = x -> isapprox(x, y; kwargs...)

const ≈ = isapprox
"""
    x ≉ y

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{<: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)))
    return atol > 0 ? zero(rtol) : rtol
end

# fused multiply-add

"""
    fma(x, y, z)

Computes `x*y+z` without rounding the intermediate result `x*y`. On some systems this is
significantly more expensive than `x*y+z`. `fma` is used to improve accuracy in certain
algorithms. See [`muladd`](@ref).
"""
function fma end
function fma_emulated(a::Float16, b::Float16, c::Float16)
    Float16(muladd(Float32(a), Float32(b), Float32(c))) #don't use fma if the hardware doesn't have it.
end
function fma_emulated(a::Float32, b::Float32, c::Float32)::Float32
    ab = Float64(a) * b
    res = ab+c
    reinterpret(UInt64, res)&0x1fff_ffff!=0x1000_0000 && return res
    # yes error compensation is necessary. It sucks
    reslo = abs(c)>abs(ab) ? ab-(res - c) : c-(res - ab)
    res = iszero(reslo) ? res : (signbit(reslo) ? prevfloat(res) : nextfloat(res))
    return res
end

""" Splits a Float64 into a hi bit and a low bit where the high bit has 27 trailing 0s and the low bit has 26 trailing 0s"""
@inline function splitbits(x::Float64)
    hi = reinterpret(Float64, reinterpret(UInt64, x) & 0xffff_ffff_f800_0000)
    return hi, x-hi
end

function twomul(a::Float64, b::Float64)
    ahi, alo = splitbits(a)
    bhi, blo = splitbits(b)
    abhi = a*b
    blohi, blolo = splitbits(blo)
    ablo = alo*blohi - (((abhi - ahi*bhi) - alo*bhi) - ahi*blo) + blolo*alo
    return abhi, ablo
end

function fma_emulated(a::Float64, b::Float64,c::Float64)
    abhi, ablo = @inline twomul(a,b)
    if !isfinite(abhi+c) || isless(abs(abhi), nextfloat(0x1p-969)) || issubnormal(a) || issubnormal(b)
        aandbfinite = isfinite(a) && isfinite(b)
        if !(isfinite(c) && aandbfinite)
            return aandbfinite ? c : abhi+c
        end
        (iszero(a) || iszero(b)) && return abhi+c
        # The checks above satisfy exponent's nothrow precondition
        bias = Math._exponent_finite_nonzero(a) + Math._exponent_finite_nonzero(b)
        c_denorm = ldexp(c, -bias)
        if isfinite(c_denorm)
            # rescale a and b to [1,2), equivalent to ldexp(a, -exponent(a))
            issubnormal(a) && (a *= 0x1p52)
            issubnormal(b) && (b *= 0x1p52)
            a = reinterpret(Float64, (reinterpret(UInt64, a) & ~Base.exponent_mask(Float64)) | Base.exponent_one(Float64))
            b = reinterpret(Float64, (reinterpret(UInt64, b) & ~Base.exponent_mask(Float64)) | Base.exponent_one(Float64))
            c = c_denorm
            abhi, ablo = twomul(a,b)
            # abhi <= 4 -> isfinite(r)      (α)
            r = abhi+c
            # s ≈ 0                         (β)
            s = (abs(abhi) > abs(c)) ? (abhi-r+c+ablo) : (c-r+abhi+ablo)
            # α ⩓ β -> isfinite(sumhi)      (γ)
            sumhi = r+s
            # If result is subnormal, ldexp will cause double rounding because subnormals have fewer mantisa bits.
            # As such, we need to check whether round to even would lead to double rounding and manually round sumhi to avoid it.
            if issubnormal(ldexp(sumhi, bias))
                sumlo = r-sumhi+s
                # finite: See γ
                # non-zero: If sumhi == ±0., then ldexp(sumhi, bias) == ±0,
                # so we don't take this branch.
                bits_lost = -bias-Math._exponent_finite_nonzero(sumhi)-1022
                sumhiInt = reinterpret(UInt64, sumhi)
                if (bits_lost != 1) ⊻ (sumhiInt&1 == 1)
                    sumhi = nextfloat(sumhi, cmp(sumlo,0))
                end
            end
            return ldexp(sumhi, bias)
        end
        isinf(abhi) && signbit(c) == signbit(a*b) && return abhi
        # fall through
    end
    r = abhi+c
    s = (abs(abhi) > abs(c)) ? (abhi-r+c+ablo) : (c-r+abhi+ablo)
    return r+s
end

# Disable LLVM's fma if it is incorrect, e.g. because LLVM falls back
# onto a broken system libm; if so, use a software emulated fma
@assume_effects :consistent function fma(x::T, y::T, z::T) where {T<:IEEEFloat}
    Core.Intrinsics.have_fma(T) ? fma_float(x,y,z) : fma_emulated(x,y,z)
end

# This is necessary at least on 32-bit Intel Linux, since fma_float may
# have called glibc, and some broken glibc fma implementations don't
# properly restore the rounding mode
Rounding.setrounding_raw(Float32, Rounding.JL_FE_TONEAREST)
Rounding.setrounding_raw(Float64, Rounding.JL_FE_TONEAREST)

1	# This file is a part of Julia. License is MIT: https://julialang.org/license
2
3	## floating-point functions ##
4
5	copysign(x::Float64, y::Float64) = copysign_float(x, y)	12✔
UNCOV 6	copysign(x::Float32, y::Float32) = copysign_float(x, y)	×
UNCOV 7	copysign(x::Float32, y::Real) = copysign(x, Float32(y))	×
UNCOV 8	copysign(x::Float64, y::Real) = copysign(x, Float64(y))	×
9
UNCOV 10	flipsign(x::Float64, y::Float64) = bitcast(Float64, xor_int(bitcast(UInt64, x), and_int(bitcast(UInt64, y), 0x8000000000000000)))	×
UNCOV 11	flipsign(x::Float32, y::Float32) = bitcast(Float32, xor_int(bitcast(UInt32, x), and_int(bitcast(UInt32, y), 0x80000000)))	×
UNCOV 12	flipsign(x::Float32, y::Real) = flipsign(x, Float32(y))	×
UNCOV 13	flipsign(x::Float64, y::Real) = flipsign(x, Float64(y))	×
14
UNCOV 15	signbit(x::Float64) = signbit(bitcast(Int64, x))	×
UNCOV 16	signbit(x::Float32) = signbit(bitcast(Int32, x))	×
UNCOV 17	signbit(x::Float16) = signbit(bitcast(Int16, x))	×
18
19	"""
20	maxintfloat(T=Float64)
21
22	The largest consecutive integer-valued floating-point number that is exactly represented in
23	the given floating-point type `T` (which defaults to `Float64`).
24
25	That is, `maxintfloat` returns the smallest positive integer-valued floating-point number
26	`n` such that `n+1` is not exactly representable in the type `T`.
27
28	When an `Integer`-type value is needed, use `Integer(maxintfloat(T))`.
29	"""
30	maxintfloat(::Type{Float64}) = 9007199254740992.	×
31	maxintfloat(::Type{Float32}) = Float32(16777216.)	×
32	maxintfloat(::Type{Float16}) = Float16(2048f0)	×
UNCOV 33	maxintfloat(x::T) where {T<:AbstractFloat} = maxintfloat(T)	×
34
35	"""
36	maxintfloat(T, S)
37
38	The largest consecutive integer representable in the given floating-point type `T` that
39	also does not exceed the maximum integer representable by the integer type `S`. Equivalently,
40	it is the minimum of `maxintfloat(T)` and [`typemax(S)`](@ref).
41	"""
UNCOV 42	maxintfloat(::Type{S}, ::Type{T}) where {S<:AbstractFloat, T<:Integer} = min(maxintfloat(S), S(typemax(T)))	×
43	maxintfloat() = maxintfloat(Float64)	×
44
45	isinteger(x::AbstractFloat) = iszero(x - trunc(x)) # note: x == trunc(x) would be incorrect for x=Inf	27✔
46
47	# See rounding.jl for docstring.
48
49	# NOTE: this relies on the current keyword dispatch behaviour (#9498).
50	function round(x::Real, r::RoundingMode=RoundNearest;	154✔
51	digits::Union{Nothing,Integer}=nothing, sigdigits::Union{Nothing,Integer}=nothing, base::Union{Nothing,Integer}=nothing)
52	if digits === nothing
53	if sigdigits === nothing
54	if base === nothing
55	# avoid recursive calls
56	throw(MethodError(round, (x,r)))
57	else
58	round(x,r)
59	# or throw(ArgumentError("`round` cannot use `base` argument without `digits` or `sigdigits` arguments."))
60	end
61	else
62	isfinite(x) \|\| return float(x)
63	_round_sigdigits(x, r, sigdigits, base === nothing ? 10 : base)
64	end
65	else
66	if sigdigits === nothing
67	isfinite(x) \|\| return float(x)	74✔
68	_round_digits(x, r, digits, base === nothing ? 10 : base)	74✔
69	else
70	throw(ArgumentError("`round` cannot use both `digits` and `sigdigits` arguments."))
71	end
72	end
73	end
74
75	# round x to multiples of 1/invstep
UNCOV 76	function _round_invstep(x, invstep, r::RoundingMode)	×
UNCOV 77	y = round(x * invstep, r) / invstep	×
UNCOV 78	if !isfinite(y)	×
79	return x	×
80	end
UNCOV 81	return y	×
82	end
83
84	# round x to multiples of 1/(invstepsqrt^2)
85	# Using square root of step prevents overflowing
86	function _round_invstepsqrt(x, invstepsqrt, r::RoundingMode)	×
UNCOV 87	y = round((x * invstepsqrt) * invstepsqrt, r) / invstepsqrt / invstepsqrt	×
UNCOV 88	if !isfinite(y)	×
UNCOV 89	return x	×
90	end
UNCOV 91	return y	×
92	end
93
94	# round x to multiples of step
95	function _round_step(x, step, r::RoundingMode)	×
96	# TODO: use div with rounding mode
UNCOV 97	y = round(x / step, r) * step	×
UNCOV 98	if !isfinite(y)	×
UNCOV 99	if x > 0	×
UNCOV 100	return (r == RoundUp ? oftype(x, Inf) : zero(x))	×
UNCOV 101	elseif x < 0	×
UNCOV 102	return (r == RoundDown ? -oftype(x, Inf) : -zero(x))	×
103	else
UNCOV 104	return x	×
105	end
106	end
UNCOV 107	return y	×
108	end
109
UNCOV 110	function _round_digits(x, r::RoundingMode, digits::Integer, base)	×
UNCOV 111	fx = float(x)	×
UNCOV 112	if digits >= 0	×
UNCOV 113	invstep = oftype(fx, base)^digits	×
UNCOV 114	if isfinite(invstep)	×
UNCOV 115	return _round_invstep(fx, invstep, r)	×
116	else
UNCOV 117	invstepsqrt = oftype(fx, base)^oftype(fx, digits/2)	×
UNCOV 118	return _round_invstepsqrt(fx, invstepsqrt, r)	×
119	end
120	else
UNCOV 121	step = oftype(fx, base)^-digits	×
UNCOV 122	return _round_step(fx, step, r)	×
123	end
124	end
125
126	hidigit(x::Integer, base) = ndigits0z(x, base)	×
UNCOV 127	function hidigit(x::AbstractFloat, base)	×
UNCOV 128	iszero(x) && return 0	×
UNCOV 129	if base == 10	×
UNCOV 130	return 1 + floor(Int, log10(abs(x)))	×
UNCOV 131	elseif base == 2	×
UNCOV 132	return 1 + exponent(x)	×
133	else
UNCOV 134	return 1 + floor(Int, log(base, abs(x)))	×
135	end
136	end
UNCOV 137	hidigit(x::Real, base) = hidigit(float(x), base)	×
138
UNCOV 139	function _round_sigdigits(x, r::RoundingMode, sigdigits::Integer, base)	×
UNCOV 140	h = hidigit(x, base)	×
UNCOV 141	_round_digits(x, r, sigdigits-h, base)	×
142	end
143
144	# C-style round
UNCOV 145	function round(x::AbstractFloat, ::RoundingMode{:NearestTiesAway})	×
UNCOV 146	y = trunc(x)	×
UNCOV 147	ifelse(x==y,y,trunc(2*x-y))	×
148	end
149	# Java-style round
UNCOV 150	function round(x::T, ::RoundingMode{:NearestTiesUp}) where {T <: AbstractFloat}	×
UNCOV 151	copysign(floor((x + (T(0.25) - eps(T(0.5)))) + (T(0.25) + eps(T(0.5)))), x)	×
152	end
153
UNCOV 154	function Base.round(x::AbstractFloat, ::typeof(RoundFromZero))	×
UNCOV 155	signbit(x) ? round(x, RoundDown) : round(x, RoundUp)	×
156	end
157
158	# isapprox: approximate equality of numbers
159	"""
160	isapprox(x, y; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps, nans::Bool=false[, norm::Function])
161
162	Inexact equality comparison. Two numbers compare equal if their relative distance or their
163	absolute distance is within tolerance bounds: `isapprox` returns `true` if
164	`norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))`. The default `atol` (absolute tolerance) is zero and the
165	default `rtol` (relative tolerance) depends on the types of `x` and `y`. The keyword argument `nans` determines
166	whether or not NaN values are considered equal (defaults to false).
167
168	For real or complex floating-point values, if an `atol > 0` is not specified, `rtol` defaults to
169	the square root of [`eps`](@ref) of the type of `x` or `y`, whichever is bigger (least precise).
170	This corresponds to requiring equality of about half of the significant digits. Otherwise,
171	e.g. for integer arguments or if an `atol > 0` is supplied, `rtol` defaults to zero.
172
173	The `norm` keyword defaults to `abs` for numeric `(x,y)` and to `LinearAlgebra.norm` for
174	arrays (where an alternative `norm` choice is sometimes useful).
175	When `x` and `y` are arrays, if `norm(x-y)` is not finite (i.e. `±Inf`
176	or `NaN`), the comparison falls back to checking whether all elements of `x` and `y` are
177	approximately equal component-wise.
178
179	The binary operator `≈` is equivalent to `isapprox` with the default arguments, and `x ≉ y`
180	is equivalent to `!isapprox(x,y)`.
181
182	Note that `x ≈ 0` (i.e., comparing to zero with the default tolerances) is
183	equivalent to `x == 0` since the default `atol` is `0`. In such cases, you should either
184	supply an appropriate `atol` (or use `norm(x) ≤ atol`) or rearrange your code (e.g.
185	use `x ≈ y` rather than `x - y ≈ 0`). It is not possible to pick a nonzero `atol`
186	automatically because it depends on the overall scaling (the "units") of your problem:
187	for example, in `x - y ≈ 0`, `atol=1e-9` is an absurdly small tolerance if `x` is the
188	[radius of the Earth](https://en.wikipedia.org/wiki/Earth_radius) in meters,
189	but an absurdly large tolerance if `x` is the
190	[radius of a Hydrogen atom](https://en.wikipedia.org/wiki/Bohr_radius) in meters.
191
192	!!! compat "Julia 1.6"
193	Passing the `norm` keyword argument when comparing numeric (non-array) arguments
194	requires Julia 1.6 or later.
195
196	# Examples
197	```jldoctest
198	julia> isapprox(0.1, 0.15; atol=0.05)
199	true
200
201	julia> isapprox(0.1, 0.15; rtol=0.34)
202	true
203
204	julia> isapprox(0.1, 0.15; rtol=0.33)
205	false
206
207	julia> 0.1 + 1e-10 ≈ 0.1
208	true
209
210	julia> 1e-10 ≈ 0
211	false
212
213	julia> isapprox(1e-10, 0, atol=1e-8)
214	true
215
216	julia> isapprox([10.0^9, 1.0], [10.0^9, 2.0]) # using `norm`
217	true
218	```
219	"""
UNCOV 220	function isapprox(x::Number, y::Number;	×
221	atol::Real=0, rtol::Real=rtoldefault(x,y,atol),
222	nans::Bool=false, norm::Function=abs)
UNCOV 223	x′, y′ = promote(x, y) # to avoid integer overflow	×
UNCOV 224	x == y \|\|	×
225	(isfinite(x) && isfinite(y) && norm(x-y) <= max(atol, rtol*max(norm(x′), norm(y′)))) \|\|
226	(nans && isnan(x) && isnan(y))
227	end
228
UNCOV 229	function isapprox(x::Integer, y::Integer;	×
230	atol::Real=0, rtol::Real=rtoldefault(x,y,atol),
231	nans::Bool=false, norm::Function=abs)
UNCOV 232	if norm === abs && atol < 1 && rtol == 0	×
UNCOV 233	return x == y	×
234	else
235	# We need to take the difference `max` - `min` when comparing unsigned integers.
UNCOV 236	_x, _y = x < y ? (x, y) : (y, x)	×
UNCOV 237	return norm(_y - _x) <= max(atol, rtol*max(norm(_x), norm(_y)))	×
238	end
239	end
240
241	"""
242	isapprox(x; kwargs...) / ≈(x; kwargs...)
243
244	Create a function that compares its argument to `x` using `≈`, i.e. a function equivalent to `y -> y ≈ x`.
245
246	The keyword arguments supported here are the same as those in the 2-argument `isapprox`.
247
248	!!! compat "Julia 1.5"
249	This method requires Julia 1.5 or later.
250	"""
UNCOV 251	isapprox(y; kwargs...) = x -> isapprox(x, y; kwargs...)	×
252
253	const ≈ = isapprox
254	"""
255	x ≉ y
256
257	This is equivalent to `!isapprox(x,y)` (see [`isapprox`](@ref)).
258	"""
UNCOV 259	≉(args...; kws...) = !≈(args...; kws...)	×
260
261	# default tolerance arguments
UNCOV 262	rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))	×
263	rtoldefault(::Type{<:Real}) = 0	×
UNCOV 264	function rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}, atol::Real) where {T<:Number,S<:Number}	×
UNCOV 265	rtol = max(rtoldefault(real(T)),rtoldefault(real(S)))	×
UNCOV 266	return atol > 0 ? zero(rtol) : rtol	×
267	end
268
269	# fused multiply-add
270
271	"""
272	fma(x, y, z)
273
274	Computes `xy+z` without rounding the intermediate result `xy`. On some systems this is
275	significantly more expensive than `x*y+z`. `fma` is used to improve accuracy in certain
276	algorithms. See [`muladd`](@ref).
277	"""
278	function fma end
UNCOV 279	function fma_emulated(a::Float16, b::Float16, c::Float16)	×
UNCOV 280	Float16(muladd(Float32(a), Float32(b), Float32(c))) #don't use fma if the hardware doesn't have it.	×
281	end
UNCOV 282	function fma_emulated(a::Float32, b::Float32, c::Float32)::Float32	×
UNCOV 283	ab = Float64(a) * b	×
UNCOV 284	res = ab+c	×
UNCOV 285	reinterpret(UInt64, res)&0x1fff_ffff!=0x1000_0000 && return res	×
286	# yes error compensation is necessary. It sucks
UNCOV 287	reslo = abs(c)>abs(ab) ? ab-(res - c) : c-(res - ab)	×
UNCOV 288	res = iszero(reslo) ? res : (signbit(reslo) ? prevfloat(res) : nextfloat(res))	×
UNCOV 289	return res	×
290	end
291
292	""" Splits a Float64 into a hi bit and a low bit where the high bit has 27 trailing 0s and the low bit has 26 trailing 0s"""
UNCOV 293	@inline function splitbits(x::Float64)	×
UNCOV 294	hi = reinterpret(Float64, reinterpret(UInt64, x) & 0xffff_ffff_f800_0000)	×
UNCOV 295	return hi, x-hi	×
296	end
297
UNCOV 298	function twomul(a::Float64, b::Float64)	×
UNCOV 299	ahi, alo = splitbits(a)	×
UNCOV 300	bhi, blo = splitbits(b)	×
UNCOV 301	abhi = a*b	×
UNCOV 302	blohi, blolo = splitbits(blo)	×
UNCOV 303	ablo = aloblohi - (((abhi - ahibhi) - alobhi) - ahiblo) + blolo*alo	×
304	return abhi, ablo	×
305	end
306
UNCOV 307	function fma_emulated(a::Float64, b::Float64,c::Float64)	×
UNCOV 308	abhi, ablo = @inline twomul(a,b)	×
UNCOV 309	if !isfinite(abhi+c) \|\| isless(abs(abhi), nextfloat(0x1p-969)) \|\| issubnormal(a) \|\| issubnormal(b)	×
UNCOV 310	aandbfinite = isfinite(a) && isfinite(b)	×
UNCOV 311	if !(isfinite(c) && aandbfinite)	×
UNCOV 312	return aandbfinite ? c : abhi+c	×
313	end
UNCOV 314	(iszero(a) \|\| iszero(b)) && return abhi+c	×
315	# The checks above satisfy exponent's nothrow precondition
UNCOV 316	bias = Math._exponent_finite_nonzero(a) + Math._exponent_finite_nonzero(b)	×
UNCOV 317	c_denorm = ldexp(c, -bias)	×
UNCOV 318	if isfinite(c_denorm)	×
319	# rescale a and b to [1,2), equivalent to ldexp(a, -exponent(a))
UNCOV 320	issubnormal(a) && (a *= 0x1p52)	×
UNCOV 321	issubnormal(b) && (b *= 0x1p52)	×
UNCOV 322	a = reinterpret(Float64, (reinterpret(UInt64, a) & ~Base.exponent_mask(Float64)) \| Base.exponent_one(Float64))	×
UNCOV 323	b = reinterpret(Float64, (reinterpret(UInt64, b) & ~Base.exponent_mask(Float64)) \| Base.exponent_one(Float64))	×
324	c = c_denorm	×
UNCOV 325	abhi, ablo = twomul(a,b)	×
326	# abhi <= 4 -> isfinite(r) (α)
UNCOV 327	r = abhi+c	×
328	# s ≈ 0 (β)
UNCOV 329	s = (abs(abhi) > abs(c)) ? (abhi-r+c+ablo) : (c-r+abhi+ablo)	×
330	# α ⩓ β -> isfinite(sumhi) (γ)
UNCOV 331	sumhi = r+s	×
332	# If result is subnormal, ldexp will cause double rounding because subnormals have fewer mantisa bits.
333	# As such, we need to check whether round to even would lead to double rounding and manually round sumhi to avoid it.
UNCOV 334	if issubnormal(ldexp(sumhi, bias))	×
UNCOV 335	sumlo = r-sumhi+s	×
336	# finite: See γ
337	# non-zero: If sumhi == ±0., then ldexp(sumhi, bias) == ±0,
338	# so we don't take this branch.
UNCOV 339	bits_lost = -bias-Math._exponent_finite_nonzero(sumhi)-1022	×
UNCOV 340	sumhiInt = reinterpret(UInt64, sumhi)	×
UNCOV 341	if (bits_lost != 1) ⊻ (sumhiInt&1 == 1)	×
UNCOV 342	sumhi = nextfloat(sumhi, cmp(sumlo,0))	×
343	end
344	end
UNCOV 345	return ldexp(sumhi, bias)	×
346	end
UNCOV 347	isinf(abhi) && signbit(c) == signbit(a*b) && return abhi	×
348	# fall through
349	end
UNCOV 350	r = abhi+c	×
UNCOV 351	s = (abs(abhi) > abs(c)) ? (abhi-r+c+ablo) : (c-r+abhi+ablo)	×
UNCOV 352	return r+s	×
353	end
354
355	# Disable LLVM's fma if it is incorrect, e.g. because LLVM falls back
356	# onto a broken system libm; if so, use a software emulated fma
UNCOV 357	@assume_effects :consistent function fma(x::T, y::T, z::T) where {T<:IEEEFloat}	×
UNCOV 358	Core.Intrinsics.have_fma(T) ? fma_float(x,y,z) : fma_emulated(x,y,z)	×
359	end
360
361	# This is necessary at least on 32-bit Intel Linux, since fma_float may
362	# have called glibc, and some broken glibc fma implementations don't
363	# properly restore the rounding mode
364	Rounding.setrounding_raw(Float32, Rounding.JL_FE_TONEAREST)
365	Rounding.setrounding_raw(Float64, Rounding.JL_FE_TONEAREST)

JuliaLang / julia / #37997

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous