integration

[kalmarek]

Arb has acb_calc, where acb_calc_integrate is provided.
I tried to wrap the function as follows

#src/float/integrate.jl
struct acb_calc_integrate_opt_struct
    deg_limit::Cint #slong
    eval_limit::Cint #slong
    depth_limit::Cint #slong
    use_heap::Cint #int
    verbose::Cint #int

    function acb_calc_integrate_opt_struct(deg_limit::Integer, eval_limit::Integer,
      depth_limit::Integer, use_heap::Integer=0, verbose::Integer=0)
      return new(deg_limit, eval_limit, depth_limit, use_heap, verbose)
    end

    function acb_calc_integrate_opt_struct()
      opts = new()
      ccall(@libarb(acb_calc_integrate_opt_init), Cvoid,
        (Ref{acb_calc_integrate_opt_struct},), opts)
      return opts
    end
end

function acb_calc_func(out::Ptr{ArbComplex}, inp::Ptr{ArbComplex},
    param::Ptr{Cvoid}, order::Cint, prec::Cint)
    @assert iszero(order) # || isone(order) ← we'd need to verify holomorphicity
    x = unsafe_load(inp)
    f = unsafe_pointer_to_objref(param)
    ccall(@libarb(acb_set), Cvoid, (Ptr{ArbComplex}, Ref{ArbComplex}), out, f(x))
    return zero(Cint)
end

acb_calc_func_cfun() = @cfunction(acb_calc_func, Cint,
        (Ptr{ArbComplex}, Ptr{ArbComplex}, Ptr{Cvoid}, Cint, Cint))

function acb_calc_integrate(res::ArbComplex, cfun, param,
    a::ArbComplex, b::ArbComplex, rel_goal::Integer, abs_tol::Mag,
    options::acb_calc_integrate_opt_struct, prec::Integer)
    status = ccall(@libarb(acb_calc_integrate), Cint,
        (Ref{ArbComplex}, # res
        Ptr{Cvoid}, # cfun
        Any, # param
        Ref{ArbComplex}, # a
        Ref{ArbComplex}, # b
        Cint, # rel_goal
        Ref{Mag}, # abs_tol
        Ref{acb_calc_integrate_opt_struct}, # options
        Cint, # prec
        ),
        res, cfun, param, a, b, rel_goal, abs_tol, options, prec)
    return status
end

function integrate(f, a::Number, b::Number;
    rel_tol=0.0, abs_tol=precision(ArbComplex),
    opts::acb_calc_integrate_opt_struct = acb_calc_integrate_opt_struct())
    res = zero(ArbComplex)

    status = integrate!(res, f, ArbComplex(a), ArbComplex(b);
        rel_tol=rel_tol, abs_tol=abs_tol, opts=opts)

    # status:
    # ARB_CALC_SUCCESS = 0
    # ARB_CALC_NO_CONVERGENCE = 2
    if status == 2
        @warn "Arb integrate did not achived convergence, the result might be incorrect"
    end
    return res
end

function integrate!(res::ArbComplex, f, a::ArbNumber, b::ArbNumber;
    rel_tol=0.0, abs_tol=abs_tol=max(precision(a), precision(b)),
    opts::acb_calc_integrate_opt_struct=acb_calc_integrate_opt_struct())

    prec = max(workingprecision(a), workingprecision(b))

    if rel_tol <= zero(rel_tol)
        crel_goal = prec
    else
        # crel_goal = r where rel_tol ~2^-r
        crel_goal = -(ArbFloat(rel_tol).exp)
    end

    return acb_calc_integrate(res, acb_calc_func_cfun(),
        f, ArbComplex(a), ArbComplex(b),
        crel_goal, Mag(abs_tol), opts, prec)
end

but the function evaluation segfaults at line y = f(x)
modifying the function to

function acb_calc_func(out::Ptr{ArbComplex}, inp::Ptr{ArbComplex},
    param::Ptr{Cvoid}, order::Cint, prec::Cint)
    @assert iszero(order) # || isone(order) ← we'd need to verify holomorphicity
    x = unsafe_load(inp)
    f = unsafe_pointer_to_objref(param)
    println("f(2) = $(f(ArbComplex(2)))")
    println("Loaded x from $inp")
    println("x = $x")
    y = f(x)
   println("y = $y")
    ccall(@libarb(acb_set), Cvoid, (Ptr{ArbComplex}, Ref{ArbComplex}), out, y)
    return zero(Cint)
