#37918

Committed 28 Sep 2024 11:02PM UTC coverage: 86.484% (-0.8%) from 87.243%

Build # #37918

Build Type

push

local

Committed by

web-flow

Commit Message

add --trim option for generating smaller binaries (#55047)

This adds a command line option `--trim` that builds images where code
is only included if it is statically reachable from methods marked using
the new function `entrypoint`. Compile-time errors are given for call
sites that are too dynamic to allow trimming the call graph (however
there is an `unsafe` option if you want to try building anyway to see
what happens).

The PR has two other components. One is changes to Base that generally
allow more code to be compiled in this mode. These changes will either
be merged in separate PRs or moved to a separate part of the workflow
(where we will build a custom system image for this purpose). The branch
is set up this way to make it easy to check out and try the
functionality.

The other component is everything in the `juliac/` directory, which
implements a compiler driver script based on this new option, along with
some examples and tests. This will eventually become a package "app"
that depends on PackageCompiler and provides a CLI for all of this
stuff, so it will not be merged here. To try an example:

```
julia contrib/juliac.jl --output-exe hello --trim test/trimming/hello.jl
```

When stripped the resulting executable is currently about 900kb on my
machine.

Also includes a lot of work by @topolarity

---------

Co-authored-by: Gabriel Baraldi <baraldigabriel@gmail.com>
Co-authored-by: Tim Holy <tim.holy@gmail.com>
Co-authored-by: Cody Tapscott <topolarity@tapscott.me>

Run Details

2 of 10 new or added lines in 4 files covered. (20.0%)

1269 existing lines in 41 files now uncovered.

77375 of 89467 relevant lines covered (86.48%)

15490193.44 hits per line

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

88.57

/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::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,774,246✔
6	copysign(x::Float32, y::Float32) = copysign_float(x, y)	1,568,627✔
7	copysign(x::Float32, y::Real) = copysign(x, Float32(y))	51✔
8	copysign(x::Float64, y::Real) = copysign(x, Float64(y))	173,653✔
9
10	flipsign(x::Float64, y::Float64) = bitcast(Float64, xor_int(bitcast(UInt64, x), and_int(bitcast(UInt64, y), 0x8000000000000000)))	212,997✔
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,467✔
14
15	signbit(x::Float64) = signbit(bitcast(Int64, x))	33,174,964✔
16	signbit(x::Float32) = signbit(bitcast(Int32, x))	18,743,876✔
17	signbit(x::Float16) = signbit(bitcast(Int16, x))	733,957✔
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,124,482✔
32	maxintfloat(::Type{Float16}) = Float16(2048f0)	2,283,301✔
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,259,501✔
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,608,814✔
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,024,551✔
51	digits::Union{Nothing,Integer}=nothing, sigdigits::Union{Nothing,Integer}=nothing, base::Union{Nothing,Integer}=nothing)
52	if digits === nothing	46✔
53	if sigdigits === nothing	40✔
UNCOV 54	if base === nothing	×
55	# avoid recursive calls
56	throw(MethodError(round, (x,r)))	×
57	else
UNCOV 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)	42✔
63	_round_sigdigits(x, r, sigdigits, base === nothing ? 10 : base)	36✔
64	end
65	else
66	if sigdigits === nothing	6✔
67	isfinite(x) \|\| return float(x)	43,819✔
68	_round_digits(x, r, digits, base === nothing ? 10 : base)	43,817✔
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)
77	y = round(x * invstep, r) / invstep	43,766✔
78	if !isfinite(y)	43,766✔
79	return x	×
80	end
81	return y	43,766✔
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	12✔
88	if !isfinite(y)	12✔
89	return x	2✔
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	2✔
98	if !isfinite(y)	2✔
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
107	return y	2✔
108	end
109
110	function _round_digits(x, r::RoundingMode, digits::Integer, base)	43,780✔
111	fx = float(x)	12✔
112	if digits >= 0	43,780✔
113	invstep = oftype(fx, base)^digits	87,545✔
114	if isfinite(invstep)	43,778✔
115	return _round_invstep(fx, invstep, r)	43,766✔
116	else
117	invstepsqrt = oftype(fx, base)^oftype(fx, digits/2)	12✔
118	return _round_invstepsqrt(fx, invstepsqrt, r)	12✔
119	end
120	else
121	step = oftype(fx, base)^-digits	4✔
122	return _round_step(fx, step, r)	2✔
123	end
124	end
125
126	hidigit(x::Integer, base) = ndigits0z(x, base)	×
127	function hidigit(x::AbstractFloat, base)	36✔
128	iszero(x) && return 0	36✔
129	if base == 10	34✔
130	return 1 + floor(Int, log10(abs(x)))	32✔
131	elseif base == 2	2✔
132	return 1 + exponent(x)	1✔
133	else
134	return 1 + floor(Int, log(base, abs(x)))	1✔
135	end
136	end
UNCOV 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)	36✔
141	_round_digits(x, r, sigdigits-h, base)	36✔
142	end
143
144	# C-style round
145	function round(x::AbstractFloat, ::RoundingMode{:NearestTiesAway})	40,200✔
146	y = trunc(x)	40,200✔
147	ifelse(x==y,y,trunc(2*x-y))	40,200✔
148	end
149	# Java-style round
150	function round(x::T, ::RoundingMode{:NearestTiesUp}) where {T <: AbstractFloat}	40,200✔
151	copysign(floor((x + (T(0.25) - eps(T(0.5)))) + (T(0.25) + eps(T(0.5)))), x)	40,200✔
152	end
153
154	function Base.round(x::AbstractFloat, ::typeof(RoundFromZero))	40,200✔
155	signbit(x) ? round(x, RoundDown) : round(x, RoundUp)	40,200✔
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,236✔
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,850✔
224	x == y \|\|	386,009✔
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;	808✔
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	591✔
233	return x == y	467✔
234	else
235	# We need to take the difference `max` - `min` when comparing unsigned integers.
236	_x, _y = x < y ? (x, y) : (y, x)	136✔
237	return norm(_y - _x) <= max(atol, rtol*max(norm(_x), norm(_y)))	124✔
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...)	197✔
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...)	58✔
260
261	# default tolerance arguments
262	rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))	431,767✔
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}	19,182✔
265	rtol = max(rtoldefault(real(T)),rtoldefault(real(S)))	428,928✔
266	return atol > 0 ? zero(rtol) : rtol	413,652✔
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::Float32, b::Float32, c::Float32)::Float32	262,155✔
280	ab = Float64(a) * b	262,155✔
281	res = ab+c	262,155✔
282	reinterpret(UInt64, res)&0x1fff_ffff!=0x1000_0000 && return res	262,155✔
283	# yes error compensation is necessary. It sucks
284	reslo = abs(c)>abs(ab) ? ab-(res - c) : c-(res - ab)	4✔
285	res = iszero(reslo) ? res : (signbit(reslo) ? prevfloat(res) : nextfloat(res))	3✔
286	return res	2✔
287	end
288
289	""" 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"""
290	@inline function splitbits(x::Float64)
291	hi = reinterpret(Float64, reinterpret(UInt64, x) & 0xffff_ffff_f800_0000)	925,038✔
292	return hi, x-hi	925,038✔
293	end
294
295	function twomul(a::Float64, b::Float64)
296	ahi, alo = splitbits(a)	308,346✔
297	bhi, blo = splitbits(b)	308,346✔
298	abhi = a*b	308,346✔
299	blohi, blolo = splitbits(blo)	308,346✔
300	ablo = aloblohi - (((abhi - ahibhi) - alobhi) - ahiblo) + blolo*alo	308,346✔
301	return abhi, ablo	×
302	end
303
304	function fma_emulated(a::Float64, b::Float64,c::Float64)	262,456✔
305	abhi, ablo = @inline twomul(a,b)	262,456✔
306	if !isfinite(abhi+c) \|\| isless(abs(abhi), nextfloat(0x1p-969)) \|\| issubnormal(a) \|\| issubnormal(b)	491,458✔
307	aandbfinite = isfinite(a) && isfinite(b)	69,879✔
308	if !(isfinite(c) && aandbfinite)	69,879✔
309	return aandbfinite ? c : abhi+c	419✔
310	end
311	(iszero(a) \|\| iszero(b)) && return abhi+c	69,460✔
312	# The checks above satisfy exponent's nothrow precondition
313	bias = Math._exponent_finite_nonzero(a) + Math._exponent_finite_nonzero(b)	69,459✔
314	c_denorm = ldexp(c, -bias)	138,918✔
315	if isfinite(c_denorm)	69,459✔
316	# rescale a and b to [1,2), equivalent to ldexp(a, -exponent(a))
317	issubnormal(a) && (a *= 0x1p52)	45,890✔
318	issubnormal(b) && (b *= 0x1p52)	45,890✔
319	a = reinterpret(Float64, (reinterpret(UInt64, a) & ~Base.exponent_mask(Float64)) \| Base.exponent_one(Float64))	45,890✔
320	b = reinterpret(Float64, (reinterpret(UInt64, b) & ~Base.exponent_mask(Float64)) \| Base.exponent_one(Float64))	45,890✔
321	c = c_denorm	×
322	abhi, ablo = twomul(a,b)	45,890✔
323	# abhi <= 4 -> isfinite(r) (α)
324	r = abhi+c	45,890✔
325	# s ≈ 0 (β)
326	s = (abs(abhi) > abs(c)) ? (abhi-r+c+ablo) : (c-r+abhi+ablo)	58,666✔
327	# α ⩓ β -> isfinite(sumhi) (γ)
328	sumhi = r+s	45,890✔
329	# If result is subnormal, ldexp will cause double rounding because subnormals have fewer mantisa bits.
330	# As such, we need to check whether round to even would lead to double rounding and manually round sumhi to avoid it.
331	if issubnormal(ldexp(sumhi, bias))	91,780✔
332	sumlo = r-sumhi+s	15✔
333	# finite: See γ
334	# non-zero: If sumhi == ±0., then ldexp(sumhi, bias) == ±0,
335	# so we don't take this branch.
336	bits_lost = -bias-Math._exponent_finite_nonzero(sumhi)-1022	15✔
337	sumhiInt = reinterpret(UInt64, sumhi)	15✔
338	if (bits_lost != 1) ⊻ (sumhiInt&1 == 1)	15✔
339	sumhi = nextfloat(sumhi, cmp(sumlo,0))	10✔
340	end
341	end
342	return ldexp(sumhi, bias)	45,890✔
343	end
344	isinf(abhi) && signbit(c) == signbit(a*b) && return abhi	23,569✔
345	# fall through
346	end
347	r = abhi+c	216,146✔
348	s = (abs(abhi) > abs(c)) ? (abhi-r+c+ablo) : (c-r+abhi+ablo)	334,245✔
349	return r+s	216,146✔
350	end
351	fma_llvm(x::Float32, y::Float32, z::Float32) = fma_float(x, y, z)	2,448,598✔
352	fma_llvm(x::Float64, y::Float64, z::Float64) = fma_float(x, y, z)	408,659,024✔
353
354	# Disable LLVM's fma if it is incorrect, e.g. because LLVM falls back
355	# onto a broken system libm; if so, use a software emulated fma
356	@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✔
357	@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,659,024✔
358
359	function fma(a::Float16, b::Float16, c::Float16)	1✔
360	Float16(muladd(Float32(a), Float32(b), Float32(c))) #don't use fma if the hardware doesn't have it.	499,724✔
361	end
362
363	# This is necessary at least on 32-bit Intel Linux, since fma_llvm may
364	# have called glibc, and some broken glibc fma implementations don't
365	# properly restore the rounding mode
366	Rounding.setrounding_raw(Float32, Rounding.JL_FE_TONEAREST)
367	Rounding.setrounding_raw(Float64, Rounding.JL_FE_TONEAREST)

JuliaLang / julia / #37918

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