Merge pull request #158 from kalmarek/serialize-hash

Add hash and serialization methods

.github/workflows/format.yml

@@ -18,7 +18,7 @@ jobs:
       - uses: actions/checkout@v2
       - name: Install JuliaFormatter and format
         run: |
-          julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.22.10"))'
+          julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="1.0.9"))'
           julia  -e 'using JuliaFormatter; format(["./src", "./test"], verbose=true)'
       - name: Format check
         run: |

Project.toml

@@ -1,13 +1,14 @@
 name = "Arblib"
 uuid = "fb37089c-8514-4489-9461-98f9c8763369"
 authors = ["Marek Kaluba <[email protected]>", "Sascha Timme <Sascha Timme <[email protected]>", "Joel Dahne <[email protected]>"]
-version = "0.8.0"
+version = "0.8.1"
 
 [deps]
 Arb_jll = "d9960996-1013-53c9-9ba4-74a4155039c3"
 FLINT_jll = "e134572f-a0d5-539d-bddf-3cad8db41a82"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 

src/Arblib.jl

@@ -2,6 +2,7 @@ module Arblib
 
 using Arb_jll
 import LinearAlgebra
+import Serialization
 import SpecialFunctions
 
 # So that the parsed contains method extends the base function
@@ -52,8 +53,11 @@ function flint_set_num_threads(a::Integer)
 end
 
 include("arb_types.jl")
+include("fmpz.jl")
 include("rounding.jl")
 include("types.jl")
+include("hash.jl")
+include("serialize.jl")
 
 include("ArbCall/ArbCall.jl")
 import .ArbCall: @arbcall_str, @arbfpwrapcall_str

src/fmpz.jl

@@ -0,0 +1,25 @@
+"""
+    fmpz_struct()
+
+Low level wrapper of `fmpz_t`. Not part of the Arblib interface but
+only used internally by a few methods for conversion from `fmpz_t` to
+`BigInt`.
+"""
+mutable struct fmpz_struct
+    d::Int
+
+    function fmpz_struct()
+        z = new()
+        ccall(@libflint(fmpz_init), Nothing, (Ref{fmpz_struct},), z)
+        finalizer(fmpz_clear!, z)
+        return z
+    end
+end
+
+fmpz_clear!(x::fmpz_struct) = ccall(@libflint(fmpz_clear), Nothing, (Ref{fmpz_struct},), x)
+
+function Base.BigInt(x::fmpz_struct)
+    res = BigInt()
+    ccall(@libflint(fmpz_get_mpz), Nothing, (Ref{BigInt}, Ref{fmpz_struct}), res, x)
+    return res
+end

src/hash.jl