end

now segfaults at first access of the value of x, i.e. at println("x = $x") (f works correctly).

I'm quite confused, since similar code in Nemo https://github.com/wbhart/Nemo.jl/blob/5f96afc4af415df5985432b56bd1a972a0025d4b/src/arb/acb_calc.jl#L2 works as expected.

Probably crel_goal and the default values here and there need adjusting as well.

[JeffreySarnoff]

Anything not imported from the Arb C Library is likely not to work. Provision of additional functionality requires deep integration with the extant design. Just providing a julia-side struct to go along with the C constructs and using that with calls to the library certainly will segfault. While ArbNumerics does rely on some technology that Nemo has refined, it is not intended to be just like Nemo (there would be no need for both). That said, I am happy to consider your issue a request to introduce an Arb driven integration facility. It will take me some time to see what is really involved. Every structured type that does not yet exist in ArbNumerics usually involves a fair amount of work; and when there is interplay .. representing a mathematical expression to integrate and providing for Arb's integration .. it is nontrivial at best. For a glimpse, see (for example) #27 and #30

[kalmarek]

thanks for the quick answer; I'm happy to work-out a PR with this functionality.

As far as I understand ArbComplex is binary compatible with acb_t on Arb side?
Then http://arblib.org/acb_calc.html#c.acb_calc_func_t suggests that the signature for acb_calc_func should be acb_calc_func_cfun() = @cfunction(acb_calc_func, Cint, (Ptr{ArbComplex}, Ref{ArbComplex}, Ptr{Cvoid}, Cint, Cint)), however since acb_t inp is owned by the arb library I changed this the type of the second one to Ptr{ArbComplex}. Is this correct?

[JeffreySarnoff]

ArbComplex is the Julian perspective on acb_t. Each of the numeric types, ArbFloat, ArbReal, ArbComplex are implemented through mutable structs whose fields coincide with the the roles and the bitwidths and type-like intent of the C structs that the Arb C Library uses in implementing arf_t, arb_t, and acb_t values. To say the are binary compatible may be a bit too strong .. they are certainly binary conversant, though. One example is the Arb C Library matrix support and the associated specific functionality, a subset of which is available with ArbNumerics. C Matricies are stored in row major order while Julia using column major order -- so in this case the structures are not binary compatible (and a great deal of detail management is required for that facility to work).

I do not know anything about Arb's internal representation of functions or how best to provide a Julia perspective thereupon [I have not looked at that at all]. There are some very good Julia packages or parts of packages that work with math expressions, and one of them probably is appropriate for this use. Getting that and the Arb C side to just work is almost certainly going to require good planning and some quick efforts exploring efficiency and effectiveness. I suggest that before attending the Arb functions that operation on expressions (i.e. integration), first settle on how expressions are to be implemented so that they do not try to mirror how Arb does it, rather use how Arb does it (subfunctions) as scaffolding from which a clean design for memory-content mapping to/from Julia is architected.

The greatest efforts in developing ArbNumerics have been to become free of the Arb C Library manner of working and to find appropriate and much more naturally Julia ways to proceed.

The greatest time in developing ArbNumerics has been to make it so anything that works for an ArbComplex also works for ArbReal and ArbFloat and vice versas. The Arb C Library is not at all symmetric in that regard. Yet this is a most important pillar, and one of the reasons why very different sorts of expertise and need do find value in this Package.

As questions arise, do ask -- preferably with some context, and I will attempt to be of help.

[kalmarek]

By binary compatibility I meant if ArbComplex and acb_t have the same memory layout/are aligned in memory, so that if I'm given acb_t*, I can turn it into Ptr{ArbComplex} and directly unsafe_loadanArbComplex`.

---

Ok, so a brief explanation on what I'm attempting to do here:
acb_calc_integrate takes

