Merge pull request #173 from kalmarek/extend-interval-methods

Add midpoint for Acb and union and intersection for polynomials

src/interval.jl

@@ -22,6 +22,18 @@ midpoint(::Type{Arf}, x::ArbOrRef) = Arf(midref(x))
 midpoint(::Type{Arb}, x::ArbOrRef) = Arb(midref(x))
 midpoint(x::ArbOrRef) = midpoint(Arf, x)
 
+"""
+    midpoint([T, ] z::AcbOrRef)
+
+Returns the midpoint of `z` as a `Complex{Arf}`. If `T` is given and
+equal to `Arf` or `Arb`, convert to `Complex{T}`. If `T` is `Acb` then
+convert to that.
+"""
+midpoint(::Type{Acb}, z::AcbOrRef) = Acb(midref(realref(z)), midref(imagref(z)))
+midpoint(T::Type{<:Union{Arf,Arb}}, z::AcbOrRef) =
+    Complex(midpoint(T, realref(z)), midpoint(T, imagref(z)))
+midpoint(z::AcbOrRef) = midpoint(Arf, z)
+
 """
     lbound([T, ] x::ArbOrRef)
 
@@ -131,15 +143,26 @@ getball(::Type{Arb}, x::ArbOrRef) = (Arb(midref(x)), Arb(radref(x), prec = preci
 """
     union(x::ArbOrRef, y::ArbOrRef)
     union(x::AcbOrRef, y::AcbOrRef)
+    union(x::T, y::T) where {T<:Union{ArbPoly,AcbPoly,ArbSeries,AcbSeries}}
     union(x, y, z...)
 
-`union(x, y)` returns a ball containing the union of `x` and `y`.
+Returns a ball containing the union of `x` and `y`. For polynomials
+and series the union is taken coefficient-wise.
 
 `union(x, y, z...)` returns a ball containing the union of all given
 balls.
 """
 union(x::ArbOrRef, y::ArbOrRef) = union!(Arb(prec = _precision(x, y)), x, y)
 union(x::AcbOrRef, y::AcbOrRef) = union!(Acb(prec = _precision(x, y)), x, y)
+union(x::T, y::T) where {T<:Union{ArbPoly,AcbPoly}} =
+    _union!(T(prec = _precision(x, y)), x, y)
+function union(x::T, y::T) where {T<:Union{ArbSeries,AcbSeries}}
+    degree(x) == degree(y) || throw(ArgumentError("union of series requires same degree"))
+    res = T(degree = degree(x), prec = _precision(x, y))
+    _union!(res.poly, x.poly, y.poly)
+    return res
+end
+
 function union(x::ArbOrRef, y::ArbOrRef, z::ArbOrRef, xs...)
     res = union(y, z, xs...)
     return union!(res, res, x)
@@ -148,18 +171,56 @@ function union(x::AcbOrRef, y::AcbOrRef, z::AcbOrRef, xs...)
     res = union(y, z, xs...)
     return union!(res, res, x)
 end
+function union(x::T, y::T, z::T, xs...) where {T<:Union{ArbPoly,AcbPoly}}
+    res = union(y, z, xs...)
+    return _union!(res, res, x)
+end
+function union(x::T, y::T, z::T, xs...) where {T<:Union{ArbSeries,AcbSeries}}
+    res = union(y, z, xs...)
+    degree(res) == degree(x) || throw(ArgumentError("union of series requires same degree"))
+    _union!(res.poly, res.poly, x.poly)
+    return res
+end
+
+# Used internally by union
+function _union!(res::T, x::T, y::T) where {T<:Union{ArbPoly,AcbPoly}}
+    res_length = max(length(x), length(y))
+    common_degree = min(degree(x), degree(y))
+
+    fit_length!(res, res_length)
+
+    for i = 0:common_degree
+        union!(ref(res, i), ref(x, i), ref(y, i))
+    end
+
+    if common_degree + 1 < res_length
+        z = zero(eltype(T))
+        # At most one of the below loops will run
+        for i = common_degree+1:degree(x)
+            union!(ref(res, i), ref(x, i), z)
+        end
+        for i = common_degree+1:degree(y)
+            union!(ref(res, i), ref(y, i), z)
+        end
+    end
+
+    set_length!(res, res_length)
+
+    return res
+end
 
 """
     intersection(x::ArbOrRef, y::ArbOrRef)
