Merge pull request #114 from kalmarek/enh/improve-arithmetic-support

Enh/improve arithmetic support

README.md

@@ -140,9 +140,72 @@ x = Arblib.set!(Arb(), (0, π))
 y = Arb((0, π))
 ```
 
+## Pitfalls when interacting with the Julia ecosystem
+Arb is made for rigorous numerics and any functions which do not
+produce rigorous results are clearly marked as such. This is not the
+case with Julia in general and you therefore have to be careful when
+interacting with the ecosystem if you want your results to be
+completely rigorous. Below are three things which you have to be extra
+careful with.
+
+### Implicit promotion
+Julia automatically promotes types in many cases and in particular you
+have to watch out for temporary non-rigorous values. For example
+`2(π*(Arb(ℯ)))` is okay, but not `2π*Arb(ℯ)`
+
+``` julia
+julia> 2(π*(Arb(ℯ)))
+[17.079468445347134130927101739093148990069777071530229923759202260358457222314 +/- 9.19e-76]
+
+julia> 2π*Arb(ℯ)
+[17.079468445347133465140073658536286170170195258393831755094914544308087031794 +/- 7.93e-76]
+
+julia> Arblib.overlaps(2(π*(Arb(ℯ))), 2π*Arb(ℯ))
+false
+```
+
+### Non-rigorous approximations
+In many cases this is obvious, for example Julias built in methods for
+solving linear systems will not produce rigorous results.
+
+TODO: Come up with more examples
+
+### Implementation details
+In some cases the implementation in Julia implicitly makes certain
+assumptions to improve performance and this can lead to issues. For
+example the `maximum` method in Julia checks for `NaN` results (on
+which is short fuses) using `x == x`, which works for most numerical
+types but not for `Arb` (`x == x` is only true if the radius is zero).
+See <https://github.com/JuliaLang/julia/issues/36287> for some more
+details. Arblib implements its own `maximum` method which gives
+rigorous results, but it only covers the case
+`maximum(AbstractFloat{Arb})`.
+
+``` julia
+julia> f = i -> Arb((i, i + 1));
+
+julia> A = f.(0:1000);
+
+julia> maximum(A)
+[1.00e+3 +/- 1.01]
+
+julia> maximum(A, dims = 1)
+1-element Array{Arb,1}:
+ [+/- 1.01]
+
+julia> maximum(f, 0:1000)
+[+/- 1.01]
+```
+
+These types of problems are the hardest to find since they are not
+clear from the documentation but you have to read the implementation,
+`@which` and `@less` are your friends in these cases.
+
 ## Example
 
-Here is the naive sine compuation example form the [Arb documentation](http://arblib.org/using.html#a-worked-example-the-sine-function) in Julia:
+Here is the naive sine compuation example form the [Arb
+documentation](http://arblib.org/using.html#a-worked-example-the-sine-function)
+in Julia:
 
 ```julia
 using Arblib

src/arbcalls/acb.jl

@@ -194,7 +194,7 @@ arbcall"void acb_sinh_cosh(acb_t s, acb_t c, const acb_t z, slong prec)"
 arbcall"void acb_tanh(acb_t s, const acb_t z, slong prec)"
 arbcall"void acb_coth(acb_t s, const acb_t z, slong prec)"
 arbcall"void acb_sech(acb_t res, const acb_t z, slong prec)"
-arbcall"void acb_csch(acb_t res, const arb_t z, slong prec)"
+arbcall"void acb_csch(acb_t res, const acb_t z, slong prec)"
 
 ### Inverse hyperbolic functions
 arbcall"void acb_asinh(acb_t res, const acb_t z, slong prec)"

src/arithmetic.jl

@@ -1,3 +1,7 @@
+const _BitInteger = (Int == Int64 ? Base.BitInteger64 : Base.BitInteger32)
+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)
@@ -18,63 +22,210 @@ 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) = min!(zero(x), x, y)
-Base.max(x::MagOrRef, y::MagOrRef) = 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
 
-### Arf, Arb and Acb
+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::T, y::T) where {T<:Union{Arf,ArfRef,Arb,ArbRef,Acb,AcbRef}}
-        z = T(prec = max(precision(x), precision(y)))
+    @eval function Base.$jf(x::ArfOrRef, y::Union{ArfOrRef,_BitInteger})
+        z = Arf(prec = _precision((x, y)))
         $af(z, x, y)
-        z
+        return z
     end
 end
