Enh/acb integrate

[kalmarek]

Please take a careful look, I am really not experienced in C.

With this pull You can e.g.:

const FF = AcbField()
Nemo.integrate(x -> x^2, FF(0), FF(1))
Nemo.integrate(sin, FF(0), 2const_pi(FF))

TODO:

• documentation
• tests

Note: this does not allow to pass additional parameters to the integrand. Should it?

Also: should I rewrite history (i.e. prepare a new branch with clean commits)?

[wbhart]

It looks like it tries to include a file acb_struct, which I don't see (I assume it is no longer needed).

[thofma]

Thanks. I want to play with it a bit before I merge it.

[kalmarek]

Sure, that's what I asked for.

I've already found that this is not GC-safe??:

for i in 1:10000
    @time Nemo.integrate(x -> x^2 ,CC(-1), CC(1))
end

0.003962 seconds (981 allocations: 57.355 KiB)
  0.000092 seconds (212 allocations: 8.922 KiB)
  0.000048 seconds (212 allocations: 8.922 KiB)
  0.000044 seconds (212 allocations: 8.922 KiB)
  0.000045 seconds (212 allocations: 8.922 KiB)
  0.000043 seconds (212 allocations: 8.922 KiB)
  0.000042 seconds (212 allocations: 8.922 KiB)
  0.000042 seconds (212 allocations: 8.922 KiB)

.....

signal (11): Segmentation fault
while loading no file, in expression starting on line 0
unknown function (ip: 0x7f78ba8f96e1)
unknown function (ip: 0x7f78ba8fa5f2)
unknown function (ip: 0x7f78ba8fa16e)
unknown function (ip: 0x7f78ba8fe500)
jl_gc_collect at /usr/lib/libjulia.so.0.6 (unknown line)
jl_gc_pool_alloc at /usr/lib/libjulia.so.0.6 (unknown line)
Type at /home/kalmar/.julia/v0.6/Nemo/src/arb/ArbTypes.jl:80
Type at /home/kalmar/.julia/v0.6/Nemo/src/arb/acb.jl:1630
jl_apply_generic at /usr/lib/libjulia.so.0.6 (unknown line)
acb_calc_func_wrap at /home/kalmar/.julia/v0.6/Nemo/src/arb/acb_calc.jl:7
unknown function (ip: 0x7f789467eb76)
acb_calc_integrate_gl_auto_deg at /home/kalmar/.julia/v0.6/Nemo/deps/arb/acb_calc/integrate_gl_auto_deg.c:191
acb_calc_integrate at /home/kalmar/.julia/v0.6/Nemo/deps/arb/acb_calc/integrate.c:220
#integrate#38 at /home/kalmar/.julia/v0.6/Nemo/src/arb/acb_calc.jl:28
unknown function (ip: 0x7f789467e9a6)
jl_apply_generic at /usr/lib/libjulia.

[thofma]

Yes, the handling of the pointers etc is not quite correct. It is too complicated to explain here. Do you mind if I open a new PR based on the commits here? For this it would be helpful if you could squash the commits here into one (rewrite the history).

[fredrik-johansson]

I would suggest extending the interface right away to support non-analytic/branch cut-containing integrands (we did it by passing a flag to the object function in Sage) as well as wrapping the optional settings (tolerances etc.).

[kalmarek]

ok, I tried really hard to understand the troubles with memory management and julias GC destroying objects that arb has already a pointer to, but this still eludes me.
I tried introducing

struct acb_calc_func
    f::Function
    domain::AcbField
end

function acb_calc_func_wrap(res::acb_struct, x::acb_struct, param::Ptr{Void}, order::Int, precision::Int)
    acbF = unsafe_pointer_to_objref(param)::acb_calc_func
    @assert precision == prec(acbF.domain)
    X = acb()
    X.parent = acbF.domain
    _acb_set(X, x)

    z = acbF.f(X)
    _acb_set(res, acb_struct(z))
    return Cint(0)::Cint
end

(following closely int acb_calc_func_callback of Sage, and/or this -- "Taking out the trash (or not)") but arrived at no solution.
@thofma If You could shed some light it'd be greatly appreciated!

[wbhart]

Julia can't know what Arb still has pointers to. It only knows about the objects it created and that it still has references to. Also, unsafe_blah means it is unsafe from the point of view of garbage collection. The object needs to be kept alive manually, since Julia doesn't know to keep it alive automatically. Having said that, I've no idea what the specifics of your code are here. I'll let Tommy and Fredrik respond.

[thofma]

@fredrik-johansson

        RR = RealField()
        if rel_tol is None:
            cgoal = self._prec
        else:
            tmp = <RealNumber> RR(rel_tol)
            mpfr_get_d_2exp(&cgoal, tmp.value, MPFR_RNDD)
            cgoal = -cgoal

Stupid question. Ignoring the sign change at the end, we want 2^cgoal <= rel_tol, right? After looking at the documentation of mpfr_get_d_2exp, I think we have rel_tol >= tmp * 2^cgoal >= 2^(cgoal - 1) (because 0.5 <= tmp <= 1). Is there a +/- 1 missing?

[wbhart]

Sorry there is a minor rebase of this PR that is needed after a particularly large and unwieldy patch was merged. We have now completely split Nemo.jl and AbstractAlgebra.jl and updated all the docs accordingly. Apart from minor bugs that might manifest through use, there are no similarly large merges planned in the near future (famous last words).

[thofma]

I have a local branch where I have this thing here almost working. I will finish it soon.

[thofma]

Sorry for the delay. I think it is best to wait until 0.7 is released. Then cfunction will also support closures, see JuliaLang/julia#26486. This will make everything much easier.

[vtjnash]

If the C library follows modern conventions and provides some sort of void* environment pointer, you're much better off going with that than using the closure support. Otherwise, where you once had one bug with your handling of pointers, now you'll have three.

[thofma]

Thanks for chimming in @vtjnash. Is https://julialang.org/blog/2013/05/callback still the recommended way? Is there an example somewhere in Base?

[vtjnash]

Yes

There are many examples. They aren't much changed from the blog post.

[kalmarek]

@thofma any chances of this happening soon? Could You publish Your branch so that I can play with/study Your solution?

[thofma]

My solution had some other flaws, so I took yours and fixed the memory problem. See
#407

[kalmarek]

fixed by #407