Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Arithmetic cleanup and additions #164

Merged
merged 6 commits into from
Aug 21, 2023
Merged

Arithmetic cleanup and additions #164

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
a6b3ba9
Split arithmetic.jl into three files
Joel-Dahne Aug 7, 2023
07cf9a9
Minor cleanup of elementary.jl
Joel-Dahne Aug 7, 2023
a06d9f8
Cleanup of arithmethic.jl and a few new methods
Joel-Dahne Aug 7, 2023
db6745f
Cleanup of minmax.jl
Joel-Dahne Aug 8, 2023
9a10de8
Factor out multi-argument methods into separate file
Joel-Dahne Aug 19, 2023
dc93618
Add optimized arithmethic with rationals and integers
Joel-Dahne Aug 20, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions src/Arblib.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,8 @@ include("predicates.jl")
include("show.jl")
include("promotion.jl")
include("arithmetic.jl")
include("elementary.jl")
include("minmax.jl")
include("rand.jl")
include("float.jl")
include("interval.jl")
Expand Down
329 changes: 119 additions & 210 deletions src/arithmetic.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,55 +2,79 @@
const _BitSigned = (Int == Int64 ? Base.BitSigned64 : Base.BitSigned32)
const _BitUnsigned = (Int == Int64 ? Base.BitUnsigned64 : Base.BitUnsigned32)

### Mag
### +, -, *, /
for (jf, af) in [(:+, :add!), (:-, :sub!), (:*, :mul!), (:/, :div!)]
@eval Base.$jf(x::MagOrRef, y::MagOrRef) = $af(zero(x), x, y)
end
Base.:+(x::MagOrRef, y::Integer) = add!(zero(x), x, convert(UInt, y))
Base.:+(x::Integer, y::MagOrRef) = add!(zero(y), y, convert(UInt, x))
Base.:*(x::MagOrRef, y::Integer) = mul!(zero(x), x, convert(UInt, y))
Base.:*(x::Integer, y::MagOrRef) = mul!(zero(y), y, convert(UInt, x))
Base.:/(x::MagOrRef, y::Integer) = div!(zero(x), x, convert(UInt, y))

Base.:(^)(x::MagOrRef, e::Integer) = pow!(zero(x), x, convert(UInt, e))
rsqrt(x::MagOrRef) = rsqrt!(zero(x), x)
Base.hypot(x::MagOrRef, y::MagOrRef) = hypot!(zero(x), x, y)
root(x::MagOrRef, n::Integer) = root!(zero(x), x, convert(UInt, n))
neglog(x::MagOrRef) = neg_log!(zero(x), x)
expinv(x::MagOrRef) = expinv!(zero(x), x)
for f in [:inv, :sqrt, :log, :log1p, :exp, :expm1, :atan, :cosh, :sinh]
@eval Base.$f(x::MagOrRef) = $(Symbol(f, :!))(zero(x), x)
end

Base.min(x::MagOrRef, y::MagOrRef) = Arblib.min!(zero(x), x, y)
Base.max(x::MagOrRef, y::MagOrRef) = Arblib.max!(zero(x), x, y)
Base.minmax(x::MagOrRef, y::MagOrRef) = (min(x, y), max(x, y))

### Arf
function Base.sign(x::ArfOrRef)
isnan(x) && return Arf(NaN) # Follow Julia and return NaN
return Arf(sgn(x))
end

Base.abs(x::ArfOrRef) = abs!(zero(x), x)
Base.:(-)(x::ArfOrRef) = neg!(zero(x), x)
for (jf, af) in [(:+, :add!), (:-, :sub!), (:*, :mul!), (:/, :div!)]
@eval function Base.$jf(x::ArfOrRef, y::Union{ArfOrRef,_BitInteger})
z = Arf(prec = _precision(x, y))
$af(z, x, y)
return z
end

@eval Base.$jf(x::ArbOrRef, y::Union{ArbOrRef,ArfOrRef,_BitInteger}) =
$af(Arb(prec = _precision(x, y)), x, y)