@@ -0,0 +1,98 @@
+#=
+The Base implementation of hash(::Real) is based on Base.decompose
+to guarantee that all numbers that compare equal are given the same
+hash. We implement Base.decompose for Mag and Arf.
+
+It should return a, not necessarily canonical, decomposition of
+rational values as `num*2^pow/den`. For Mag and Arf we always have den
+= 1 and we hence only need to find num and pow, corresponding to the
+mantissa end exponent. For Arf this is straight forward using
+arf_get_fmpz_2exp. For Mag the mantissa is stored directly in the
+struct as a UInt and the exponent as a fmpz.
+=#
+function Base.decompose(x::Union{mag_struct,Ptr{mag_struct}})::Tuple{UInt,BigInt,Int}
+    isinf(x) && return 1, 0, 0
+
+    if x isa Ptr{mag_struct}
+        x = unsafe_load(x)
+    end
+
+    pow = BigInt()
+    ccall(@libflint(fmpz_get_mpz), Nothing, (Ref{BigInt}, Ref{UInt}), pow, x.exponent)
+    # There is an implicit factor 2^30 for the exponent, coming from
+    # the number of bits of the mantissa
+    pow -= 30
+
+    return x.mantissa, pow, 1
+end
+
+function Base.decompose(x::Union{arf_struct,Ptr{arf_struct}})::Tuple{BigInt,BigInt,Int}
+    isnan(x) && return 0, 0, 0
+    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
+
+    num = fmpz_struct()
+    pow = fmpz_struct()
+    ccall(
+        @libarb(arf_get_fmpz_2exp),
+        Cvoid,
+        (Ref{fmpz_struct}, Ref{fmpz_struct}, Ref{arf_struct}),
+        num,
+        pow,
+        x,
+    )
+
+    return BigInt(num), BigInt(pow), 1
+end
+
+Base.decompose(x::Union{MagOrRef,ArfOrRef}) = Base.decompose(cstruct(x))
+
+# Hashes of structs are computed using the method for the wrapping
+# type
+Base.hash(x::mag_struct, h::UInt) = hash(Mag(x), h)
+Base.hash(x::arf_struct, h::UInt) = hash(Arf(x), h)
+Base.hash(x::arb_struct, h::UInt) = hash(Arb(x), h)
+Base.hash(x::acb_struct, h::UInt) = hash(Acb(x), h)
+
+# Hashes of Mag and Arf are computed using the Base implementation
+# which used Base.decompose defined above.
+
+function Base.hash(x::ArbLike, h::UInt)
+    # If the radius is zero we compute the hash using only the
+    # midpoint, so that we get identical hashes as for the
+    # corresponding Arf
+    if !isexact(x)
+        h = hash(Arblib.radref(x), h)
+    end
+    return hash(Arblib.midref(x), h)
+end
+
+function Base.hash(z::AcbLike, h::UInt)
+    # Same as for Complex{T}
+    hash(realref(z), h ⊻ hash(imagref(z), Base.h_imag) ⊻ Base.hash_0_imag)
+end
+
+# Compare with Base.h_imag
+if UInt === UInt64
+    const h_poly = 0xfd6de1a6c0e66975
+else
+    const h_poly = 0xa0617887
+end
+# arb_poly_struct and acb_poly_struct use default hash implementation,
+# this is okay since they don't implement an isequal method.
+function Base.hash(p::Union{ArbPoly,AcbPoly}, h::UInt)
+    h = hash(h_poly, h)
+    for i = 0:degree(p)
+        h = hash(ref(p, i), h)
+    end
+    return h
+end
+
+function Base.hash(p::Union{ArbSeries,AcbSeries}, h::UInt)
+    # Conversion of Number to series gives a degree 0 series, we want
+    # the hashes to match in this case
+    degree(p) == 0 && return hash(ref(p, 0), h)
+
+    hash(p.poly, hash(degree(p), h))
+end
+
+# Vectors and Matrices have an implementation in Base that works well

src/matrix.jl

@@ -106,7 +106,7 @@ end
 
 # General constructor
 for T in [:ArbMatrix, :ArbRefMatrix, :AcbMatrix, :AcbRefMatrix]
-    @eval function $T(A::AbstractMatrix; prec::Integer = _precision(first(A)))
+    @eval function $T(A::AbstractMatrix; prec::Integer = _precision(A))
         B = $T(size(A)...; prec = prec)
         # ensure to handle all kind of indices
         ax1, ax2 = axes(A)
@@ -116,7 +116,7 @@ for T in [:ArbMatrix, :ArbRefMatrix, :AcbMatrix, :AcbRefMatrix]
         return B
     end
 
-    @eval function $T(v::AbstractVector; prec::Integer = _precision(first(v)))
+    @eval function $T(v::AbstractVector; prec::Integer = _precision(v))
         A = $T(length(v), 1; prec = prec)
         for (i, vᵢ) in enumerate(v)
             A[i, 1] = vᵢ

src/poly.jl

@@ -129,7 +129,7 @@ for (TPoly, TSeries) in [(:ArbPoly, :ArbSeries), (:AcbPoly, :AcbSeries)]
 
     @eval function $TPoly(
         coeffs::Union{Tuple,AbstractVector};
-        prec::Integer = _precision(first(coeffs)),
+        prec::Integer = _precision(coeffs),
     )
         p = fit_length!($TPoly(; prec), length(coeffs))
         @inbounds for (i, c) in enumerate(coeffs)