• res::PtrToArbComplex, where it stores the result,
• a function func of signature acb_calc_func_t to be integrated,
• param parameters to be passed to func
• boundaries of integration and a bunch of non-essential parameters (not important for now)

On the julia side I create func via the appropriate call acb_calc_func = @cfunction( ... ), and pass the actual julia function f to be integrated as param. Therefore acb_calc_func(out, inp, param, ... ) works by

1. "unloading" f from param and
2. storing f(inp) inside out.

I can confirm that f is "unloaded" correctly by f = unsafe_pointer_to_objref(param).
however I'm having trouble "unloading" the argument inp: it is declared as PtrToArbComplex in the signature, so I thought that unsafe_load(inp) should do the trick, however this is the place where julia segfaults.

Is there a special procedure/function for creating an ArbComplex from acb_t* ?

[JeffreySarnoff]

Is there a special procedure/function for creating an ArbComplex from acb_t* ?
We never create a Julia-side type from an Arb-side value. We always create the Arb-side value in the process of constructing a Julia-side type. So don't think in terms of ArbComplex from an acb_t*, rather think about constructing an ArbComplex value and just know that the acb_t* becomes made properly underneath the need for us to peer at it or muck with its generation or with its removal (the construction of an Arb_ var provides for the initialization and for the finalization of the concomitant a__t* value. You cannot build an ArbComplex value from a given acb_t* value -- at least not with any assurance that the result is properly usable.

[kalmarek]

I understand, but here we have the following situation:
if a julia function is to be evaluated by the Arb-side , the former can work only ArbComplex, but the latter can supply only acb_t. So there must be a bridge in the other direction to make this work, at least when we want to call acb_calc_integrate.

What I understand now is that the segfault happens, since Ptr{ArbComplex} can not be unsafe_loaded, as ArbComplex needs the precision parameter. Similar segfault can be obtained here:

w = ArbComplex(1+im)
ptrw = Base.cconvert(Ref{ArbComplex}, w)
pptrw = Base.unsafe_convert(Ref{ArbComplex}, ptrw)
# exactly what ccall does to `w::ArbComplex{P}` when specifying `Ref{ArbComplex}` as target type
# ww = unsafe_load(pptrw) # this will sefgault

what works is

function acb_calc_func(out::Ptr{ArbComplex}, inp::Ptr{ArbComplex},
        param::Ptr{Cvoid}, order::Cint, prec::Cint)
    @assert iszero(order) # || isone(order) ← we'd need to verify holomorphicity
    x = unsafe_load(convert(Ptr{ArbComplex{Int(prec)}}, inp))
    f = unsafe_pointer_to_objref(param)
    @debug "Evaluating at" x f(x)
    ccall(@libarb(acb_set), Cvoid, (Ptr{ArbComplex}, Ref{ArbComplex}), out, f(x))
    return zero(Cint)
end

[JeffreySarnoff]

For illustrative clarification, here is the source code for evaluating an elliptic integral of the first kind (which is a function in Arb C Library that takes a complex parameter).

# void acb_elliptic_k(acb_t res, const acb_t m, slong prec)

"""
    elliptic_k(modulus)
Computes the complete elliptic integral of the first kind  K(m)
"""
function elliptic_k(modulus::ArbComplex{P}) where {P}
    isone(modulus) && return ArbComplex{P}(Inf,0)
    result = ArbComplex{P}()
    ccall(@libarb(acb_elliptic_k), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Cint), result, modulus, P)
    return result
end

Note that, at least in this sort of use, it is unnecessary to pass Ptr{ArbComplex} typed args. It works well (much better: more robust in fact) passing args of type ArbComplex{P}. Similarly, it is ill-advised to call acb_set directly because there is no assurance of good Julia-level garbage collection management [I understand that Ptr may be pointing to some well initialized and constructed and finalized variable, but there is nothing that precludes misuse of the function signature as you sketched.]

And, as a practical matter, avoiding unsafe functions unless their use is required or demanded for reasons of throughput makes maintaining and upgrading the source much easier.

[kalmarek]

I don't quite understand how this example is connected to the code for arb-based integration.

elliptic_k takes an argument, allocates return value (result), calls appropriate arb function passing the argument and the reference to result and after the ccall is finished returns the result (which was modified by the ccall).