@eval Base.$jf(x::AcbOrRef, y::Union{AcbOrRef,ArbOrRef,_BitInteger}) =
$af(Acb(prec = _precision(x, y)), x, y)

# Avoid one allocation for operations on irrationals
@eval function Base.$jf(x::Union{ArbOrRef,AcbOrRef}, y::Irrational)
z = zero(x)
z[] = y
return $af(z, x, z)

Check warning on line 25 in src/arithmetic.jl

View check run for this annotation

Codecov / codecov/patch

src/arithmetic.jl#L22-L25 

Added lines #L22 - L25 were not covered by tests
end

if jf == :+ || jf == :*
# Addition and multiplication is commutative
@eval Base.$jf(x::_BitInteger, y::ArfOrRef) = $jf(y, x)

@eval Base.$jf(x::Union{ArfOrRef,_BitInteger,Irrational}, y::ArbOrRef) = $jf(y, x)

@eval Base.$jf(x::Union{ArbOrRef,_BitInteger,Irrational}, y::AcbOrRef) = $jf(y, x)

# Define more efficient multi argument versions, similar to BigFloat
# TODO: Now precision is determined by first two arguments
for T in (MagOrRef, ArfOrRef, ArbOrRef, AcbOrRef)
@eval function Base.$jf(a::$T, b::$T, c::$T)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe redefine _precision to handle multiple arguments:

_precision(a,b,args...) = 
    max(_precision(a), _precision(b, args...))

and then

function Base.$jf(a::$T, b::$T, args...)
    z = $T(prec=_precision(a, b, args...))
    z = $af(z, a, b)
    for x in args
        z = $af(z, z, x)
    end
    return z
end

?
we could also abuse this pattern and define

$af(z, a, b, args...) = (z = $af(z, z, a); $af(z, b, args...))

but that changes the semantics of $af, so better not do this directly (function wrapper in the middle would be ok?)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Leaving aside the precision question for now. I got interested in the different ways to write a method handling any numbers of arguments. The results were somewhat surprising and differs a lot between Julia 1.9.2 and 1.10-beta. The methods I tried were

function mul_test_1(a::Arb, b::Arb, c::Arb, d::Arb, e::Arb)
    z = a * b
    mul!(z, z, c)
    mul!(z, z, d)
    mul!(z, z, e)
    return z
end

function mul_test_2(a::Arb, b::Arb, c::Arb, rest::Vararg{Arb})
    z = a * b
    mul!(z, z, c)
    for x in rest
        mul!(z, z, x)
    end
    return z
end

mul_test_3(a::Arb, b::Arb) = a * b
function mul_test_3(a::Arb, b::Arb, c::Arb, rest::Vararg{Arb})
    z = mul_test_3(b, c, rest...)
    mul!(z, a, z)
    return z
end

The tests I ran are

a = Arb(1)
b1 = @benchmarkable Arblib.mul_test_1($a, $a, $a, $a, $a)
b2 = @benchmarkable Arblib.mul_test_2($a, $a, $a, $a, $a)
b3 = @benchmarkable Arblib.mul_test_3($a, $a, $a, $a, $a)

In Julia 1.9.2 I get the following results

julia> run(b1, evals=1000)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):   68.075 ns … 54.788 μs  ┊ GC (min … max):  0.00% … 58.35%
 Time  (median):      70.518 ns              ┊ GC (median):     0.00%
 Time  (mean ± σ):   147.215 ns ±  1.955 μs  ┊ GC (mean ± σ):  29.33% ±  2.20%

  ▁▅▇█▇▅▄▃▃▃▂▃▄▄▃▂▁                  ▁▂▃▃▃▂▁▁                  ▂
  ██████████████████▇█▇▇▇▆▆▆▆▅▆▅▅▄▄▄▆██████████▇▇▇▆▆▆▆▅▅▅▆▆▅▄▅ █
  68.1 ns       Histogram: log(frequency) by time       102 ns <

 Memory estimate: 64 bytes, allocs estimate: 1.