-    intersection(x::AcbOrRef, y::AcbOrRef)
+    intersection(x::T, y::T) where {T<:Union{ArbPoly,ArbSeries}}
     intersection(x, y, z...)
 
 `intersection(x, y)` returns a ball containing the intersection of `x`
 and `y`. If `x` and `y` do not overlap (as given by `overlaps(a, b)`)
-throws an `ArgumentError`.
+throws an `ArgumentError`. For polynomials and series the intersection
+is taken coefficient-wise.
 
-`intersection(x, y, z...)` returns a ball containing the intersection of
-all given balls. If all the balls do not overlap throws an
+`intersection(x, y, z...)` returns a ball containing the intersection
+of all given balls. If all the balls do not overlap throws an
 `ArgumentError`.
 """
 function intersection(x::ArbOrRef, y::ArbOrRef)
@@ -169,13 +230,72 @@ function intersection(x::ArbOrRef, y::ArbOrRef)
         throw(ArgumentError("intersection of non-intersecting balls not allowed"))
     return res
 end
+intersection(x::ArbPoly, y::ArbPoly) =
+    _intersection!(ArbPoly((prec = _precision(x, y))), x, y)
+function intersection(x::ArbSeries, y::ArbSeries)
+    degree(x) == degree(y) ||
+        throw(ArgumentError("intersection of series requires same degree"))
+    res = ArbSeries(degree = degree(x), prec = _precision(x, y))
+    _intersection!(res.poly, x.poly, y.poly)
+    return res
+end
+
 function intersection(x::ArbOrRef, y::ArbOrRef, z::ArbOrRef, xs...)
     res = intersection(y, z, xs...)
     sucess = intersection!(res, res, x)
     iszero(sucess) &&
         throw(ArgumentError("intersection of non-intersecting balls not allowed"))
     return res
 end
+function intersection(x::ArbPoly, y::ArbPoly, z::ArbPoly, xs...)
+    res = intersection(y, z, xs...)
+    return _intersection!(res, res, x)
+end
+function intersection(x::ArbSeries, y::ArbSeries, z::ArbSeries, xs...)
+    res = intersection(y, z, xs...)
+    degree(res) == degree(x) ||
+        throw(ArgumentError("intersection of series requires same degree"))
+    _intersection!(res.poly, res.poly, x.poly)
+    return res
+end
+
+# Used internally by intersection
+function _intersection!(res::ArbPoly, x::ArbPoly, y::ArbPoly)
+    res_length = max(length(x), length(y))
+    common_degree = min(degree(x), degree(y))
+
+    fit_length!(res, res_length)
+
+    for i = 0:common_degree
+        sucess = intersection!(ref(res, i), ref(x, i), ref(y, i))
+        if iszero(sucess)
+            set_length!(res, res_length)
+            normalise!(res)
+            throw(ArgumentError("intersection of non-intersecting balls not allowed"))
+        end
+    end
+
+    if common_degree + 1 < res_length
+        # At most one of the below loops will run
+        for i = common_degree+1:degree(x)
+            xi = ref(x, i)
+            contains_zero(xi) ||
+                throw(ArgumentError("intersection of non-intersecting balls not allowed"))
+            isnan(midref(xi)) && indeterminate!(ref(res, i))
+        end
+        for i = common_degree+1:degree(y)
+            yi = ref(y, i)
+            contains_zero(yi) ||
+                throw(ArgumentError("intersection of non-intersecting balls not allowed"))
+            isnan(midref(yi)) && indeterminate!(ref(res, i))
+        end
+    end
+
+    set_length!(res, res_length)
+    normalise!(res)
+
+    return res
+end
 
 """
     add_error(x, err)

src/poly.jl

@@ -80,35 +80,92 @@ end
     ref(p::Union{ArbPoly,ArbSeries,AcbPoly,AcbSeries}, i)
 
 Similar to `p[i]` but instead of an `Arb` or `Acb` returns an `ArbRef`