@@ -142,7 +142,7 @@ for (TPoly, TSeries) in [(:ArbPoly, :ArbSeries), (:AcbPoly, :AcbSeries)]
     # with two elements. This would for example be used when
     # constructing a polynomial with a constant plus x, e.g
     # ArbPoly((x, 1))
-    @eval function $TPoly(coeffs::Tuple{Any,Any}; prec::Integer = _precision(first(coeffs)))
+    @eval function $TPoly(coeffs::Tuple{Any,Any}; prec::Integer = _precision(coeffs))
         p = fit_length!($TPoly(; prec), length(coeffs))
         @inbounds p[0] = coeffs[1]
         @inbounds p[1] = coeffs[2]
@@ -175,8 +175,8 @@ for (TSeries, TPoly) in [(:ArbSeries, :ArbPoly), (:AcbSeries, :AcbPoly)]
 
     @eval function $TSeries(
         coeffs::Union{Tuple,AbstractVector};
-        degree::Integer = length(coeffs) - 1,
-        prec::Integer = _precision(first(coeffs)),
+        degree::Integer = max(length(coeffs) - 1, 0),
+        prec::Integer = _precision(coeffs),
     )
         p = $TSeries(; degree, prec)
         @inbounds for (i, c) in enumerate(coeffs)
@@ -193,7 +193,7 @@ for (TSeries, TPoly) in [(:ArbSeries, :ArbPoly), (:AcbSeries, :AcbPoly)]
     @eval function $TSeries(
         coeffs::Tuple{Any,Any};
         degree::Integer = length(coeffs) - 1,
-        prec::Integer = _precision(first(coeffs)),
+        prec::Integer = _precision(coeffs),
     )
         p = $TSeries(; degree, prec)
         @inbounds p[0] = coeffs[1]

src/precision.jl

@@ -22,6 +22,8 @@ Base.precision(x::Union{ArbSeries,AcbSeries}) = precision(x.poly)
 
 @inline _precision(x::Union{ArbTypes,BigFloat}) = precision(x)
 @inline _precision(z::Complex) = max(_precision(real(z)), _precision(imag(z)))