julia> run(b2, evals=1000)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):   97.290 ns … 15.608 μs  ┊ GC (min … max):  0.00% … 61.43%
 Time  (median):     100.489 ns              ┊ GC (median):     0.00%
 Time  (mean ± σ):   184.264 ns ±  1.084 μs  ┊ GC (mean ± σ):  27.50% ±  4.62%

  ▃▇█▆▄▄▃▄▅▄▃▂▂▂▂▂▁▁                                           ▂
  ███████████████████▇█▇▇▇▆▆▆▆▄▄▄▅▅▄▄▄▁▁▃▄▄▃▁▃▁▁▃▄▆▆▅▆▇▇▇▆▃▆▆▆ █
  97.3 ns       Histogram: log(frequency) by time       168 ns <

 Memory estimate: 256 bytes, allocs estimate: 7.

julia> run(b3, evals=1000)
BenchmarkTools.Trial: 8418 samples with 1000 evaluations.
 Range (min … max):  465.435 ns …   8.043 μs  ┊ GC (min … max):  0.00% … 67.40%
 Time  (median):     486.168 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   592.378 ns ± 847.952 ns  ┊ GC (mean ± σ):  12.32% ±  7.95%

  █                                                             ▁
  ██▆▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▆▇ █
  465 ns        Histogram: log(frequency) by time       7.54 μs <

 Memory estimate: 656 bytes, allocs estimate: 18.

So version 1 is the fastest and the only one that doesn't do extra allocations. The third one allocates a lot and is much slower!

With Julia 1.10-beta I instead get

julia> run(b1, evals=1000)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):   76.832 ns … 378.336 μs  ┊ GC (min … max):  0.00% … 56.04%
 Time  (median):      98.380 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   136.132 ns ±   3.782 μs  ┊ GC (mean ± σ):  15.58% ±  0.56%

                                    ▂▅██▃                        
  ▃▄▅▄▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▁▁▂▁▂▂▁▂▂▄▆██████▆▅▅▅▄▄▃▃▃▃▃▃▃▃▃▂▂▂▂▂▂ ▃
  76.8 ns          Histogram: frequency by time          112 ns <

 Memory estimate: 64 bytes, allocs estimate: 1.

julia> run(b2, evals=1000)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  108.310 ns … 324.235 μs  ┊ GC (min … max):  0.00% … 67.75%
 Time  (median):     194.815 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   287.768 ns ±   5.343 μs  ┊ GC (mean ± σ):  24.01% ±  1.35%

   ▄                                        ▁▆█▆▆▆▇▄             
  ▄█▃▃▄▃▂▂▂▂▂▂▂▂▂▂▂▂▂▁▁▁▂▁▂▁▁▁▂▂▃▄▃▄▅▄▃▄▄▄▄▅████████▇▇▅▄▃▃▂▂▂▂▂ ▄
  108 ns           Histogram: frequency by time          225 ns <

 Memory estimate: 256 bytes, allocs estimate: 7.

julia> run(b3, evals=1000)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):   78.538 ns … 379.233 μs  ┊ GC (min … max):  0.00% … 55.86%
 Time  (median):     100.804 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   139.085 ns ±   3.791 μs  ┊ GC (mean ± σ):  15.23% ±  0.56%

                                   ▃▆█▄                          
  ▂▄▄▄▂▂▂▂▂▂▂▂▂▂▂▁▂▁▁▂▂▂▁▁▂▁▂▂▂▂▃▄██████▅▅▄▄▄▃▃▃▂▃▃▃▄▄▃▂▂▂▂▂▂▂▂ ▃
  78.5 ns          Histogram: frequency by time          117 ns <

 Memory estimate: 64 bytes, allocs estimate: 1.

So now version 1 and 3 are more or less the same, in particular version 3 is not doing any extra allocations anymore. The second version is still slightly slower. Note also that version 1 is slightly faster in 1.9.2 than in 1.10-beta, I'm not sure why this is the case (can you reproduce it?).

My conclusion from this is that it is maybe not the best idea to write a generic multi argument version, but stick with the manual implementation of 3, 4 and 5 arguments. This is also the same as BigFloat does.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is the issue of inlining; on julia-1.9.2 with

