#37447

pending completion

Build # #37447

Build Type

push

local

Committed by

web-flow

Commit Message

`@inbounds` in `copyto!` for structured broadcasting (#48437)

* `@inbounds` in diagonal broadcasting

* inbounds copyto for bi/tridiag and triangular

* Move inbounds to broadcast getindex

---------

Co-authored-by: Daniel Karrasch <daniel.karrasch@posteo.de>

Run Details

12 of 12 new or added lines in 1 file covered. (100.0%)

71468 of 82456 relevant lines covered (86.67%)

32670681.69 hits per line

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

88.44

/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) = (x - trunc(x) == 0)

"""
    round([T,] x, [r::RoundingMode])
    round(x, [r::RoundingMode]; digits::Integer=0, base = 10)
    round(x, [r::RoundingMode]; sigdigits::Integer, base = 10)

Rounds the number `x`.

Without keyword arguments, `x` is rounded to an integer value, returning a value of type
`T`, or of the same type of `x` if no `T` is provided. An [`InexactError`](@ref) will be
thrown if the value is not representable by `T`, similar to [`convert`](@ref).

If the `digits` keyword argument is provided, it rounds to the specified number of digits
after the decimal place (or before if negative), in base `base`.

If the `sigdigits` keyword argument is provided, it rounds to the specified number of
significant digits, in base `base`.

The [`RoundingMode`](@ref) `r` controls the direction of the rounding; the default is
[`RoundNearest`](@ref), which rounds to the nearest integer, with ties (fractional values
of 0.5) being rounded to the nearest even integer. Note that `round` may give incorrect
results if the global rounding mode is changed (see [`rounding`](@ref)).

# Examples
```jldoctest
julia> round(1.7)
2.0

julia> round(Int, 1.7)
2

julia> round(1.5)
2.0

julia> round(2.5)
2.0

julia> round(pi; digits=2)
3.14

julia> round(pi; digits=3, base=2)
3.125

julia> round(123.456; sigdigits=2)
120.0

julia> round(357.913; sigdigits=4, base=2)
352.0
```

!!! note
    Rounding to specified digits in bases other than 2 can be inexact when
    operating on binary floating point numbers. For example, the [`Float64`](@ref)
    value represented by `1.15` is actually *less* than 1.15, yet will be
    rounded to 1.2. For example:

    ```jldoctest
    julia> x = 1.15
    1.15

    julia> big(1.15)
    1.149999999999999911182158029987476766109466552734375

    julia> x < 115//100
    true

    julia> round(x, digits=1)
    1.2
    ```

# Extensions

To extend `round` to new numeric types, it is typically sufficient to define `Base.round(x::NewType, r::RoundingMode)`.
"""
round(T::Type, x)

function round(::Type{T}, x::AbstractFloat, r::RoundingMode) where {T<:Integer}
    r != RoundToZero && (x = round(x,r))
    trunc(T, x)
end

# 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

trunc(x::Real; kwargs...) = round(x, RoundToZero; kwargs...)
floor(x::Real; kwargs...) = round(x, RoundDown; kwargs...)
ceil(x::Real; kwargs...)  = round(x, RoundUp; kwargs...)

round(x::Integer, r::RoundingMode) = x

