cfunction acb_calc_func_t

[kalmarek]

(If You think that this is not the place to ask for help in wrapping Arb I can move the discussion e.g. to discourse.julialang.org)

I would like to use brand new and cool acb_calc_integrate to rigorously integrate functions defined in julia (passed via cfunction). So far I've been unable to do so, (this is an excuse for me to learn some lower level stuff, so please be patient:)

looking at acb_calc.h I see

int
acb_calc_integrate(acb_t res, acb_calc_func_t f, void * param,
    const acb_t a, const acb_t b,
    slong goal, const mag_t tol,
    const acb_calc_integrate_opt_t options,
    slong prec);

1. Reading arguments list in reverse we need a struct acb_calc_integrate_opt:

   struct acb_calc_integrate_opts
   deg_limit::Int
   eval_limit::Int
   depth_limit::Int
   use_heap::Int32
   verbose::Int32
   
   function acb_calc_integrate_opts(deg_limit::Int, eval_limit::Int, depth_limit::Int,
       use_heap::Int32, verbose::Int32)
       return new(deg_limit, eval_limit, depth_limit, use_heap, verbose)
   end
   
   function acb_calc_integrate_opts()
       opts = new()
       ccall((:acb_calc_integrate_opt_init, "/usr/lib/libarb.so"),
           Void, (Ref{acb_calc_integrate_opts}, ), opts)
       return opts
   end
   end

2. Leaving out acb_calc_func_t, I translated call to :acb_calc_integrate as

   ccall((:acb_calc_integrate, "/usr/lib/libarb.so"), UInt, 
       (Ref{acb}, #res 
       Ptr{Void}, #f
       Ptr{Void}, #params
       Ref{acb},  #a
       Ref{acb},  #b
       Int,       #rel_goal
       Ref{Nemo.mag_struct},  #abs_tol
       Ref{acb_calc_integrate_opts}, #opts
       Int),
       res, f, C_NULL, FF(-1), FF(1), prec(FF), mag, opts, prec(FF))

   where, e.g.

   FF = AcbField(400)
   res = Ref{acb}(FF())
   mag = Nemo.mag_struct(prec(FF), UInt(2))
   opts = acb_calc_integrate_opts()