mul_test_3(a::Arb, b::Arb) = a * b
function mul_test_3(a::Arb, b::Arb, c::Arb, rest::Vararg{Arb})
    z = mul_test_3(b, c, rest...)
    Arblib.mul!(z, a, z)
    return z
end

mul_test_4(a::Arb, b::Arb) = a * b
@inline function mul_test_4(a::Arb, b::Arb, c::Arb, rest::Vararg{Arb})
    z = mul_test_4(b, c, rest...)
    Arblib.mul!(z, a, z)
    return z
end

a = Arb(1)
b3 = @benchmarkable mul_test_3($a, $a, $a, $a, $a)
b4 = @benchmarkable mul_test_4($a, $a, $a, $a, $a)

run(b3, evals=1000)
run(b4, evals=1000)

i get

julia> run(b3, evals=1000)
BenchmarkTools.Trial: 5755 samples with 1000 evaluations.
 Range (min … max):  668.101 ns … 16.772 μs  ┊ GC (min … max):  0.00% … 74.94%
 Time  (median):     695.130 ns              ┊ GC (median):     0.00%
 Time  (mean ± σ):   852.314 ns ±  1.251 μs  ┊ GC (mean ± σ):  11.42% ±  7.25%

  █                                                            ▁
  █▇█▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▆▆ █
  668 ns        Histogram: log(frequency) by time      10.6 μs <

 Memory estimate: 656 bytes, allocs estimate: 18.

julia> run(b4, evals=1000)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):   88.698 ns … 87.804 μs  ┊ GC (min … max):  0.00% … 59.67%
 Time  (median):      91.143 ns              ┊ GC (median):     0.00%
 Time  (mean ± σ):   202.562 ns ±  2.812 μs  ┊ GC (mean ± σ):  28.17% ±  2.02%

  ▃██▆▄▂                                      ▁▁▂▃▂▂▁▁         ▂
  ███████▆▅▆▆▄▆██▇█▇▇▇▅▆▆▆▆▆▅▅▇▇▆▅▅▄▅▆▅▄▄▅▆▆▇██████████▇▆▆▇▆▇▇ █
  88.7 ns       Histogram: log(frequency) by time       150 ns <

 Memory estimate: 64 bytes, allocs estimate: 1.

z = $jf(a, b)
$af(z, z, c)
return z
end

@eval function Base.$jf(a::$T, b::$T, c::$T, d::$T)
z = $jf(a, b)
$af(z, z, c)
$af(z, z, d)
return z
end

@eval function Base.$jf(a::$T, b::$T, c::$T, d::$T, e::$T)
z = $jf(a, b)
$af(z, z, c)
$af(z, z, d)
$af(z, z, e)
return z
end
end
else
@eval function Base.$jf(x::Irrational, y::Union{ArbOrRef,AcbOrRef})
z = zero(y)
z[] = x
return $af(z, z, y)

Check warning on line 64 in src/arithmetic.jl

View check run for this annotation

Codecov / codecov/patch

src/arithmetic.jl#L61-L64 

Added lines #L61 - L64 were not covered by tests
end
end
end
function Base.:+(x::_BitInteger, y::ArfOrRef)
z = zero(y)
add!(z, y, x)
return z
end
function Base.:*(x::_BitInteger, y::ArfOrRef)
z = zero(y)
mul!(z, y, x)
return z
end

Base.:-(x::Union{ArfOrRef,ArbOrRef,AcbOrRef}) = neg!(zero(x), x)
Base.inv(x::Union{MagOrRef,ArbOrRef,AcbOrRef}) = inv!(zero(x), x)

Base.:+(x::MagOrRef, y::Integer) = add!(zero(x), x, convert(UInt, y))
Base.:+(x::Integer, y::MagOrRef) = y + x
Base.:*(x::MagOrRef, y::Integer) = mul!(zero(x), x, convert(UInt, y))
Base.:*(x::Integer, y::MagOrRef) = y * x
Base.:/(x::MagOrRef, y::Integer) = div!(zero(x), x, convert(UInt, y))