+@inline _precision(v::Union{Tuple,AbstractVector}) =
+    isempty(v) ? DEFAULT_PRECISION[] : _precision(first(v))
 @inline _precision(
     a::Union{ArbTypes,BigFloat,Complex{<:Union{ArbTypes,BigFloat}}},
     b::Union{ArbTypes,BigFloat,Complex{<:Union{ArbTypes,BigFloat}}},

src/serialize.jl

@@ -0,0 +1,125 @@
+# Compare with BigInt in julia/stdlib/v1.7/Serialization/src/Serialization.jl
+function Serialization.serialize(
+    s::Serialization.AbstractSerializer,
+    x::Union{mag_struct,arf_struct,arb_struct},
+)
+    Serialization.serialize_type(s, typeof(x))
+    Serialization.serialize(s, dump_string(x))
+end
+
+function Serialization.serialize(s::Serialization.AbstractSerializer, x::acb_struct)
+    Serialization.serialize_type(s, typeof(x))
+    str = dump_string(Arblib.realref(x)) * " " * dump_string(Arblib.imagref(x))
+    Serialization.serialize(s, str)
+end
+
+function Serialization.serialize(
+    s::Serialization.AbstractSerializer,
+    v::Union{arb_vec_struct,acb_vec_struct},
+)
+    Serialization.serialize_type(s, typeof(v))
+    Serialization.serialize(s, size(v)[1])
+    for i = 1:size(v)[1]
+        Serialization.serialize(s, unsafe_load(v[i]))
+    end
+end
+
+function Serialization.serialize(
+    s::Serialization.AbstractSerializer,
+    v::Union{arb_mat_struct,acb_mat_struct},
+)
+    Serialization.serialize_type(s, typeof(v))
+    Serialization.serialize(s, size(v))
+    for i = 1:size(v)[1]
+        for j = 1:size(v)[2]
+            Serialization.serialize(s, unsafe_load(v[i, j]))
+        end
+    end
+end
+
+function Serialization.serialize(
+    s::Serialization.AbstractSerializer,
+    p::Union{arb_poly_struct,acb_poly_struct},
+)
+    Serialization.serialize_type(s, typeof(p))
+    Serialization.serialize(s, length(p))
+    T = p isa arb_poly_struct ? arb_struct : acb_struct
+    for i = 0:length(p)-1
+        Serialization.serialize(s, unsafe_load(p.coeffs + i * sizeof(T)))
+    end
+end
+
+
+Serialization.deserialize(
+    s::Serialization.AbstractSerializer,
+    T::Type{<:Union{mag_struct,arf_struct,arb_struct}},
+) = Arblib.load_string!(T(), Serialization.deserialize(s))
+
+function Serialization.deserialize(s::Serialization.AbstractSerializer, T::Type{acb_struct})
+    str = Serialization.deserialize(s)
+    # Three spaces in the real part, so we are looking for
+    spaces = findall(" ", str)
+    @assert length(spaces) == 7
+
+    real_str = str[1:spaces[4].start-1]
+    imag_str = str[spaces[4].stop+1:end]
+
+    res = acb_struct()
+    Arblib.load_string!(Arblib.realref(res), real_str)
+    Arblib.load_string!(Arblib.imagref(res), imag_str)
+    return res
+end
+
+function Serialization.deserialize(
+    s::Serialization.AbstractSerializer,
+    T::Type{<:Union{arb_vec_struct,acb_vec_struct}},
+)
+    n = Serialization.deserialize(s)
+    res = T(n)
+    for i = 1:n
+        res[i] = Serialization.deserialize(s)
+    end
+    return res
+end
+
+function Serialization.deserialize(
+    s::Serialization.AbstractSerializer,
+    T::Type{<:Union{arb_mat_struct,acb_mat_struct}},
+)
+    r, c = Serialization.deserialize(s)
+    res = T(r, c)
+    for i = 1:r
+        for j = 1:c
+            res[i, j] = Serialization.deserialize(s)
+        end
+    end
+    return res
+end
+
+function Serialization.deserialize(
+    s::Serialization.AbstractSerializer,
+    T::Type{<:Union{arb_poly_struct,acb_poly_struct}},
+)
+    n = Serialization.deserialize(s)
+    res = T()
+    for i = 0:n-1
+        set_coeff!(res, i, Serialization.deserialize(s))
+    end
+    return res
+end
+
+# For series we want to make sure that we have allocated the right
+# number of coefficients
+function Serialization.deserialize(
+    s::Serialization.AbstractSerializer,
+    T::Type{<:Union{ArbSeries,AcbSeries}},
+)
+    # This uses the base implementation of serialize so if that
+    # changes this might need to be updated
+    poly = Serialization.deserialize(s)
+    degree = Serialization.deserialize(s)
+    # The series constructor assures us that the number of allocated
+    # coefficients is determined from the degree
+    res = T(poly; degree) # TODO: Here an inplace constructor could be used
+    return res
+end

src/show.jl

@@ -106,21 +106,17 @@ for ArbT in (Mag, MagRef, Arf, ArfRef, Arb, ArbRef, Acb, AcbRef)
     end
 end
 
-for ArbT in (Mag, MagRef, Arf, ArfRef, Arb, ArbRef)
-    @eval begin
-        function load_string!(x::$ArbT, str::AbstractString)
-            res = load!(x, str)
-            iszero(res) || throw(ArgumentError("could not load $str as " * $(string(ArbT))))
-            return x
-        end
+function load_string!(x::Union{MagLike,ArfLike,ArbLike}, str::AbstractString)
+    res = load!(x, str)
+    iszero(res) || throw(ArgumentError("could not load $str as $(string(typeof(x)))"))
+    return x
+end
 
-        function dump_string(x::$ArbT)
-            char_ptr = dump(x)
-            str = unsafe_string(char_ptr)
-            ccall(@libflint(flint_free), Cvoid, (Cstring,), char_ptr)
-            return str
-        end
-    end
+function dump_string(x::Union{MagLike,ArfLike,ArbLike})
+    char_ptr = dump(x)
+    str = unsafe_string(char_ptr)
+    ccall(@libflint(flint_free), Cvoid, (Cstring,), char_ptr)
+    return str
 end
 
 for T in [