3. According to (manual)[https://docs.julialang.org/en/stable/manual/calling-c-and-fortran-code/#Passing-Pointers-for-Modifying-Inputs-1] I should pass Ref{acb}(...), for the ccall above will modify it. Am I correct so far?

4. Now comes the acb_calc_func_t f: it needs to be a pointer to a

   typedef int (*acb_calc_func_t)(acb_ptr out,
   const acb_t inp, void * param, slong order, slong prec);

   I tried

   
   square(x::acb) = x*x
   function square_acb_func(res::Ref{acb}, x::acb, param::Ptr{Void}, order::Int, prec::Int)
   for fn in fieldnames(x)
       setfield!(res[], fn, convert(fieldtype(res[], fn), getfield(square(x), fn)))
   end
   return convert(Cint, 0)::Cint
   end

csquare_acb_func = cfunction(square_acb_func, Cint, (Ref{acb}, Ref{acb}, Ptr{Void}, Int, Int))

(btw, can You modify an `acb` in-place?). 

Such `csquare_acb_func`  works fine when `ccalled`:
```julia
FF = AcbField(400)
x=FF(3)
res = Ref(FF(2))

@show res[]
@time ccall(cf_acb, Cint, (Ref{acb}, Ref{acb}, Ptr{Void}, Int, Int), 
                     res,      x,        C_NULL,    1,   prec(parent(x)))
@show res[];

res[] = 2.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 + i0 0.000063 seconds (21 allocations: 608 bytes) res[] = 9.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 + i0

but segfaults when placed into ccall to :acb_calc_integrate. I guess this is because I don't understand how acb_ptr and acb_t interact with julia types (I translated both to Ref{acb}). Could You give me a tip?

[wbhart]

On 5 February 2018 at 14:11, kalmarek notifications@github.com wrote:

(If You think that this is not the place to ask for help in wrapping Arb I can move the discussion e.g. to discourse.julialang.org)

I would like to use brand new and cool acb_calc_integrate to rigorously integrate functions defined in julia (passed via cfunction). So far I've been unable to do so, (this is an excuse for me to learn some lower level stuff, so please be patient:)

looking at acb_calc.h I see

intacb_calc_integrate(acb_t res, acb_calc_func_t f, void * param, const acb_t a, const acb_t b, slong goal, const mag_t tol, const acb_calc_integrate_opt_t options, slong prec);

1. Reading arguments list in reverse we need a struct acb_calc_integrate_opt:

struct acb_calc_integrate_opts

You will almost certainly need "mutable struct", not "struct".

deg_limit::Int
eval_limit::Int
depth_limit::Int
use_heap::Int32
verbose::Int32

function acb_calc_integrate_opts(deg_limit::Int, eval_limit::Int, depth_limit::Int,
    use_heap::Int32, verbose::Int32)
    return new(deg_limit, eval_limit, depth_limit, use_heap, verbose)
end

function acb_calc_integrate_opts()
    opts = new()
    ccall((:acb_calc_integrate_opt_init, "/usr/lib/libarb.so"),

You can use libarb in Nemo, instead of

"/usr/lib/libarb.so"

        Void, (Ref{acb_calc_integrate_opts}, ), opts)

    return opts
endend

1. Leaving out acb_calc_func_t, I translated call to :acb_calc_integrate as

ccall((:acb_calc_integrate, "/usr/lib/libarb.so"), UInt, (Ref{acb}, #res Ptr{Void}, #f Ptr{Void}, #params Ref{acb}, #a Ref{acb}, #b Int, #rel_goal Ref{Nemo.mag_struct}, #abs_tol Ref{acb_calc_integrate_opts}, #opts Int), res, f, C_NULL, FF(-1), FF(1), prec(FF), mag, opts, prec(FF))

where, e.g.

FF = AcbField(400) res = Ref{acb}(FF()) mag = Nemo.mag_struct(prec(FF), UInt(2)) opts = acb_calc_integrate_opts()

1.

According to (manual)[https://docs.julialang.org/en/stable/ manual/calling-c-and-fortran-code/#Passing-Pointers-for- Modifying-Inputs-1 https://docs.julialang.org/en/stable/manual/calling-c-and-fortran-code/#Passing-Pointers-for-Modifying-Inputs-1] I should pass Ref{acb}(...), for the ccall above will modify it. Am I correct so far?

Yes, I think so. We now use Ref instead of the old syntax.

1. 2.

Now comes the acb_calc_func_t f: it needs to be a pointer to a

typedef int (acb_calc_func_t)(acb_ptr out, const acb_t inp, void param, slong order, slong prec);

I tried

square(x::acb) = x*xfunction square_acb_func(res::Ref{acb}, x::acb, param::Ptr{Void}, order::Int, prec::Int) for fn in fieldnames(x) setfield!(res[], fn, convert(fieldtype(res[], fn), getfield(square(x), fn))) end return convert(Cint, 0)::Cintend

csquare_acb_func = cfunction(square_acb_func, Cint, (Ref{acb}, Ref{acb}, Ptr{Void}, Int, Int))

(btw, can You modify an acb in-place?).

Yes, you can. Of course you need to make sure you own the acb, i.e. that there are no other references to it in the system, otherwise all other references will be modified too.

Usually this means creating a fresh acb that can be modified. The exception is if the arb is known to be the result of a nontrivial coercion or that it is the result of an arithmetic operation. These both create fresh arbs with no other references.

Such csquare_acb_func works fine when ccalled:

FF = AcbField(400) x=FF(3) res = Ref(FF(2)) @show res[]@time ccall(cf_acb, Cint, (Ref{acb}, Ref{acb}, Ptr{Void}, Int, Int), res, x, C_NULL, 1, prec(parent(x)))@show res[];

res[] = 2.000000000000000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000000000000000 + i

*0 0.000063 seconds (21 allocations: 608 bytes) res[] = 9.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

• i*0

but segfaults when placed into ccall to :acb_calc_integrate. I guess this is because I don't understand how acb_ptr and acb_t interact with julia types (I translated both to Ref{acb}). Could You give me a tip?

I don't know what an acb_ptr is defined as in Arb, so I can't answer this. But as a guess, it is interchangeable with an acb_t from the point of view of Julia.

Perhaps @fredrik-johansson can help if the hints above don't help.

[fredrik-johansson]

Yes, acb_ptr should be interchangeable with acb_t in this context (in general, the output could be a vector of length > 1, but the integration code only uses length 1).

I can't really help much with the details of the Julia/C interface, but I would certainly like to see this working!

[fredrik-johansson]

By the way, mag = Nemo.mag_struct(prec(FF), UInt(2)) can't be correct. You need to call mag_init and use some function like mag_set_d, mag_set_ui_2exp_si, arb_get_mag (and presumably set it to something like 2^-prec, not 2^prec).

[kalmarek]

@wbhart,

You can use libarb in Nemo, instead of "/usr/lib/libarb.so"

Can I? I moved acb_calc_integrate_opt into ArbTypes.jl, but then I get

ccall: could not find function acb_calc_integrate_opt_init in library libarb

(I've just updated Nemo and build it).

@fredrik-johansson thanks for noting the issue with mag; does the following fully initializes mag?

mag = Nemo.mag_struct(0, UInt(0.0))
ccall((:mag_one, "/usr/lib/libarb.so"), Void, (Ref{Nemo.mag_struct},), mag) 

(the code still segfaults)

Now it seems to me that I probably need to pass res::Ptr{Nemo.acb_struct} to :acb_calc_integrate, otherwise, how arb is going to know what to do with the memory where parent is stored (when passing just acb....?

[fredrik-johansson]

Yes, zeroing the struct is correct as initialization (though 1 then is probably not the value you want to use!).

[wbhart]

On 5 February 2018 at 16:17, kalmarek notifications@github.com wrote:

@wbhart https://github.com/wbhart,

You can use libarb in Nemo, instead of "/usr/lib/libarb.so" Can I? I moved acb_calc_integrate_opt into ArbTypes.jl, but then I get ccall: could not find function acb_calc_integrate_opt_init in library libarb (I've just updated Nemo and build it).

Yes, you certainly can. Take a look at how ccalls to Arb are made in src/arb/arb.jl

If it says it can't find the symbol, you probably need to update the sha of Arb in deps/build.jl, otherwise it checks out the wrong commit of Arb.

We ( @thofma ) will need to rebuild the Windows binary when this gets merged.

@fredrik-johansson https://github.com/fredrik-johansson thanks for noting the issue with mag; does the following fully initializes mag?

mag = Nemo.mag_struct(0, UInt(0.0))ccall((:mag_one, "/usr/lib/libarb.so"), Void, (Ref{Nemo.mag_struct},), mag)

(the code still segfaults)

Now it seems to me that I probably need to pass res::Ptr{Nemo.acb_struct} to :acb_calc_integrate, otherwise, how arb is going to know what to do with the memory where parent is stored (when passing just acb....?

You can't use C structs in Julia directly. It always creates a struct on the heap and mallocs it (the equivalent of an arb_t). Using Ref then passes that pointer to the C side.

[kalmarek]

Is there an easy way to convert between Nemo.acb_struct and Nemo.acb? I have moved everything to acb.structs and defined basic arithmetic using calls to (:libarb)[arblib.org/acb.html]: no segfaults, but I have no idea if the results are correct ;-) I tried

function (FF::AcbField)(x::Nemo.acb_struct)
    res = FF()
    ccall((:acb_set, :libarb), Void,
        (Ref{acb}, Ref{Nemo.acb_struct}), res, x)
    return res
end

function (T::Type{Nemo.acb_struct})(x::acb)
    res = Nemo.acb_struct()
    ccall((:acb_set, :libarb), Void,
        (Ref{Nemo.acb_struct}, Ref{acb}), res, x)
    return res
end

are these correct? then my code is not ;-) Note: I changed Nemo.acb_struct() to

mutable struct acb_struct
  # ...
  # all the fields

  function acb_struct()
    x = new()
    ccall((:acb_init, :libarb), Void, (Ref{Nemo.acb_struct},), x)
    return x
  end
end

[wbhart]

On 5 February 2018 at 19:40, kalmarek notifications@github.com wrote:

Is there an easy way to convert between Nemo.acb_struct and Nemo.acb?

I don't understand the difference between the two. Perhaps @thofma can help here?

I have moved everything to acb.structs and defined basic arithmetic using calls to (:libarb)[arblib.org/acb.html]: no segfaults, but I have no idea if the results are correct ;-) I tried

function (FF::AcbField)(x::Nemo.acb_struct) res = FF() ccall((:acb_set, :libarb), Void, (Ref{acb}, Ref{Nemo.acb_struct}), res, x) return resend function (T::Type{Nemo.acb_struct})(x::acb) res = Nemo.acb_struct() ccall((:acb_set, :libarb), Void, (Ref{Nemo.acb_struct}, Ref{acb}), res, x) return resend

are these correct? then my code is not ;-) Note: I changed Nemo.acb_struct() to

mutable struct acb_struct

...

all the fields

function acb_struct() x = new() ccall((:acb_init, :libarb), Void, (Ref{Nemo.acb_struct},), x) return x endend

It's ok to add another inner constructor, if that is what has happened here.

But if you add a call to init, presumably arb also requires a call to clear? If so, you need to set a finalizer in Nemo, or it will memory leak.

I actually have no idea what the current code is supposed to be doing, so I'm afraid @thofma and @fredrik-johansson are going to have to enlighten us.

[kalmarek]

@wbhart, Yes, I just added an inner constructor to allow easy creation of acb_structs from julia.

The difference between acb and acb_struct is that acb has parent fieldname and some of its fields (namely {real,imag}_mid_{d1,d2}) are of type UInt instead of Int. It seems that it (acb_struct) follows the memory layout of (acb_struct)[https://github.com/fredrik-johansson/arb/blob/master/acb.h#L29] of arb C-library exactly. Looking at the few examples of usage of Nemos acb_struct we can see that its size is used for pointer calculation and acb_vec(n) allocates vector, and returns Ptr{acb_struct}.

According to clear and finalize - well I have very little idea how interfacing arb works and who does what part of memory management - which should be pretty clear by reading my comments so far ;-)

[thofma]

I have some ideas what could be wrong, but I want to try it out first.

Is this on a branch somewhere?

[kalmarek]

Ok, I got the basic integration working here: https://github.com/kalmarek/Nemo.jl/tree/enh/acb_integrate You can do simple things such as

using Nemo
using Base.Test

const PREC = 200
const FF = AcbField(PREC)
const acb_struct = Nemo.acb_struct

@testset "integrate simple function" begin
    @test isa(Nemo.mag_struct(), Nemo.mag_struct)
    @test isa(Nemo.mag_struct(2.0^-20), Nemo.mag_struct)
    mag = Nemo.mag_struct(2.0^-prec(FF))

    function square_acbs(res::acb_struct, x::acb_struct, param, order::Int, prec::Int)
        mul!(res, x, x, prec)
        return convert(Cint, 0)::Cint
    end

    square_acbs_cf = cfunction(square_acbs, Cint,
        (Ref{acb_struct}, Ref{acb_struct}, Ptr{Void}, Int, Int))

    opts = Nemo.acb_calc_integrate_opts()
    @test isa(opts, Nemo.acb_calc_integrate_opts)

    res = acb_struct()
    status = ccall((:acb_calc_integrate, Nemo.libarb), UInt,
        (Ref{acb_struct}, #res
        Ptr{Void},        #func
        Ptr{Void},        #params
        Ref{acb_struct},  #a
        Ref{acb_struct},  #b
        Int,              #rel_goal
        Ref{Nemo.mag_struct},  #abs_tol
        Ref{Nemo.acb_calc_integrate_opts}, #opts
        Int),
        res, square_acbs_cf, C_NULL, acb_struct(0), acb_struct(2), prec(FF), mag, opts, prec(FF))

    @test status == Nemo.ARB_CALC_SUCCESS
    @test contains(real(FF(res)), 8//3)
    @test imag(FF(res)) == 0
end

I tried however to make this work with user (run-time) defined functions. According to this You can do that through a wrapper interface as closures are not (yet) c-callable. Here it is: (already included in the branch above)

struct acb_calc_func{T<:Function}
  domain::AcbField
  func::T
end

(F::acb_calc_func)(x::acb) = F.func(x)
(F::acb_calc_func)(x::acb_struct) = F(F.domain(x))

function acb_calc_func_wrap(res::acb_struct, x::acb_struct, param::Ptr{Void}, order::Int, prec::Int)
    acbFunc = unsafe_pointer_to_objref(param)::acb_calc_func
    _acb_set(res, acb_struct(acbFunc(x)))
    return Cint(0)::Cint
end

Base.cfunction(F::Nemo.acb_calc_func) =
    cfunction(acb_calc_func_wrap, Cint,
        (Ref{acb_struct}, Ref{acb_struct}, Ptr{Void}, Int, Int))

With this I can successfully evaluate the created cfunction via ccall: we evaluate acb_calc_func_wrap at x passing the actual function (acb_calc_func object) through void * param, and storing the result in res. The created cfunction follows the declaration of acb_calc_func_t (or at least I intended so):

# header as above

@testset "cfunction" begin

    f(x) = (x+1)^2

    acbFunc = Nemo.acb_calc_func(FF, f)
    @test isa(acbFunc, Nemo.acb_calc_func)

    @test acbFunc(FF(3)) == FF(f(3))
    @test acbFunc(acb_struct(FF(3))) == FF(f(3))

    let acbFunc_c = cfunction(acbFunc)
        # (3+1)^2
        res = acb_struct()
        @test ccall(acbFunc_c, Cint,
        (Ref{acb_struct}, Ref{acb_struct}, Ptr{Void}, Int, Int),
        res, acb_struct(3), pointer_from_objref(acbFunc), 1, prec(FF)) == Cint(0)
        @test FF(res) == FF(f(3))

        res = acb_struct()
        @test Nemo.acb_calc_func_wrap(res, acb_struct(3), pointer_from_objref(acbFunc), 1, prec(FF)) == Cint(0)
        @test FF(res) == FF(f(3))

        # (2im+1)^2
        res = acb_struct()
        @test ccall(acbFunc_c, Cint,
        (Ref{acb_struct}, Ref{acb_struct}, Ptr{Void}, Int, Int),
        res, 2onei(acb_struct), pointer_from_objref(acbFunc), 1, prec(FF)) == Cint(0)
        @test FF(res) == f(FF(0,2))

        res = acb_struct()
        @test Nemo.acb_calc_func_wrap(res, 2onei(acb_struct), pointer_from_objref(acbFunc), 1, prec(FF)) == Cint(0)
        @test FF(res) == f(FF(0,2))

    end
end

However I failed to pass such construct to acb_calc_integrate (the following results in a segfault)

    # header as above
    opts = Nemo.acb_calc_integrate_opts()
    mag = Nemo.mag_struct(2.0^-prec(FF))

    f(x) = (x+1)^2

    acb_f = Nemo.acb_calc_func(FF, f)
    acb_cf = cfunction(acb_f)

    res = zero(acb_struct)
    status = ccall((:acb_calc_integrate, Nemo.libarb), UInt,
            (Ref{acb_struct}, #res  -> store the result
            Ptr{Void},        #func -> cfunction(acb_calc_func)
            Ptr{Void},        #params-> pointer to actual function to evaluate
            Ref{acb_struct},  #a
            Ref{acb_struct},  #b
            Int,              #rel_goal
            Ref{Nemo.mag_struct},              #abs_tol
            Ref{Nemo.acb_calc_integrate_opts}, #opts
            Int               # prec
            ),
            res, acb_cf, pointer_from_objref(acb_f), acb_struct(-1), acb_struct(1), prec(FF), mag, opts, prec(FF))
    @show res
    @test status == Nemo.ARB_CALC_SUCCESS

    @test contains(real(FF(res)), 8//3)
    @test imag(FF(res)) == 0

Edit: when running in terminal I get

integration: GC error (probable corruption) :
Allocations: 4928131 (Pool: 4926685; Big: 1446); GC: 7
<?#0x7f4d1d74ad10::
!!! ERROR in jl_ -- ABORTING !!!

Did I botch memory management and an object is read from pointer information on the C side after being GCed?!

[fredrik-johansson]

Meanwhile here is also a WIP Sage wrapper for comparison: https://trac.sagemath.org/ticket/24686

[thofma]

OK, I got it to work. There is actually no need for acb_struct. acb is sufficient.

At the moment I don't see an easy way to get a nice API. The biggest problem is that on the C side, new elements are created, which together with prec are passed to the function we provide. We are then expected to evaluate our function with prec precision (what does this actually mean @fredrik-johansson). But this does not really fit well into our arb interface. Since we have parent objects, all our functions operate on acb elements and not on acb + prec as on the C side. So passing x -> x^2 to the integrator is quite cumbersome as seen in your example. In the end I really want an interface of the form integrate(f, a, b), where f can be anything like f = sin or f = x -> x^2.

Here is one idea to make this work. Assume the integrator wants us to evaluate an acb object z at f with precision prec. The result should be written into res. We can do this with any function f defined in Nemo by wrapping it as

(z, res, prec) -> begin
  z.parent = AcbField(prec)
  w = f(z)
  acb_set!(res, w)
end

Do you think this would work @fredrik-johansson ? Does this "evaluate f with precision prec"?

[fredrik-johansson]

Something like that; the wrapper should pass z that is a Nemo acb with parent AcbField(prec) to the Julia integrand, which can then be a normal function. It's exactly the same as the Sage code.

[kalmarek]

@thofma could You share a gist of Your changes? since I can not get it to work with or without acb_struct;

@fredrik-johansson If I understand correctly acb_calc_func_wrap (i.e. the first argument to acb_calc_integrate) will be called with arguments (res::acb_struct, x::acb_struct, param::Ptr{Void}, order::Int, prec::Int)? Suppose that I want to pass in param just a pointer to the actual function to be evaluated, then

function acb_calc_func_wrap(res::acb_struct, x::acb_struct, param::Ptr{Void}, order::Int, prec::Int)
    f = unsafe_pointer_to_objref(param)::Function
    x = AcbField(prec)(x)
    z = f(x)
    _acb_set(res, z)
    return Cint(0)::Cint
end

should do the job following @thofma remark, right?

[kalmarek]

ah finally, this works:

function acb_calc_func_wrap(res::acb_struct, x::acb, param::Ptr{Void}, order::Int, prec::Int)

    acbF = unsafe_pointer_to_objref(param)::Function
    FF = AcbField(prec)
    z = acbF(FF(x))

    _acb_set(res, acb_struct(z))
    return Cint(0)::Cint
end

acb_calc_func_wrap_c() = cfunction(acb_calc_func_wrap, Cint,
        (Ref{acb_struct}, Ref{acb}, Ptr{Void}, Int, Int))

julia> FF =AcbField();
julia> Nemo.integrate(sin, FF(0), 2FF(big(π)))
[+/- 1.11e-75] + i*0

[kalmarek]

fixed by https://github.com/Nemocas/Nemo.jl/pull/407