#37999

Committed 02 Feb 2025 07:22AM UTC coverage: 17.218% (-8.3%) from 25.515%

Build # #37999

Build Type

push

local

Committed by

web-flow

Commit Message

bpart: Start tracking backedges for bindings (#57213)

This PR adds limited backedge support for Bindings. There are two
classes of bindings that get backedges:

1. Cross-module `GlobalRef` bindings (new in this PR)
2. Any globals accesses through intrinsics (i.e. those with forward
edges from #57009)

This is a time/space trade-off for invalidation. As a result of the
first category, invalidating a binding now only needs to scan all the
methods defined in the same module as the binding. At the same time, it
is anticipated that most binding references are to bindings in the same
module, keeping the list of bindings that need explicit (back)edges
small.

Run Details

7 of 30 new or added lines in 3 files covered. (23.33%)

4235 existing lines in 124 files now uncovered.

7882 of 45779 relevant lines covered (17.22%)

98289.89 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✔
6	copysign(x::Float32, y::Float32) = copysign_float(x, y)	×
7	copysign(x::Float32, y::Real) = copysign(x, Float32(y))	×
8	copysign(x::Float64, y::Real) = copysign(x, Float64(y))	×
9
10	flipsign(x::Float64, y::Float64) = bitcast(Float64, xor_int(bitcast(UInt64, x), and_int(bitcast(UInt64, y), 0x8000000000000000)))	×
11	flipsign(x::Float32, y::Float32) = bitcast(Float32, xor_int(bitcast(UInt32, x), and_int(bitcast(UInt32, y), 0x80000000)))	×
12	flipsign(x::Float32, y::Real) = flipsign(x, Float32(y))	×
13	flipsign(x::Float64, y::Real) = flipsign(x, Float64(y))	×
14
UNCOV 15	signbit(x::Float64) = signbit(bitcast(Int64, x))	×
16	signbit(x::Float32) = signbit(bitcast(Int32, x))	×
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)	×
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	"""
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)	×
77	y = round(x * invstep, r) / invstep	×
78	if !isfinite(y)	×
79	return x	×
80	end
81	return y	×
82	end
83
84	# round x to multiples of 1/(invstepsqrt^2)
85	# Using square root of step prevents overflowing
UNCOV 86	function _round_invstepsqrt(x, invstepsqrt, r::RoundingMode)	×
87	y = round((x * invstepsqrt) * invstepsqrt, r) / invstepsqrt / invstepsqrt	×
88	if !isfinite(y)	×
89	return x	×
90	end
91	return y	×
92	end
93
94	# round x to multiples of step
UNCOV 95	function _round_step(x, step, r::RoundingMode)	×
96	# TODO: use div with rounding mode
97	y = round(x / step, r) * step	×
98	if !isfinite(y)	×
99	if x > 0	×
100	return (r == RoundUp ? oftype(x, Inf) : zero(x))	×
101	elseif x < 0	×
102	return (r == RoundDown ? -oftype(x, Inf) : -zero(x))	×
103	else
104	return x	×
105	end
106	end
107	return y	×
108	end
109
110	function _round_digits(x, r::RoundingMode, digits::Integer, base)	×
111	fx = float(x)	×
112	if digits >= 0	×
113	invstep = oftype(fx, base)^digits	×
114	if isfinite(invstep)	×
115	return _round_invstep(fx, invstep, r)	×
116	else
117	invstepsqrt = oftype(fx, base)^oftype(fx, digits/2)	×
118	return _round_invstepsqrt(fx, invstepsqrt, r)	×
119	end
120	else
121	step = oftype(fx, base)^-digits	×
122	return _round_step(fx, step, r)	×
123	end
124	end
125
126	hidigit(x::Integer, base) = ndigits0z(x, base)	×
127	function hidigit(x::AbstractFloat, base)	×
128	iszero(x) && return 0	×
129	if base == 10	×
130	return 1 + floor(Int, log10(abs(x)))	×
131	elseif base == 2	×
132	return 1 + exponent(x)	×
133	else
134	return 1 + floor(Int, log(base, abs(x)))	×
135	end
136	end
137	hidigit(x::Real, base) = hidigit(float(x), base)	×
138
139	function _round_sigdigits(x, r::RoundingMode, sigdigits::Integer, base)	×
140	h = hidigit(x, base)	×
141	_round_digits(x, r, sigdigits-h, base)	×
142	end
143
144	# C-style round
145	function round(x::AbstractFloat, ::RoundingMode{:NearestTiesAway})	×
146	y = trunc(x)	×
147	ifelse(x==y,y,trunc(2*x-y))	×
148	end
149	# Java-style round
150	function round(x::T, ::RoundingMode{:NearestTiesUp}) where {T <: AbstractFloat}	×
151	copysign(floor((x + (T(0.25) - eps(T(0.5)))) + (T(0.25) + eps(T(0.5)))), x)	×
152	end
153
154	function Base.round(x::AbstractFloat, ::typeof(RoundFromZero))	×
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	"""
220	function isapprox(x::Number, y::Number;	×
221	atol::Real=0, rtol::Real=rtoldefault(x,y,atol),
222	nans::Bool=false, norm::Function=abs)
223	x′, y′ = promote(x, y) # to avoid integer overflow	×
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
229	function isapprox(x::Integer, y::Integer;	×
230	atol::Real=0, rtol::Real=rtoldefault(x,y,atol),
231	nans::Bool=false, norm::Function=abs)
232	if norm === abs && atol < 1 && rtol == 0	×
233	return x == y	×
234	else
235	# We need to take the difference `max` - `min` when comparing unsigned integers.
236	_x, _y = x < y ? (x, y) : (y, x)	×
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	"""
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	"""
259	≉(args...; kws...) = !≈(args...; kws...)	×
260
261	# default tolerance arguments
262	rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))	×
263	rtoldefault(::Type{<:Real}) = 0	×
264	function rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}, atol::Real) where {T<:Number,S<:Number}	×
265	rtol = max(rtoldefault(real(T)),rtoldefault(real(S)))	×
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
279	function fma_emulated(a::Float16, b::Float16, c::Float16)	×
280	Float16(muladd(Float32(a), Float32(b), Float32(c))) #don't use fma if the hardware doesn't have it.	×
281	end
282	function fma_emulated(a::Float32, b::Float32, c::Float32)::Float32	×
283	ab = Float64(a) * b	×
284	res = ab+c	×
285	reinterpret(UInt64, res)&0x1fff_ffff!=0x1000_0000 && return res	×
286	# yes error compensation is necessary. It sucks
287	reslo = abs(c)>abs(ab) ? ab-(res - c) : c-(res - ab)	×
288	res = iszero(reslo) ? res : (signbit(reslo) ? prevfloat(res) : nextfloat(res))	×
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)	×
294	hi = reinterpret(Float64, reinterpret(UInt64, x) & 0xffff_ffff_f800_0000)	×
295	return hi, x-hi	×
296	end
297
UNCOV 298	function twomul(a::Float64, b::Float64)	×
299	ahi, alo = splitbits(a)	×
300	bhi, blo = splitbits(b)	×
301	abhi = a*b	×
302	blohi, blolo = splitbits(blo)	×
303	ablo = aloblohi - (((abhi - ahibhi) - alobhi) - ahiblo) + blolo*alo	×
304	return abhi, ablo	×
305	end
306
307	function fma_emulated(a::Float64, b::Float64,c::Float64)	×
308	abhi, ablo = @inline twomul(a,b)	×
309	if !isfinite(abhi+c) \|\| isless(abs(abhi), nextfloat(0x1p-969)) \|\| issubnormal(a) \|\| issubnormal(b)	×
310	aandbfinite = isfinite(a) && isfinite(b)	×
311	if !(isfinite(c) && aandbfinite)	×
312	return aandbfinite ? c : abhi+c	×
313	end
314	(iszero(a) \|\| iszero(b)) && return abhi+c	×
315	# The checks above satisfy exponent's nothrow precondition
316	bias = Math._exponent_finite_nonzero(a) + Math._exponent_finite_nonzero(b)	×
317	c_denorm = ldexp(c, -bias)	×
318	if isfinite(c_denorm)	×
319	# rescale a and b to [1,2), equivalent to ldexp(a, -exponent(a))
320	issubnormal(a) && (a *= 0x1p52)	×
321	issubnormal(b) && (b *= 0x1p52)	×
322	a = reinterpret(Float64, (reinterpret(UInt64, a) & ~Base.exponent_mask(Float64)) \| Base.exponent_one(Float64))	×
323	b = reinterpret(Float64, (reinterpret(UInt64, b) & ~Base.exponent_mask(Float64)) \| Base.exponent_one(Float64))	×
324	c = c_denorm	×
325	abhi, ablo = twomul(a,b)	×
326	# abhi <= 4 -> isfinite(r) (α)
327	r = abhi+c	×
328	# s ≈ 0 (β)
329	s = (abs(abhi) > abs(c)) ? (abhi-r+c+ablo) : (c-r+abhi+ablo)	×
330	# α ⩓ β -> isfinite(sumhi) (γ)
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.
334	if issubnormal(ldexp(sumhi, bias))	×
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.
339	bits_lost = -bias-Math._exponent_finite_nonzero(sumhi)-1022	×
340	sumhiInt = reinterpret(UInt64, sumhi)	×
341	if (bits_lost != 1) ⊻ (sumhiInt&1 == 1)	×
342	sumhi = nextfloat(sumhi, cmp(sumlo,0))	×
343	end
344	end
345	return ldexp(sumhi, bias)	×
346	end
347	isinf(abhi) && signbit(c) == signbit(a*b) && return abhi	×
348	# fall through
349	end
350	r = abhi+c	×
351	s = (abs(abhi) > abs(c)) ? (abhi-r+c+ablo) : (c-r+abhi+ablo)	×
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}	×
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 / #37999

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