ArbVector (first try)

[kalmarek]

Just to let you know, I tried the following:

const libarb = Nemo.libarb
import Nemo.acb_struct

mutable struct AcbVector <: AbstractVector{acb_struct}
    ptr::Ptr{acb_struct}
    length::Int
    precision::Int

    function AcbVector(n::Integer, precision::Integer)
        v = new(
            ccall((:_acb_vec_init, libarb), Ptr{acb_struct}, (Clong,), n),
            n,
            precision,
        )
        finalizer(clear!, v)
        return v
    end
end

Base.cconvert(::Type{Ptr{acb_struct}}, acb_v::AcbVector) = acb_v.ptr
Base.size(acb_v::AcbVector) = (acb_v.length,)
Base.precision(acb_v::AcbVector) = acb_v.precision

function clear!(acb_v::AcbVector)
    ccall(
        (:_acb_vec_clear, libarb),
        Cvoid,
        (Ptr{acb_struct}, Clong),
        acb_v,
        length(acb_v),
    )
end

Base.@propagate_inbounds function Base.getindex(acb_v::AcbVector, i::Integer)
    @boundscheck checkbounds(acb_v, i)
    return unsafe_load(acb_v.ptr, i)
end

_get_ptr(acb_v::AcbVector, i::Int = 1) =
    acb_v.ptr + (i - 1) * sizeof(acb_struct)

function AcbVector(v::AbstractVector{acb}, p = prec(parent(first(v))))
    acb_v = AcbVector(length(v), p)
    for (i, val) in zip(eachindex(acb_v), v)
        ccall(
            (:acb_set, libarb),
            Cvoid,
            (Ptr{acb_struct}, Ref{acb}),
            _get_ptr(acb_v, i),
            val,
        )
    end
    return acb_v
end

function approx_eig_qr!(v::AcbVector, R::acb_mat, A::acb_mat)
    ccall(
        (:acb_mat_approx_eig_qr, libarb),
        Cint,
        (
            Ptr{acb_struct},
            Ptr{Cvoid},
            Ref{acb_mat},
            Ref{acb_mat},
            Ptr{Cvoid},
            Int,
            Int,
        ),
        v,
        C_NULL,
        R,
        A,
        C_NULL,
        0,
        prec(parent(A)),
    )
    return v
end

and it works pretty well (this is still Nemo based). Eigenvalues in Nemo are incorrect so if you need them quickly here there are:

function approx_eig_qr!(v::AcbVector, R::acb_mat, A::acb_mat)
    ccall(
        (:acb_mat_approx_eig_qr, libarb),
        Cint,
        (
            Ptr{acb_struct},
            Ptr{Cvoid},
            Ref{acb_mat},
            Ref{acb_mat},
            Ptr{Cvoid},
            Int,
            Int,
        ),
        v,
        C_NULL,
        R,
        A,
        C_NULL,
        0,
        prec(parent(A)),
    )
    return v
end

function (C::AcbField)(z::acb_struct)
    res = zero(C)
    ccall((:acb_set, libarb), Cvoid, (Ref{acb}, Ref{acb_struct}), res, z)
    return res
end

function LinearAlgebra.eigvals(A::acb_mat)
    n = nrows(A)
    CC = base_ring(A)
    p = prec(CC)
    λ_approx = AcbVector(n, p)
    R_approx = similar(A)
    v = approx_eig_qr!(λ_approx, R_approx, A)

    λ = AcbVector(n, p)
    b = ccall(
        (:acb_mat_eig_multiple, libarb),
        Cint,
        (Ptr{acb_struct}, Ref{acb_mat}, Ptr{acb_struct}, Ref{acb_mat}, Int),
        λ,
        A,
        λ_approx,
        R_approx,
        p,
    )

    return CC.(λ)
end

function (C::AcbField)(z::acb_struct)
    res = zero(C)
    ccall((:acb_set, libarb), Cvoid, (Ref{acb}, Ref{acb_struct}), res, z)
    return res
end

function LinearAlgebra.eigvals(A::acb_mat)
    n = nrows(A)
    CC = base_ring(A)
    p = prec(CC)
    λ_approx = AcbVector(n, p)
    R_approx = similar(A)
    v = approx_eig_qr!(λ_approx, R_approx, A)

    λ = AcbVector(n, p)
    b = ccall(
        (:acb_mat_eig_multiple, libarb),
        Cint,
        (Ptr{acb_struct}, Ref{acb_mat}, Ptr{acb_struct}, Ref{acb_mat}, Int),
        λ,
        A,
        λ_approx,
        R_approx,
        p,
    )

    return CC.(λ)
end

---

EDIT: Updated with precision field

[Joel-Dahne]

Nice! So this would be an example where indexing does not required allocations right? And then it would be the users responsibility to not use any variables gotten from indexing after the vector is freed.

I have some ideas for how we could do the unsafe indexing in Arblib, I'll write about it in the other issue.

[Joel-Dahne]

Never mind, I'll just write it here. If we define an ArbVector as a Ptr{arb_struct} (plus some metadata) then one way to allow inplace operations would be to rewrite our arbcall to not only work on Arb but also on arb_struct. Then you could work directly on the Ptr{arb_struct} using the same low-level methods as for Arb, one issue is that we would not know which precision to use but we could maybe just default to the global precision and if you want something else you have to be explicit about it. I think this would be fine to do for the low-level methods, if you use them you should expect some friction. So you could do something like

v = ArbVector([1, 2, 3])
Arblib.add!(Arblib.unsafe_getindex(v, 1), Arblib.unsafe_getindex(v, 1), 5, prec = precision(v))

It should also work for referencing radius and midpoint

x = Arb(pi)
r = Mag(1e-10)
Arblib.cmp(x.radius, r)

The higher level interface, v[1] and radius(x), could then be written in a safe way with allocations.

[kalmarek]

unsafe_load does copy the data, so it's not allocation free.

I just tried

v = [1,2,3]
CC = AcbField(256)
jlv = CC.(v)
arbv = AcbVector(jlv)
unsafe_wrap(Array, arbv.ptr, (3,)) # ← truly non-allocating

but it crashes instantly :D

[Joel-Dahne]

Okay, so it does allocation for the struct but does not call arb_init and thus doesn't need to allocate space for the pointers in the struct at least? I hope to have time to experiment a bit more with this in the weekend.

[kalmarek]

@Joel-Dahne that is correct: unsafe_load allocates a julia object and copies the corresponding fields to it. However the arblibs data (the fields are pointing to) are not copied, therefore the result is an unsafe object.
To be precise, this is the correct syntax to add 5 to the first coordinate:

v = ArbVector([1, 2, 3])
GC.@preserve v
    Arblib.add!(Arblib.unsafe_getindex(v, 1), Arblib.unsafe_getindex(v, 1), 5, prec = precision(v))
    b = Arblib.unsafe_getindex(v, 1)
    # do something with b, but never leak it out of @preserve block:
    (b > 1 ? true : false)
end

I'm not sure if unsafe_getindex is the right name here: in julia unsafe_* means that the arguments need to be checked for validity, but the result is a safe julia object. Here, it's the opposite...

[Joel-Dahne]

One problem I see if we copy the struct is that some of the data is stored directly in the struct and some of it where the fields are pointing to, which data is stored in the struct also depends on the data. So if we copy the struct and update the fields the original copy might or might not be updated correctly.

[Joel-Dahne]

The struct for ArbVector, AcbVector and the such should not really depend on how we choose to handle mutability and non-allocation right? What I might do is to start implementing the simpler parts, defining the types and such.