function Base.:/(x::_BitUnsigned, y::ArfOrRef)
z = zero(y)
ui_div!(z, x, y)
Expand All @@ -62,193 +86,78 @@
return z
end

function Base.sqrt(x::ArfOrRef)
y = zero(x)
sqrt!(y, x)
return y
end
function rsqrt(x::ArfOrRef)
y = zero(x)
rsqrt!(y, x)
return y
Base.:(/)(x::_BitUnsigned, y::ArbOrRef) = ui_div!(zero(y), x, y)

function Base.:*(x::AcbOrRef, y::Complex{Bool})
if real(y)
if imag(y)
z = mul_onei!(zero(x), x)
return add!(z, x, z)
else
return Acb(x)
end
end
imag(y) && return mul_onei!(zero(x), x)
return zero(x)
end
function root(x::ArfOrRef, k::Integer)
y = zero(x)
root!(y, x, convert(UInt, k))
return y
Base.:*(x::Complex{Bool}, y::AcbOrRef) = y * x

### fma and muladd
function Base.fma(x::ArfOrRef, y::ArfOrRef, z::ArfOrRef)
res = zero(x)
fma!(res, x, y, z)
return res
end
Base.fma(x::ArbOrRef, y::ArbOrRef, z::ArbOrRef) = fma!(zero(x), x, y, z)
Base.fma(x::ArbOrRef, y::Union{_BitInteger,ArfOrRef}, z::ArbOrRef) = fma!(zero(x), x, y, z)
Base.fma(x::Union{_BitInteger,ArfOrRef}, y::ArbOrRef, z::ArbOrRef) = fma(y, x, z)

Base.min(x::ArfOrRef, y::ArfOrRef) = Arblib.min!(zero(x), x, y)
Base.max(x::ArfOrRef, y::ArfOrRef) = Arblib.max!(zero(x), x, y)
Base.minmax(x::ArfOrRef, y::ArfOrRef) = (min(x, y), max(x, y))
Base.muladd(x::ArfOrRef, y::ArfOrRef, z::ArfOrRef) = fma(x, y, z)
Base.muladd(x::ArbOrRef, y::ArbOrRef, z::ArbOrRef) = fma(x, y, z)
Base.muladd(x::ArbOrRef, y::Union{_BitInteger,ArfOrRef}, z::ArbOrRef) = fma(x, y, z)
Base.muladd(x::Union{_BitInteger,ArfOrRef}, y::ArbOrRef, z::ArbOrRef) = fma(x, y, z)

### Arb and Acb
for (jf, af) in [(:+, :add!), (:-, :sub!), (:*, :mul!), (:/, :div!)]
@eval Base.$jf(x::ArbOrRef, y::Union{ArbOrRef,ArfOrRef,_BitInteger}) =
$af(Arb(prec = _precision(x, y)), x, y)
@eval Base.$jf(x::AcbOrRef, y::Union{AcbOrRef,ArbOrRef,_BitInteger}) =
$af(Acb(prec = _precision(x, y)), x, y)
if jf == :(+) || jf == :(*)
@eval Base.$jf(x::Union{ArfOrRef,_BitInteger}, y::ArbOrRef) =
$af(Arb(prec = _precision(x, y)), y, x)
@eval Base.$jf(x::Union{ArbOrRef,_BitInteger}, y::AcbOrRef) =
$af(Acb(prec = _precision(x, y)), y, x)
end
end
### signbit, sign and abs
Base.signbit(x::MagOrRef) = false
Base.signbit(x::ArfOrRef) = !isnan(x) && sgn(x) < 0
Base.signbit(x::ArbOrRef) = isnegative(x)

Base.:(-)(x::Union{ArbOrRef,AcbOrRef}) = neg!(zero(x), x)
Base.abs(x::ArbOrRef) = abs!(zero(x), x)
Base.:(/)(x::_BitUnsigned, y::ArbOrRef) = ui_div!(zero(y), x, y)
# For Arf sign of NaN is undefined, we follow Julia and return NaN
Base.sign(x::MagOrRef) = iszero(x) ? zero(x) : one(x)
Base.sign(x::ArfOrRef) = isnan(x) ? nan!(zero(x)) : Arf(sgn(x), prec = _precision(x))
Base.sign(x::Union{ArbOrRef,AcbRef}) = sgn!(zero(x), x)