# 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 || (isfinite(x) && isfinite(y) && norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))) || (nans && isnan(x) && isnan(y))
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)	2,463,512✔
6	copysign(x::Float32, y::Float32) = copysign_float(x, y)	1,547,610✔
7	copysign(x::Float32, y::Real) = copysign(x, Float32(y))	51✔
8	copysign(x::Float64, y::Real) = copysign(x, Float64(y))	96✔
9
10	flipsign(x::Float64, y::Float64) = bitcast(Float64, xor_int(bitcast(UInt64, x), and_int(bitcast(UInt64, y), 0x8000000000000000)))	115,314✔
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))	144✔
14
15	signbit(x::Float64) = signbit(bitcast(Int64, x))	16,638,916✔
16	signbit(x::Float32) = signbit(bitcast(Int32, x))	9,082,705✔
17	signbit(x::Float16) = signbit(bitcast(Int16, x))	8,814✔
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.	855,550✔
31	maxintfloat(::Type{Float32}) = Float32(16777216.)	3,125,333✔
32	maxintfloat(::Type{Float16}) = Float16(2048f0)	2,284,662✔
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,261,984✔
43	maxintfloat() = maxintfloat(Float64)	1✔
44
45	isinteger(x::AbstractFloat) = (x - trunc(x) == 0)	14,367✔
46
47	"""
48	round([T,] x, [r::RoundingMode])
49	round(x, [r::RoundingMode]; digits::Integer=0, base = 10)
50	round(x, [r::RoundingMode]; sigdigits::Integer, base = 10)
51
52	Rounds the number `x`.
53
54	Without keyword arguments, `x` is rounded to an integer value, returning a value of type
55	`T`, or of the same type of `x` if no `T` is provided. An [`InexactError`](@ref) will be
56	thrown if the value is not representable by `T`, similar to [`convert`](@ref).
57
58	If the `digits` keyword argument is provided, it rounds to the specified number of digits
59	after the decimal place (or before if negative), in base `base`.
60
61	If the `sigdigits` keyword argument is provided, it rounds to the specified number of
62	significant digits, in base `base`.
63
64	The [`RoundingMode`](@ref) `r` controls the direction of the rounding; the default is
65	[`RoundNearest`](@ref), which rounds to the nearest integer, with ties (fractional values
66	of 0.5) being rounded to the nearest even integer. Note that `round` may give incorrect
67	results if the global rounding mode is changed (see [`rounding`](@ref)).
68
69	# Examples
70	```jldoctest
71	julia> round(1.7)
72	2.0
73
74	julia> round(Int, 1.7)
75	2
76
77	julia> round(1.5)
78	2.0
79
80	julia> round(2.5)
81	2.0
82
83	julia> round(pi; digits=2)
84	3.14
85
86	julia> round(pi; digits=3, base=2)
87	3.125
88
89	julia> round(123.456; sigdigits=2)
90	120.0
91
92	julia> round(357.913; sigdigits=4, base=2)
93	352.0
94	```
95
96	!!! note
97	Rounding to specified digits in bases other than 2 can be inexact when
98	operating on binary floating point numbers. For example, the [`Float64`](@ref)
99	value represented by `1.15` is actually less than 1.15, yet will be
100	rounded to 1.2. For example:
101
102	```jldoctest
103	julia> x = 1.15
104	1.15
105
106	julia> big(1.15)
107	1.149999999999999911182158029987476766109466552734375
108
109	julia> x < 115//100
110	true
111
112	julia> round(x, digits=1)
113	1.2
114	```
115
116	# Extensions
117
118	To extend `round` to new numeric types, it is typically sufficient to define `Base.round(x::NewType, r::RoundingMode)`.
119	"""
120	round(T::Type, x)
121
122	function round(::Type{T}, x::AbstractFloat, r::RoundingMode) where {T<:Integer}	13,325✔
123	r != RoundToZero && (x = round(x,r))	13,325✔
124	trunc(T, x)	13,325✔
125	end
126
127	# NOTE: this relies on the current keyword dispatch behaviour (#9498).
128	function round(x::Real, r::RoundingMode=RoundNearest;	699,374✔
129	digits::Union{Nothing,Integer}=nothing, sigdigits::Union{Nothing,Integer}=nothing, base::Union{Nothing,Integer}=nothing)
130	if digits === nothing	46✔
131	if sigdigits === nothing	40✔
132	if base === nothing	×
133	# avoid recursive calls
134	throw(MethodError(round, (x,r)))	×
135	else
136	round(x,r)	×
137	# or throw(ArgumentError("`round` cannot use `base` argument without `digits` or `sigdigits` arguments."))
138	end
139	else
140	isfinite(x) \|\| return float(x)	42✔
141	_round_sigdigits(x, r, sigdigits, base === nothing ? 10 : base)	36✔
142	end
143	else
144	if sigdigits === nothing	6✔
145	isfinite(x) \|\| return float(x)	30,129✔
146	_round_digits(x, r, digits, base === nothing ? 10 : base)	30,127✔
147	else
148	throw(ArgumentError("`round` cannot use both `digits` and `sigdigits` arguments."))	×
149	end
150	end
151	end
152
153	trunc(x::Real; kwargs...) = round(x, RoundToZero; kwargs...)	197,690✔
154	floor(x::Real; kwargs...) = round(x, RoundDown; kwargs...)	86,729✔
155	ceil(x::Real; kwargs...) = round(x, RoundUp; kwargs...)	5,884✔
156
157	round(x::Integer, r::RoundingMode) = x	387✔
158
159	# round x to multiples of 1/invstep
160	function _round_invstep(x, invstep, r::RoundingMode)	11✔
161	y = round(x * invstep, r) / invstep	30,149✔
162	if !isfinite(y)	30,149✔
163	return x	×
164	end
165	return y	30,149✔
166	end
167
168	# round x to multiples of 1/(invstepsqrt^2)
169	# Using square root of step prevents overflowing
170	function _round_invstepsqrt(x, invstepsqrt, r::RoundingMode)	1✔
171	y = round((x * invstepsqrt) * invstepsqrt, r) / invstepsqrt / invstepsqrt	12✔
172	if !isfinite(y)	12✔
173	return x	2✔
174	end
175	return y	10✔
176	end
177
178	# round x to multiples of step
179	function _round_step(x, step, r::RoundingMode)	×
180	# TODO: use div with rounding mode
181	y = round(x / step, r) * step	2✔
182	if !isfinite(y)	2✔
183	if x > 0	×
184	return (r == RoundUp ? oftype(x, Inf) : zero(x))	×
185	elseif x < 0	×
186	return (r == RoundDown ? -oftype(x, Inf) : -zero(x))	×
187	else
188	return x	×
189	end
190	end
191	return y	2✔
192	end
193
194	function _round_digits(x, r::RoundingMode, digits::Integer, base)	30,163✔
195	fx = float(x)	12✔
196	if digits >= 0	30,163✔
197	invstep = oftype(fx, base)^digits	60,311✔
198	if isfinite(invstep)	30,161✔
199	return _round_invstep(fx, invstep, r)	30,149✔
200	else
201	invstepsqrt = oftype(fx, base)^oftype(fx, digits/2)	12✔
202	return _round_invstepsqrt(fx, invstepsqrt, r)	12✔
203	end
204	else
205	step = oftype(fx, base)^-digits	4✔
206	return _round_step(fx, step, r)	2✔
207	end
208	end
209
210	hidigit(x::Integer, base) = ndigits0z(x, base)	×
211	function hidigit(x::AbstractFloat, base)	36✔
212	iszero(x) && return 0	36✔
213	if base == 10	34✔
214	return 1 + floor(Int, log10(abs(x)))	32✔
215	elseif base == 2	2✔
216	return 1 + exponent(x)	1✔
217	else
218	return 1 + floor(Int, log(base, abs(x)))	1✔
219	end
220	end
221	hidigit(x::Real, base) = hidigit(float(x), base)	×
222
223	function _round_sigdigits(x, r::RoundingMode, sigdigits::Integer, base)	36✔
224	h = hidigit(x, base)	36✔
225	_round_digits(x, r, sigdigits-h, base)	36✔
226	end
227
228	# C-style round
229	function round(x::AbstractFloat, ::RoundingMode{:NearestTiesAway})	40,200✔
230	y = trunc(x)	40,200✔
231	ifelse(x==y,y,trunc(2*x-y))	40,200✔
232	end
233	# Java-style round
234	function round(x::T, ::RoundingMode{:NearestTiesUp}) where {T <: AbstractFloat}	40,200✔
235	copysign(floor((x + (T(0.25) - eps(T(0.5)))) + (T(0.25) + eps(T(0.5)))), x)	40,200✔
236	end
237
238	function Base.round(x::AbstractFloat, ::typeof(RoundFromZero))	40,200✔
239	signbit(x) ? round(x, RoundDown) : round(x, RoundUp)	40,200✔
240	end
241
242	# isapprox: approximate equality of numbers
243	"""
244	isapprox(x, y; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps, nans::Bool=false[, norm::Function])
245
246	Inexact equality comparison. Two numbers compare equal if their relative distance or their
247	absolute distance is within tolerance bounds: `isapprox` returns `true` if
248	`norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))`. The default `atol` (absolute tolerance) is zero and the
249	default `rtol` (relative tolerance) depends on the types of `x` and `y`. The keyword argument `nans` determines
250	whether or not NaN values are considered equal (defaults to false).
251
252	For real or complex floating-point values, if an `atol > 0` is not specified, `rtol` defaults to
253	the square root of [`eps`](@ref) of the type of `x` or `y`, whichever is bigger (least precise).
254	This corresponds to requiring equality of about half of the significant digits. Otherwise,
255	e.g. for integer arguments or if an `atol > 0` is supplied, `rtol` defaults to zero.
256
257	The `norm` keyword defaults to `abs` for numeric `(x,y)` and to `LinearAlgebra.norm` for
258	arrays (where an alternative `norm` choice is sometimes useful).
259	When `x` and `y` are arrays, if `norm(x-y)` is not finite (i.e. `±Inf`
260	or `NaN`), the comparison falls back to checking whether all elements of `x` and `y` are
261	approximately equal component-wise.
262
263	The binary operator `≈` is equivalent to `isapprox` with the default arguments, and `x ≉ y`
264	is equivalent to `!isapprox(x,y)`.
265
266	Note that `x ≈ 0` (i.e., comparing to zero with the default tolerances) is
267	equivalent to `x == 0` since the default `atol` is `0`. In such cases, you should either
268	supply an appropriate `atol` (or use `norm(x) ≤ atol`) or rearrange your code (e.g.
269	use `x ≈ y` rather than `x - y ≈ 0`). It is not possible to pick a nonzero `atol`
270	automatically because it depends on the overall scaling (the "units") of your problem:
271	for example, in `x - y ≈ 0`, `atol=1e-9` is an absurdly small tolerance if `x` is the
272	[radius of the Earth](https://en.wikipedia.org/wiki/Earth_radius) in meters,
273	but an absurdly large tolerance if `x` is the
274	[radius of a Hydrogen atom](https://en.wikipedia.org/wiki/Bohr_radius) in meters.
275
276	!!! compat "Julia 1.6"
277	Passing the `norm` keyword argument when comparing numeric (non-array) arguments
278	requires Julia 1.6 or later.
279
280	# Examples
281	```jldoctest
282	julia> isapprox(0.1, 0.15; atol=0.05)
283	true
284
285	julia> isapprox(0.1, 0.15; rtol=0.34)
286	true
287
288	julia> isapprox(0.1, 0.15; rtol=0.33)
289	false
290
291	julia> 0.1 + 1e-10 ≈ 0.1
292	true
293
294	julia> 1e-10 ≈ 0
295	false
296
297	julia> isapprox(1e-10, 0, atol=1e-8)
298	true
299
300	julia> isapprox([10.0^9, 1.0], [10.0^9, 2.0]) # using `norm`
301	true
302	```
303	"""
304	function isapprox(x::Number, y::Number;	726,370✔
305	atol::Real=0, rtol::Real=rtoldefault(x,y,atol),
306	nans::Bool=false, norm::Function=abs)
307	x == y \|\| (isfinite(x) && isfinite(y) && norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))) \|\| (nans && isnan(x) && isnan(y))	377,081✔
308	end
309
310	"""
311	isapprox(x; kwargs...) / ≈(x; kwargs...)
312
313	Create a function that compares its argument to `x` using `≈`, i.e. a function equivalent to `y -> y ≈ x`.
314
315	The keyword arguments supported here are the same as those in the 2-argument `isapprox`.
316
317	!!! compat "Julia 1.5"
318	This method requires Julia 1.5 or later.
319	"""
320	isapprox(y; kwargs...) = x -> isapprox(x, y; kwargs...)	127✔
321
322	const ≈ = isapprox
323	"""
324	x ≉ y
325
326	This is equivalent to `!isapprox(x,y)` (see [`isapprox`](@ref)).
327	"""
328	≉(args...; kws...) = !≈(args...; kws...)	58✔
329
330	# default tolerance arguments
331	rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))	782,545✔
332	rtoldefault(::Type{<:Real}) = 0	11,498✔
333	function rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}, atol::Real) where {T<:Number,S<:Number}	395,085✔
334	rtol = max(rtoldefault(real(T)),rtoldefault(real(S)))	395,085✔
335	return atol > 0 ? zero(rtol) : rtol	398,061✔
336	end
337
338	# fused multiply-add
339
340	"""
341	fma(x, y, z)
342
343	Computes `xy+z` without rounding the intermediate result `xy`. On some systems this is
344	significantly more expensive than `x*y+z`. `fma` is used to improve accuracy in certain
345	algorithms. See [`muladd`](@ref).
346	"""
347	function fma end
348	function fma_emulated(a::Float32, b::Float32, c::Float32)::Float32	262,155✔
349	ab = Float64(a) * b	262,155✔
350	res = ab+c	262,155✔
351	reinterpret(UInt64, res)&0x1fff_ffff!=0x1000_0000 && return res	262,155✔
352	# yes error compensation is necessary. It sucks
353	reslo = abs(c)>abs(ab) ? ab-(res - c) : c-(res - ab)	4✔
354	res = iszero(reslo) ? res : (signbit(reslo) ? prevfloat(res) : nextfloat(res))	3✔
355	return res	2✔
356	end
357
358	""" 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"""
359	@inline function splitbits(x::Float64)	×
360	hi = reinterpret(Float64, reinterpret(UInt64, x) & 0xffff_ffff_f800_0000)	922,482✔
361	return hi, x-hi	922,482✔
362	end
363
364	function twomul(a::Float64, b::Float64)	×
365	ahi, alo = splitbits(a)	307,494✔
366	bhi, blo = splitbits(b)	307,494✔
367	abhi = a*b	307,494✔
368	blohi, blolo = splitbits(blo)	307,494✔
369	ablo = aloblohi - (((abhi - ahibhi) - alobhi) - ahiblo) + blolo*alo	307,494✔
370	return abhi, ablo	×
371	end
372
373	function fma_emulated(a::Float64, b::Float64,c::Float64)	262,160✔
374	abhi, ablo = @inline twomul(a,b)	262,160✔
375	if !isfinite(abhi+c) \|\| isless(abs(abhi), nextfloat(0x1p-969)) \|\| issubnormal(a) \|\| issubnormal(b)	491,383✔
376	aandbfinite = isfinite(a) && isfinite(b)	69,263✔
377	if !(isfinite(c) && aandbfinite)	69,263✔
378	return aandbfinite ? c : abhi+c	404✔
379	end
380	(iszero(a) \|\| iszero(b)) && return abhi+c	68,859✔
381	# The checks above satisfy exponent's nothrow precondition
382	bias = Math._exponent_finite_nonzero(a) + Math._exponent_finite_nonzero(b)	68,858✔
383	c_denorm = ldexp(c, -bias)	137,716✔
384	if isfinite(c_denorm)	68,858✔
385	# rescale a and b to [1,2), equivalent to ldexp(a, -exponent(a))
386	issubnormal(a) && (a *= 0x1p52)	45,334✔
387	issubnormal(b) && (b *= 0x1p52)	45,334✔
388	a = reinterpret(Float64, (reinterpret(UInt64, a) & ~Base.exponent_mask(Float64)) \| Base.exponent_one(Float64))	45,334✔
389	b = reinterpret(Float64, (reinterpret(UInt64, b) & ~Base.exponent_mask(Float64)) \| Base.exponent_one(Float64))	45,334✔
390	c = c_denorm	×
391	abhi, ablo = twomul(a,b)	45,334✔
392	# abhi <= 4 -> isfinite(r) (α)
393	r = abhi+c	45,334✔
394	# s ≈ 0 (β)
395	s = (abs(abhi) > abs(c)) ? (abhi-r+c+ablo) : (c-r+abhi+ablo)	58,055✔
396	# α ⩓ β -> isfinite(sumhi) (γ)
397	sumhi = r+s	45,334✔
398	# If result is subnormal, ldexp will cause double rounding because subnormals have fewer mantisa bits.
399	# As such, we need to check whether round to even would lead to double rounding and manually round sumhi to avoid it.
400	if issubnormal(ldexp(sumhi, bias))	90,668✔
401	sumlo = r-sumhi+s	16✔
402	# finite: See γ
403	# non-zero: If sumhi == ±0., then ldexp(sumhi, bias) == ±0,
404	# so we don't take this branch.
405	bits_lost = -bias-Math._exponent_finite_nonzero(sumhi)-1022	16✔
406	sumhiInt = reinterpret(UInt64, sumhi)	16✔
407	if (bits_lost != 1) ⊻ (sumhiInt&1 == 1)	16✔
408	sumhi = nextfloat(sumhi, cmp(sumlo,0))	12✔
409	end
410	end
411	return ldexp(sumhi, bias)	45,334✔
412	end
413	isinf(abhi) && signbit(c) == signbit(a*b) && return abhi	23,524✔
414	# fall through
415	end
416	r = abhi+c	216,421✔
417	s = (abs(abhi) > abs(c)) ? (abhi-r+c+ablo) : (c-r+abhi+ablo)	334,563✔
418	return r+s	216,421✔
419	end
420	fma_llvm(x::Float32, y::Float32, z::Float32) = fma_float(x, y, z)	2,448,598✔
421	fma_llvm(x::Float64, y::Float64, z::Float64) = fma_float(x, y, z)	404,478,542✔
422
423	# Disable LLVM's fma if it is incorrect, e.g. because LLVM falls back
424	# onto a broken system libm; if so, use a software emulated fma
425	@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✔
426	@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)	404,478,542✔
427
428	function fma(a::Float16, b::Float16, c::Float16)	999,728✔
429	Float16(muladd(Float32(a), Float32(b), Float32(c))) #don't use fma if the hardware doesn't have it.	999,728✔
430	end
431
432	# This is necessary at least on 32-bit Intel Linux, since fma_llvm may
433	# have called glibc, and some broken glibc fma implementations don't
434	# properly restore the rounding mode
435	Rounding.setrounding_raw(Float32, Rounding.JL_FE_TONEAREST)
436	Rounding.setrounding_raw(Float64, Rounding.JL_FE_TONEAREST)

JuliaLang / julia / #37447

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