#37886

Committed 28 Aug 2024 07:22AM UTC coverage: 87.854% (+3.1%) from 84.78%

Build # #37886

Build Type

push

local

Committed by

web-flow

Commit Message

Use native tls model in macos for better performance (#55576)

Macos has a native tls implementation in clang since at least clang 3.7
which much older than what we require so lets enable it for some small
performance gains.

We may want to turn on the ifunc optimization that's there as well but I
haven't tested it and it's only in clang 18 and up so it would be off
for most since Apple clang is 16 on their latest beta
https://github.com/llvm/llvm-project/pull/73687

Run Details

78230 of 89045 relevant lines covered (87.85%)

16853037.77 hits per line

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

92.96

/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
        return norm(x - y) <= 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::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
fma_llvm(x::Float32, y::Float32, z::Float32) = fma_float(x, y, z)
fma_llvm(x::Float64, y::Float64, z::Float64) = fma_float(x, y, z)

# 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 fma(x::Float32, y::Float32, z::Float32) = Core.Intrinsics.have_fma(Float32) ? fma_llvm(x,y,z) : fma_emulated(x,y,z)
@assume_effects :consistent fma(x::Float64, y::Float64, z::Float64) = Core.Intrinsics.have_fma(Float64) ? fma_llvm(x,y,z) : fma_emulated(x,y,z)

function fma(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

# This is necessary at least on 32-bit Intel Linux, since fma_llvm 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)	4,787,272✔
6	copysign(x::Float32, y::Float32) = copysign_float(x, y)	1,566,447✔
7	copysign(x::Float32, y::Real) = copysign(x, Float32(y))	51✔
8	copysign(x::Float64, y::Real) = copysign(x, Float64(y))	173,722✔
9
10	flipsign(x::Float64, y::Float64) = bitcast(Float64, xor_int(bitcast(UInt64, x), and_int(bitcast(UInt64, y), 0x8000000000000000)))	212,006✔
11	flipsign(x::Float32, y::Float32) = bitcast(Float32, xor_int(bitcast(UInt32, x), and_int(bitcast(UInt32, y), 0x80000000)))	123✔
12	flipsign(x::Float32, y::Real) = flipsign(x, Float32(y))	63✔
13	flipsign(x::Float64, y::Real) = flipsign(x, Float64(y))	96,575✔
14
15	signbit(x::Float64) = signbit(bitcast(Int64, x))	50,125,497✔
16	signbit(x::Float32) = signbit(bitcast(Int32, x))	18,757,379✔
17	signbit(x::Float16) = signbit(bitcast(Int16, x))	739,218✔
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.)	3,126,142✔
32	maxintfloat(::Type{Float16}) = Float16(2048f0)	2,284,993✔
33	maxintfloat(x::T) where {T<:AbstractFloat} = maxintfloat(T)	63✔
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)))	6,263,338✔
43	maxintfloat() = maxintfloat(Float64)	1✔
44
45	isinteger(x::AbstractFloat) = iszero(x - trunc(x)) # note: x == trunc(x) would be incorrect for x=Inf	36,843,415✔
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;	2,022,393✔
51	digits::Union{Nothing,Integer}=nothing, sigdigits::Union{Nothing,Integer}=nothing, base::Union{Nothing,Integer}=nothing)
52	if digits === nothing	280✔
53	if sigdigits === nothing	47✔
54	if base === nothing	1✔
55	# avoid recursive calls
56	throw(MethodError(round, (x,r)))	×
57	else
58	round(x,r)	1✔
59	# or throw(ArgumentError("`round` cannot use `base` argument without `digits` or `sigdigits` arguments."))
60	end
61	else
62	isfinite(x) \|\| return float(x)	48✔
63	_round_sigdigits(x, r, sigdigits, base === nothing ? 10 : base)	42✔
64	end
65	else
66	if sigdigits === nothing	233✔
67	isfinite(x) \|\| return float(x)	41,291✔
68	_round_digits(x, r, digits, base === nothing ? 10 : base)	41,283✔
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
76	function _round_invstep(x, invstep, r::RoundingMode)	6✔
77	y = round(x * invstep, r) / invstep	41,173✔
78	if !isfinite(y)	41,173✔
79	return x	×
80	end
81	return y	41,173✔
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)	×
87	y = round((x * invstepsqrt) * invstepsqrt, r) / invstepsqrt / invstepsqrt	26✔
88	if !isfinite(y)	26✔
89	return x	16✔
90	end
91	return y	10✔
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
97	y = round(x / step, r) * step	15✔
98	if !isfinite(y)	15✔
99	if x > 0	12✔
100	return (r == RoundUp ? oftype(x, Inf) : zero(x))	4✔
101	elseif x < 0	8✔
102	return (r == RoundDown ? -oftype(x, Inf) : -zero(x))	4✔
103	else
104	return x	4✔
105	end
106	end
107	return y	3✔
108	end
109
110	function _round_digits(x, r::RoundingMode, digits::Integer, base)	41,208✔
111	fx = float(x)	157✔
112	if digits >= 0	41,214✔
113	invstep = oftype(fx, base)^digits	82,377✔
114	if isfinite(invstep)	41,199✔
115	return _round_invstep(fx, invstep, r)	41,173✔
116	else
117	invstepsqrt = oftype(fx, base)^oftype(fx, digits/2)	26✔
118	return _round_invstepsqrt(fx, invstepsqrt, r)	26✔
119	end
120	else
121	step = oftype(fx, base)^-digits	24✔
122	return _round_step(fx, step, r)	15✔
123	end
124	end
125
126	hidigit(x::Integer, base) = ndigits0z(x, base)	×
127	function hidigit(x::AbstractFloat, base)	40✔
128	iszero(x) && return 0	42✔
129	if base == 10	40✔
130	return 1 + floor(Int, log10(abs(x)))	36✔
131	elseif base == 2	4✔
132	return 1 + exponent(x)	3✔
133	else
134	return 1 + floor(Int, log(base, abs(x)))	1✔
135	end
136	end
137	hidigit(x::Real, base) = hidigit(float(x), base)	3✔
138
139	function _round_sigdigits(x, r::RoundingMode, sigdigits::Integer, base)
140	h = hidigit(x, base)	42✔
141	_round_digits(x, r, sigdigits-h, base)	42✔
142	end
143
144	# C-style round
145	function round(x::AbstractFloat, ::RoundingMode{:NearestTiesAway})	49,520✔
146	y = trunc(x)	50,594✔
147	ifelse(x==y,y,trunc(2*x-y))	50,594✔
148	end
149	# Java-style round
150	function round(x::T, ::RoundingMode{:NearestTiesUp}) where {T <: AbstractFloat}	49,520✔
151	copysign(floor((x + (T(0.25) - eps(T(0.5)))) + (T(0.25) + eps(T(0.5)))), x)	50,594✔
152	end
153
154	function Base.round(x::AbstractFloat, ::typeof(RoundFromZero))	49,520✔
155	signbit(x) ? round(x, RoundDown) : round(x, RoundUp)	51,076✔
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;	386,270✔
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	368,877✔
224	x == y \|\|	385,920✔
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;	701✔
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	491✔
233	return x == y	437✔
234	else
235	return norm(x - y) <= max(atol, rtol*max(norm(x), norm(y)))	54✔
236	end
237	end
238
239	"""
240	isapprox(x; kwargs...) / ≈(x; kwargs...)
241
242	Create a function that compares its argument to `x` using `≈`, i.e. a function equivalent to `y -> y ≈ x`.
243
244	The keyword arguments supported here are the same as those in the 2-argument `isapprox`.
245
246	!!! compat "Julia 1.5"
247	This method requires Julia 1.5 or later.
248	"""
249	isapprox(y; kwargs...) = x -> isapprox(x, y; kwargs...)	197✔
250
251	const ≈ = isapprox
252	"""
253	x ≉ y
254
255	This is equivalent to `!isapprox(x,y)` (see [`isapprox`](@ref)).
256	"""
257	≉(args...; kws...) = !≈(args...; kws...)	58✔
258
259	# default tolerance arguments
260	rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))	430,552✔
261	rtoldefault(::Type{<:Real}) = 0	×
262	function rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}, atol::Real) where {T<:Number,S<:Number}	19,277✔
263	rtol = max(rtoldefault(real(T)),rtoldefault(real(S)))	427,392✔
264	return atol > 0 ? zero(rtol) : rtol	411,996✔
265	end
266
267	# fused multiply-add
268
269	"""
270	fma(x, y, z)
271
272	Computes `xy+z` without rounding the intermediate result `xy`. On some systems this is
273	significantly more expensive than `x*y+z`. `fma` is used to improve accuracy in certain
274	algorithms. See [`muladd`](@ref).
275	"""
276	function fma end
277	function fma_emulated(a::Float32, b::Float32, c::Float32)::Float32	262,155✔
278	ab = Float64(a) * b	262,155✔
279	res = ab+c	262,155✔
280	reinterpret(UInt64, res)&0x1fff_ffff!=0x1000_0000 && return res	262,155✔
281	# yes error compensation is necessary. It sucks
282	reslo = abs(c)>abs(ab) ? ab-(res - c) : c-(res - ab)	4✔
283	res = iszero(reslo) ? res : (signbit(reslo) ? prevfloat(res) : nextfloat(res))	3✔
284	return res	2✔
285	end
286
287	""" 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"""
288	@inline function splitbits(x::Float64)
289	hi = reinterpret(Float64, reinterpret(UInt64, x) & 0xffff_ffff_f800_0000)	923,691✔
290	return hi, x-hi	923,691✔
291	end
292
293	function twomul(a::Float64, b::Float64)
294	ahi, alo = splitbits(a)	307,897✔
295	bhi, blo = splitbits(b)	307,897✔
296	abhi = a*b	307,897✔
297	blohi, blolo = splitbits(blo)	307,897✔
298	ablo = aloblohi - (((abhi - ahibhi) - alobhi) - ahiblo) + blolo*alo	307,897✔
299	return abhi, ablo	×
300	end
301
302	function fma_emulated(a::Float64, b::Float64,c::Float64)	262,456✔
303	abhi, ablo = @inline twomul(a,b)	262,456✔
304	if !isfinite(abhi+c) \|\| isless(abs(abhi), nextfloat(0x1p-969)) \|\| issubnormal(a) \|\| issubnormal(b)	491,970✔
305	aandbfinite = isfinite(a) && isfinite(b)	69,541✔
306	if !(isfinite(c) && aandbfinite)	69,541✔
307	return aandbfinite ? c : abhi+c	440✔
308	end
309	(iszero(a) \|\| iszero(b)) && return abhi+c	69,101✔
310	# The checks above satisfy exponent's nothrow precondition
311	bias = Math._exponent_finite_nonzero(a) + Math._exponent_finite_nonzero(b)	69,100✔
312	c_denorm = ldexp(c, -bias)	138,200✔
313	if isfinite(c_denorm)	69,100✔
314	# rescale a and b to [1,2), equivalent to ldexp(a, -exponent(a))
315	issubnormal(a) && (a *= 0x1p52)	45,441✔
316	issubnormal(b) && (b *= 0x1p52)	45,441✔
317	a = reinterpret(Float64, (reinterpret(UInt64, a) & ~Base.exponent_mask(Float64)) \| Base.exponent_one(Float64))	45,441✔
318	b = reinterpret(Float64, (reinterpret(UInt64, b) & ~Base.exponent_mask(Float64)) \| Base.exponent_one(Float64))	45,441✔
319	c = c_denorm	×
320	abhi, ablo = twomul(a,b)	45,441✔
321	# abhi <= 4 -> isfinite(r) (α)
322	r = abhi+c	45,441✔
323	# s ≈ 0 (β)
324	s = (abs(abhi) > abs(c)) ? (abhi-r+c+ablo) : (c-r+abhi+ablo)	58,296✔
325	# α ⩓ β -> isfinite(sumhi) (γ)
326	sumhi = r+s	45,441✔
327	# If result is subnormal, ldexp will cause double rounding because subnormals have fewer mantisa bits.
328	# As such, we need to check whether round to even would lead to double rounding and manually round sumhi to avoid it.
329	if issubnormal(ldexp(sumhi, bias))	90,882✔
330	sumlo = r-sumhi+s	14✔
331	# finite: See γ
332	# non-zero: If sumhi == ±0., then ldexp(sumhi, bias) == ±0,
333	# so we don't take this branch.
334	bits_lost = -bias-Math._exponent_finite_nonzero(sumhi)-1022	14✔
335	sumhiInt = reinterpret(UInt64, sumhi)	14✔
336	if (bits_lost != 1) ⊻ (sumhiInt&1 == 1)	14✔
337	sumhi = nextfloat(sumhi, cmp(sumlo,0))	16✔
338	end
339	end
340	return ldexp(sumhi, bias)	45,441✔
341	end
342	isinf(abhi) && signbit(c) == signbit(a*b) && return abhi	23,659✔
343	# fall through
344	end
345	r = abhi+c	216,574✔
346	s = (abs(abhi) > abs(c)) ? (abhi-r+c+ablo) : (c-r+abhi+ablo)	334,984✔
347	return r+s	216,574✔
348	end
349	fma_llvm(x::Float32, y::Float32, z::Float32) = fma_float(x, y, z)	2,448,598✔
350	fma_llvm(x::Float64, y::Float64, z::Float64) = fma_float(x, y, z)	408,598,855✔
351
352	# Disable LLVM's fma if it is incorrect, e.g. because LLVM falls back
353	# onto a broken system libm; if so, use a software emulated fma
354	@assume_effects :consistent fma(x::Float32, y::Float32, z::Float32) = Core.Intrinsics.have_fma(Float32) ? fma_llvm(x,y,z) : fma_emulated(x,y,z)	2,448,598✔
355	@assume_effects :consistent fma(x::Float64, y::Float64, z::Float64) = Core.Intrinsics.have_fma(Float64) ? fma_llvm(x,y,z) : fma_emulated(x,y,z)	408,598,855✔
356
357	function fma(a::Float16, b::Float16, c::Float16)	1✔
358	Float16(muladd(Float32(a), Float32(b), Float32(c))) #don't use fma if the hardware doesn't have it.	499,718✔
359	end
360
361	# This is necessary at least on 32-bit Intel Linux, since fma_llvm 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 / #37886

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