Base.abs(x::MagOrRef) = copy(x)
Base.abs(x::Union{ArfOrRef,ArbOrRef}) = abs!(zero(x), x)
Base.abs(z::AcbOrRef) = abs!(Arb(prec = _precision(z)), z)

Base.:(^)(x::ArbOrRef, y::ArbOrRef) = pow!(Arb(prec = _precision(x, y)), x, y)
function Base.:(^)(x::ArbOrRef, y::_BitInteger)
### ^

Base.:^(x::MagOrRef, e::Integer) = pow!(zero(x), x, convert(UInt, e))
Base.:^(x::ArbOrRef, y::ArbOrRef) = pow!(Arb(prec = _precision(x, y)), x, y)
function Base.:^(x::ArbOrRef, y::_BitInteger)
z = zero(x)
x, n = (y >= 0 ? (x, y) : (inv!(z, x), -y))
return pow!(z, x, convert(UInt, n))
end
Base.:(^)(x::AcbOrRef, y::Union{AcbOrRef,ArbOrRef,_BitInteger}) =
Base.:^(x::AcbOrRef, y::Union{AcbOrRef,ArbOrRef,_BitInteger}) =
pow!(Acb(prec = _precision(x, y)), x, y)

# We define the same special cases as Arb does, this avoids some
# overhead
sqr(x::Union{ArbOrRef,AcbOrRef}) = sqr!(zero(x), x)
cube(x::AcbOrRef) = cube!(zero(x), x)

# We define the same special cases as Arb does, this avoids some overhead
Base.literal_pow(::typeof(^), x::AcbOrRef, ::Val{-3}) = (y = inv(x); cube!(y, y))
Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{-2}) =
(y = inv(x); sqr!(y, y))
#Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{-1}) - implemented in Base.intfuncs.jl
Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{0}) = one(x)
Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{1}) = copy(x)
Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{2}) = sqr(x)
Base.literal_pow(::typeof(^), x::AcbOrRef, ::Val{3}) = cube(x)

Base.hypot(x::ArbOrRef, y::ArbOrRef) = hypot!(Arb(prec = _precision(x, y)), x, y)

root(x::Union{ArbOrRef,AcbOrRef}, k::Integer) = root!(zero(x), x, convert(UInt, k))

# Unary methods in Base
for f in [
:inv,
:sqrt,
:log,
:log1p,
:exp,
:expm1,
:sin,
:cos,
:tan,
:cot,
:sec,
:csc,
:atan,
:asin,
:acos,
:sinh,
:cosh,
:tanh,
:coth,
:sech,
:csch,
:atanh,
:asinh,
:acosh,
]
@eval Base.$f(x::Union{ArbOrRef,AcbOrRef}) = $(Symbol(f, :!))(zero(x), x)
end

sqrtpos(x::ArbOrRef) = sqrtpos!(zero(x), x)
sqrt1pm1(x::ArbOrRef) = sqrt1pm1!(zero(x), x)
rsqrt(x::Union{ArbOrRef,AcbOrRef}) = rsqrt!(zero(x), x)
sqr(x::Union{ArbOrRef,AcbOrRef}) = sqr!(zero(x), x)

Base.sinpi(x::Union{ArbOrRef,AcbOrRef}) = sin_pi!(zero(x), x)
Base.cospi(x::Union{ArbOrRef,AcbOrRef}) = cos_pi!(zero(x), x)
tanpi(x::Union{ArbOrRef,AcbOrRef}) = tan_pi!(zero(x), x)
cotpi(x::Union{ArbOrRef,AcbOrRef}) = cot_pi!(zero(x), x)
cscpi(x::Union{ArbOrRef,AcbOrRef}) = csc_pi!(zero(x), x)
# Julias definition of sinc is equivalent to Arbs definition of sincpi
Base.sinc(x::Union{ArbOrRef,AcbOrRef}) = sinc_pi!(zero(x), x)
Base.atan(y::ArbOrRef, x::ArbOrRef) = atan2!(Arb(prec = _precision(y, x)), y, x)