-or `AcbRef` which still shares the memory with the `i`-th entry of
-`p[i]`.
-
-For `ArbPoly` and `AcbPoly` it only allows accessing coefficients up
-to the degree of the polynomial since higher coefficients are not
-guaranteed to be allocated.
-
-For `ArbSeries` and `AcbSeries` it allows accessing coefficients up to
-the degree of the series. Note that this degree might not be the same
-as the degree of the underlying polynomial (in case higher order
-coefficients are zero) but the way they are constructed ensures that
-the coefficients will always be initialised.
-
-!!! Note: If you use this to change the coefficient in a way so that
-    the degree of the polynomial might change you need to normalise
-    the polynomial afterwards to make sure that Arb recognises the
-    possibly new degree of the polynomial. If the new degree is the
-    same or lower this can be done using `Arblib.normalise!`. If the
-    new degree is higher you need to manually set the correct value of
-    `Arblib.cstruct(p).length` to be one higher than the new degree.
+or `AcbRef` which shares the memory with the `i`-th coefficient of
+`p`.
+
+!!! Note: For reading coefficients this is always safe, but if the
+    coefficient is mutated then care has to be taken. See the comment
+    further down for how to handle this.
+
+It only allows accessing coefficients that are allocated. For
+`ArbPoly` and `AcbPoly` this is typically all coefficients up to the
+degree of the polynomial, but can be higher if for example
+`Arblib.fit_length!` is used. For `ArbSeries` and `AcbSeries` all
+coefficients up to the degree of the series are guaranteed to be
+allocated, even if the underlying polynomial has a lower degree.
+
+If the coefficient is mutated in a way that the degree of the
+polynomial changes then one needs to also update the internally stored
+length of the polynomial.
+
+- If the new degree is the same or lower this can be achieved with
+  ```
+  Arblib.normalise!(p)
+  ```
+- If the new degree is higher you need to manually set the length.
+  This can be achieved with
+  ```
+  Arblib.set_length!(p, len)
+  Arblib.normalise!(p)
+  ```
+  where `len` is one higher than (an upper bound of) the new degree.
+
+See the extended help for more details.
+
+# Extended help
+Here is an example were the leading coefficient is mutated so that the
+degree is lowered.
+```jldoctest
+julia> p = ArbPoly([1, 2], prec = 64) # Polynomial of degree 1
+1.000000000000000000 + 2.000000000000000000⋅x
+
+julia> Arblib.zero!(Arblib.ref(p, 1)) # Set leading coefficient to 0
+0
+
+julia> Arblib.degree(p) # The degree is still reported as 1
+1
+
+julia> length(p) # And the length as 2
+2
+
+julia> p # Printing gives weird results
+1.000000000000000000 +
+
+julia> Arblib.normalise!(p) # Normalising the polynomial updates the degree
+1.000000000000000000
+
+julia> Arblib.degree(p) # This is now correct
+0
+
+julia> p # And so is printing
+1.000000000000000000
+```
+Here is an example when a coefficient above the degree is mutated.
+```jldoctest
+julia> q = ArbSeries([1, 2, 0], prec = 64) # Series of degree 3 with leading coefficient zero
+1.000000000000000000 + 2.000000000000000000⋅x + 𝒪(x^3)
+
+julia> Arblib.one!(Arblib.ref(q, 2)) # Set the leading coefficient to 1
+1.000000000000000000
+
+julia> q # The leading coefficient is not picked up
+1.000000000000000000 + 2.000000000000000000⋅x + 𝒪(x^3)
+
+julia> Arblib.degree(q.poly) # The degree of the underlying polynomial is still 1
+1
+
+julia> Arblib.set_length!(q, 3) # After explicitly setting the length to 3 it is ok
+1.000000000000000000 + 2.000000000000000000⋅x + 1.000000000000000000⋅x^2 + 𝒪(x^3)
+```
 """
 Base.@propagate_inbounds function ref(p::Union{ArbPoly,ArbSeries}, i::Integer)
-    @boundscheck 0 <= i <= degree(p) || throw(BoundsError(p, i))
+    @boundscheck 0 <= i < cstruct(p).alloc || throw(BoundsError(p, i))
     ptr = cstruct(p).coeffs + i * sizeof(arb_struct)
     return ArbRef(ptr, precision(p), cstruct(p))
 end
 
 Base.@propagate_inbounds function ref(p::Union{AcbPoly,AcbSeries}, i::Integer)
-    @boundscheck 0 <= i <= degree(p) || throw(BoundsError(p, i))
+    @boundscheck 0 <= i < cstruct(p).alloc || throw(BoundsError(p, i))
     ptr = cstruct(p).coeffs + i * sizeof(acb_struct)
     return AcbRef(ptr, precision(p), cstruct(p))
 end