-function Base.:(-)(x::T) where {T<:Union{Arf,ArfRef,Arb,ArbRef,Acb,AcbRef}}
-    z = T(prec = precision(x))
-    neg!(z, x)
-    z
+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
+function Base.:/(x::_BitUnsigned, y::ArfOrRef)
+    z = zero(y)
+    ui_div!(z, x, y)
+    return z
+end
+function Base.:/(x::_BitSigned, y::ArfOrRef)
+    z = zero(y)
+    si_div!(z, x, y)
+    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
+end
+function root(x::ArfOrRef, k::Integer)
+    y = zero(x)
+    root!(y, x, convert(UInt, k))
+    return y
+end
+
+### 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
-function Base.:(^)(x::T, k::Integer) where {T<:Union{Arb,ArbRef,Acb,AcbRef}}
-    z = T(prec = precision(x))
-    pow!(z, x, convert(UInt, k))
-    z
+
+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)
+
+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}) =
+    pow!(Acb(prec = _precision((x, y))), x, y)
+
+# We define the same special cases as Arb does, this avoids some
+# overhead
+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.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,
-    :exp,
-    :expm1,
     :log,
     :log1p,
+    :exp,
+    :expm1,
     :sin,
     :cos,
     :tan,
     :cot,
+    :sec,
+    :csc,
+    :atan,
     :asin,
     :acos,
-    :atan,
-    :acot,
-    :sinc,
     :sinh,
     :cosh,
     :tanh,
+    :coth,
+    :sech,
+    :csch,
+    :atanh,
     :asinh,
     :acosh,
-    :atanh,
-    :sec,
-    :asec,
-    :sech,
-    :asech,
 ]
-    @eval function Base.$f(x::T) where {T<:Union{Arb,ArbRef,Acb,AcbRef}}
-        z = T(prec = precision(x))
-        $(Symbol(f, :!))(z, x)
-        z
+    @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 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
+
+### 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 = prec), z)
 Base.imag(z::AcbLike; prec = _precision(z)) = get_imag!(Arb(prec = prec), z)
 Base.conj(z::AcbLike) = conj!(Acb(prec = _precision(z)), z)
 Base.abs(z::AcbLike) = abs!(Arb(prec = _precision(z)), z)
+
+### min and max
+for T in [MagOrRef, ArfOrRef, ArbOrRef]
+    for op in [:min, :max]
+        # min/max
+        @eval Base.$op(x::$T, y::$T) = $(Symbol(op, :!))(zero(x), x, y)
+        # minimum/maximum
+        # The default implementation of minimum/maximum in Julia is
+        # not correct for Arb, we implemented a correct version when
+        # no specific dimension is given. The default give correct
+        # results for Mag and Arf but wrong results for Arb in some
+        # cases, be careful!
+        @eval function Base.$(Symbol(op, :imum))(A::AbstractArray{<:$T})
+            isempty(A) &&
+                throw(ArgumentError("reducing over an empty collection is not allowed"))
+            res = copy(first(A))
+            for x in Iterators.drop(A, 1)
+                $(Symbol(op, :!))(res, res, x)
+            end
+            return res
+        end
+    end
+    @eval Base.minmax(x::$T, y::$T) = (min(x, y), max(x, y))
+
+    @eval function Base.extrema(A::AbstractArray{<:$T})
+        isempty(A) &&
+            throw(ArgumentError("reducing over an empty collection is not allowed"))
+        l = copy(first(A))
+        u = copy(first(A))
+        for x in Iterators.drop(A, 1)
+            min!(l, l, x)
+            max!(u, u, x)
+        end
+        return (l, u)
+    end
+end

src/constructors.jl

@@ -7,6 +7,7 @@ Mag(x, y) = set!(Mag(), x, y)
 Arf(x; prec::Integer = _precision(x)) = set!(Arf(prec = prec), x)
 # disambiguation
 Arf(x::Arf; prec::Integer = precision(x)) = set!(Arf(prec = prec), x)
+Arf(x::Rational; prec::Integer = _precision(x)) = set!(Arf(prec = prec), x)
 
 #Arb
 Arb(x; prec::Integer = _precision(x)) = set!(Arb(prec = prec), x)

src/setters.jl

@@ -19,6 +19,12 @@ function set!(res::ArfLike, x::Int128)
     return x < 0 ? neg!(res, res) : res
 end
 
+function set!(res::ArfLike, x::Rational; prec::Integer = precision(res))
+    set!(res, numerator(x))
+    div!(res, res, denominator(x), prec = prec)
+    return res
+end
+
 # Arb
 function set!(res::ArbLike, x::Union{UInt128,Int128})
     set!(radref(res), UInt64(0))

src/types.jl

@@ -253,6 +253,7 @@ for prefix in [:mag, :arf, :arb, :acb]
         parentstruct(x::$T) = cstruct(x)
         parentstruct(x::$TRef) = x
         parentstruct(x::$TStruct) = x
+        Base.copy(x::Union{$T,$TRef}) = $T(x)
     end
 end