function Base.sincos(x::Union{ArbOrRef,AcbOrRef})
s, c = zero(x), zero(x)
sin_cos!(s, c, x)
return (s, c)
end
function Base.sincospi(x::Union{ArbOrRef,AcbOrRef})
s, c = zero(x), zero(x)
sin_cos_pi!(s, c, x)
return (s, c)
end
function sinhcosh(x::Union{ArbOrRef,AcbOrRef})
s, c = zero(x), zero(x)
sinh_cosh!(s, c, x)
return (s, c)
end

Base.min(x::ArbOrRef, y::ArbOrRef) = Arblib.min!(zero(x), x, y)
Base.max(x::ArbOrRef, y::ArbOrRef) = Arblib.max!(zero(x), x, y)
Base.minmax(x::ArbOrRef, y::ArbOrRef) = (min(x, y), max(x, y))

### Acb
function Base.:(*)(x::AcbOrRef, y::Complex{Bool})
if real(y)
if imag(y)
z = mul_onei!(zero(x), x)
return add!(z, x, z)
else
return Acb(x)
end
end
imag(y) && return mul_onei!(zero(x), x)
return zero(x)
end
Base.:(*)(x::Complex{Bool}, y::AcbOrRef) = y * x

Base.real(z::AcbLike; prec = _precision(z)) = get_real!(Arb(; prec), z)
Base.imag(z::AcbLike; prec = _precision(z)) = get_imag!(Arb(; prec), z)
Base.conj(z::AcbLike) = conj!(Acb(prec = _precision(z)), z)
Base.abs(z::AcbLike) = abs!(Arb(prec = _precision(z)), z)

### minimum and maximum
# The default implemented in Julia have several issues for Arb types.
# See https://github.com/JuliaLang/julia/issues/45932.
# Note that it works fine for Mag and Arf.

# Before 1.8.0 there is no way to fix the implementation in Base.
# Instead we define a new method. Note that this doesn't fully fix the
# problem, there is no way to dispatch on for example
# minimum(x -> Arb(x), [1, 2, 3, 4]) or an array with only some Arb.
if VERSION < v"1.8.0-rc3"
function Base.minimum(A::AbstractArray{<:ArbOrRef})
isempty(A) &&
throw(ArgumentError("reducing over an empty collection is not allowed"))
res = copy(first(A))
for x in Iterators.drop(A, 1)
Arblib.min!(res, res, x)
end
return res
end

function Base.maximum(A::AbstractArray{<:ArbOrRef})
isempty(A) &&
throw(ArgumentError("reducing over an empty collection is not allowed"))
res = copy(first(A))
for x in Iterators.drop(A, 1)
Arblib.max!(res, res, x)
end
return res
end
end
### real, imag, conj

# Since 1.8.0 it is possible to fix the Base implementation by
# overloading some internal methods. This also works before 1.8.0 but
# doesn't solve the full problem.

# The default implementation in Base is not correct for Arb
Base._fast(::typeof(min), x::Arb, y::Arb) = min(x, y)
Base._fast(::typeof(min), x::Arb, y) = min(x, y)
Base._fast(::typeof(min), x, y::Arb) = min(x, y)
Base._fast(::typeof(max), x::Arb, y::Arb) = max(x, y)
Base._fast(::typeof(max), x::Arb, y) = max(x, y)
Base._fast(::typeof(max), x, y::Arb) = max(x, y)

# Mag, Arf and Arb don't have signed zeros
Base.isbadzero(::typeof(min), x::Union{Mag,Arf,Arb}) = false
Base.isbadzero(::typeof(max), x::Union{Mag,Arf,Arb}) = false
Base.real(z::AcbOrRef) = get_real!(Arb(prec = _precision(z)), z)
Base.imag(z::AcbOrRef) = get_imag!(Arb(prec = _precision(z)), z)
Base.conj(z::AcbOrRef) = conj!(zero(z), z)
Loading
Loading