In acb_calc_integrate we pass a pointer to arb_calc_function which is then called internally by the arb C-library, so necessarily its signature contains Ptr{ArbComplex}, I don't think there is any way around this. The actual julia function is passed as param to arb_calc_function as described above.

in acb_calc_func both of the arguments are allocated by the C library, so I think these should stay as pointers; None of them should be touched by (and in fact be visible to) julia GC. in x = unsafe_load(convert(Ptr{ArbComplex{Int(prec)}}, inp)) the data from inp is copied to julia managed variable and GC is free to collect it after the end of ccall.

I don't know how to read a piece of memory and reinterpret it as ArbComplex{P} other than unsafe_load (I'm not experienced in interconnecting C and julia). I looked thoroughly through the documentation, but couldn't find anything more appropriate. Do you have any suggestion?

[JeffreySarnoff]

Are you intending that the function to be integrated is developed as an arb_calc_function directly by the user (presumably by chaining aor nesting calls to math functions internal to the Arb C Library) rather than allowing the function to be integrated be itself a Julia function subject to constraints in typing and the use of math functions supported directly by ArbNumerics?

[JeffreySarnoff]

(continued)
The "Julia function" that becomes the integrand might be given as an expression string and converted into calls to the Arb C library where the additional speedup would be noticeable. Requiring that the user deal directly with elements of the Arb C Library is fully at odds with the reason for writing ArbFloats.jl originally and morphing it into the current ArbNumerics.jl.
--- You may be (or become) comfortable with doing that, however, apart from testing you would be writing for a very small subset of those who would really appreciate the provision of integration.

[kalmarek]

Ok, I see that I did not express the intent clearly:

• no user shall use acb_calc_func or acb_calc_func_cfun directly
• there is no string parsing/passing, julia functions are passed directly (and retrieved from param by arb_calc_func internally)

example syntax (which already works)

• integrate(sin, 0, \pi) (passing the native julia sin function)
• integrate(x->x^2-1, -1, 1) (passing an anonymous function)

[JeffreySarnoff]

oh now I see what you mean
I'll take it (as these parts) for a spin in a few different ways

[JeffreySarnoff]

please show me a code block that provides the supportive functions you wrote and include two simple examples of their use [I understand there is your PR, also having it this way simplifies my focus].

[kalmarek]

in #33 there are a couple of tests added with intended use:
https://github.com/JeffreySarnoff/ArbNumerics.jl/pull/33/files#diff-bfccaf87ae7e3d080517d8bc105803f7R13

User facing functions are:

• integrate(f, a::Number, b::Number; [rtol [, atol [, opts]]]), (https://github.com/JeffreySarnoff/ArbNumerics.jl/pull/33/files#diff-2d8d17589430d597b53a1777f1a8526aR77)
• integrate!(res::ArbComplex, f, a::Number, b::Number; [rtol [, atol [, opts]]])
  each computing the integral of f (where f is any(?) callable julia object) on a straight line between (possibly complex) numbers a and b (the latter storing the result in its first argument).

The internal (low-level) function is acb_calc_integrate (https://github.com/JeffreySarnoff/ArbNumerics.jl/pull/33/files#diff-2d8d17589430d597b53a1777f1a8526aR35) which follows exactly the signature of the arb function of the same name. acb_calc_func and acb_calc_func_cfun (latter creates C-compatible function pointer of the former) only prepare the appropriate arguments for arb's acb_calc_integrate and are intended to be never used outside of julia acb_calc_integrate. Indeed acb_calc_func_cfun() is instantiated only once in the call to abc_calc_integrate (https://github.com/JeffreySarnoff/ArbNumerics.jl/pull/33/files#diff-2d8d17589430d597b53a1777f1a8526aR114)

Julia function f (the integrand) is passed as param to julia acb_calc_integrate function; it is then cconverted to Any in the ccall (essentially becomes jl_value_t*); it resurfaces in acb_calc_func called by arb C-library (during the integration process) as the param argument and needs to be unpacked (f = unsafe_pointer_to_objref(param)) there.

I'm not sure if this explanation is what You asked for?

[JeffreySarnoff]

That is helpful.